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, .uid 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: Miha Brojan Co-Editor-in-Chief: Matevž Zupancic Section Editors: Domen Šeruga, Structural Design Matej Borovinšek, Mechanics Dominik Kozjek, Mechatronics Simon Klancnik, Production Engineering Jaka Tušek, Process Engineering Luka Lešnik, Power Engineering Joško Valentincic, Additive Manufacturing Editorial O.ce University of Ljubljana, Faculty of Mechanical Engineering SV-JME, Aškerceva 6, 1000 Ljubljana, Slovenia Phone: +386 (0)1 4771 137 info@sv-jme.eu, http://www.sv-jme.eu Technical Editor: Pika Škraba Print: Gra.ka Gracer d.o.o. printed in 190 copies President of Publishing Council Mihael Sekavcnik University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Vice-President of Publishing Council Matej Vesenjak University of Maribor, Faculty of Mechanical Engineering, Slovenia 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 Founding Editor Bojan Kraut University of Ljubljana, Faculty of Mechanical Engineering, Slovenia ISSN 0039-2480, ISSN 2536-2948 (online) © 2025 with Authors. General information Strojniški vestnik – Journal of Mechanical Engineering is published in 6 double issues per year. International Editorial Board Ha.z Muhammad Ali, King Fahd U. of Petroleum & Minerals, Saudi Arabia Josep M. Bergada, Politechnical University of Catalonia, Spain Anton Bergant, Litostroj Power, Slovenia Matej Borovinšek, University of Maribor, Slovenia Filippo Cianetti, University of Perugia, Italy Peng Cheng, Virginia State University, USA Franco Concli, University of Bolzano, Italy J.Paulo Davim, University of Aveiro, Portugal Igor Emri, University of Ljubljana, Slovenia Imre Felde, Obuda University, Faculty of Informatics, Hungary Soichi Ibaraki, Kyoto University, Department of Micro Engineering, Japan Julius Kaplunov, Brunel University, West London, UK Iyas Khader, Fraunhofer Institute for Mechanics of Materials, Germany Simon Klancnik, University of Maribor, Slovenia Jernej Klemenc, University of Ljubljana, Slovenia Milan Kljajin, J.J. Strossmayer University of Osijek, Croatia Dominik Kozjek, University of Ljubljana, Slovenia Peter Krajnik, Chalmers University of Technology, Sweden Janez Kušar, University of Ljubljana, Slovenia Luka Lešnik, University of Maribor, Slovenia Edgar Lopez, University of Istmo, Mexico Trung-Thanh Nguyen, Le Quy Don Technical University, Vietnam Vladimir Popovic, University of Belgrade, Serbia Franci Pušavec, University of Ljubljana, Slovenia Mohammad Reza Safaei, Florida International University, USA Silvio Simani, University of Ferrara, Italy Marco Sortino, University of Udine, ItalyDomen Šeruga, University of Ljubljana, Slovenia Jaka Tušek, University of Ljubljana, Slovenia Branko Vasic, University of Belgrade, Serbia Arkady Voloshin, Lehigh University, Bethlehem, USA Institutional prices include print & online access: institutional subscription price and foreign subscription €100,00 (the price of a single issue is €20,00); general public subscription and student subscription €50,00 (the price of a single issue is €10,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. Every manuscript submitted to the SV-JME undergoes a peer-review process. We would like to thank the reviewers who have taken part in the peer-review process. 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 and Inovation Agency. Strojniški vestnik - Journal of Mechanical Engineering is available on https://www.sv-jme.eu. Contents Strojniški vestnik - Journal of Mechanical Engineering Volume 71, (2025), Number 5-6 Ljubljana, May-June 2025 ISSN 0039-2480 Published every two months 149 Analytical, Numerical 1D and 3D Water Hammer 192 Tooth Contact Analysis of Involute Beveloid Gear Based on Investigations in a Simple Pipeline Apparatus Higher-Order Curve Axial Modification Anton Bergant, Zlatko Rek, Kamil Urbanowicz Yongping Liu, Qi Chen, Changbin Dong 157 Study on the Influence of the Splitter Blade Length on 199 Geometric Design Method of Lightweight Line Gear Radial and Axial Force of a Centrifugal Pump Mechanism Qijiang Ma, Zhenbo Liu, Sen Jiang Chao He, Yangzhi Chen, Xiaoxiao Ping, Zhen Chen, Maoxi Zheng, Qinsong Zhang 169 Multi-Objective Optimization Design of the Ejector Plate for Rear-Loader Garbage Trucks 207 Simulation Analysis and Experimental Study on Vibration Fu-sheng Ding, Hong-ming Lyu, Jun Chen, Hao-ran-Cao, Reduction Performance of Groove-Textured Friction Pair Lan-xiang Zhang Surfaces Wusheng Tang, Yufei Nie, Zhuo Zhang, Wei Lin, YanKai Rong, Yaochen Shi, Ning Ding 179 Identification Method of Tire-Road Adhesion Coefficient Based on Tire Physical Model and Strain Signal for Pure Longitudinal Slip Jintao Zhang, Zhecheng Jing, Haichao Zhou, Yu Zhang, Guolin Wang ON THE COVER The leading paper Analytical, Numerical 1D and 3D Water Hammer Investigations in a Simple Pipeline Apparatus authored by Anton Bergant, Zlatko Rek, and Kamil Urbanowicz, presents advanced analytical and numerical methods for simulating water hammer phenomena. It compares 1D and 3D simulations to experimental data, highlighting the accuracy of computational models in predicting complex hydraulic behaviors essential for engineering pipelines and systems. The cover graphic shows a dynamic water splash against a clear blue sky, symbolizing rapid pressure fluctuations studied in water hammer research. Image Courtesy: iStock.com SV-JME . VOL 71 . NO 1-2 . Y 2025 . 147 Analytical, Numerical 1D and 3D Water Hammer Investigations in a Simple Pipeline Apparatus Anton Bergant1,2 – Zlatko Rek2 – Kamil Urbanowicz3 1 Litostroj Power d.o.o., Slovenia 2 University of Ljubljana, Faculty of Mechanical Engineering, Slovenia 3 West Pomeranian University of Technology in Szczecin, Faculty of Mechanical Engineering and Mechatronics, Poland anton.bergant@litostrojpower.eu Abstract This paper deals with the analytical and numerical simulation of a water hammer in the reservoir-pipeline-valve (RPV) system. An analytical solution of the water hammer equations with unsteady friction term was derived for transient laminar pipe flow in a RPV system. For simulation of an arbitrary flow situation a number of one-dimensional (1D) numerical methods have been developed. The physically based method of characteristics proved to be computationally efficient and can handle complex boundary conditions. The accuracy of the 1D numerical model is increased by introducing terms that take into account 3D effects (example of an unsteady skin friction). The 3D model predicts these influences directly and represents an excellent tool for researching multidimensional properties of fluids (numerical laboratory). Calculation results based on 1D and 3D numerical models are in good agreement with results of measurements taking into account adequate prediction and modelling of influential physical parameters during laminar and low-Reynolds number turbulent water hammer events. Quantitative comparison analysis yields up to 2 % difference in maximum head at the valve and up to 5 % relative difference in pressure head drop at the midpoint of the pipeline monitored over the first four positive pressure pulses. Keywords pipeline, water hammer, analytical solution, method of characteristics, computational fluid dynamics, unsteady skin friction Highlights . Convolution-based unsteady friction model accurately captures water hammer wave mechanics. . Analytical 1D laminar model effectively applies to low-Re number turbulent water hammer. . Indirect 1D and direct 3D friction approaches yield similar responses to water hammer disturbances. . 1D and 3D models successfully validated against rapid water hammer experiments. 1 INTRODUCTION Hydraulic piping systems, such as those found in hydroelectric power plants, water supply networks, operate under different flow conditions. A change in flow rate causes a pressure rise and drop in the system. Hydraulic transients in piping systems (water hammer, hydraulic vibrations) can induce extreme pressures, formation of large local vapour cavities and distributed cavitating flow zones (zones of negative pressure), liquid and structural vibrations, and mass oscillations [1-3]. Transient loads are kept within the permissible limits by water hammer control means that include (i) alteration of operational regimes (valve closure time, limitation of hydraulic turbine maximum output, decrease of maximum flow rate), (ii) installation of surge control devices in the system (additional flywheel, air vessel, surge tank, pressure regulating valve, air valve) and (iii) redesign of the flow-passage system layout (pipeline diameter, pipe-wall material, change in elevation) [4-6]. In this paper, we investigate a water hammer, which is induced by an aperiodic change in the pipe flow velocity. The theoretical analysis of pressure changes in systems by solving the classic 1D water hammer equations gives reliable results as long as the conditions in the derivation of these equations are valid. In practice, the conditions in piping system may be far from the idealized situation described by the classic water hammer equations. For example, the steady skin friction model does not produce sufficient damping for rapid transient events [7,8]. Consequently, we use more advanced models of unsteady skin friction, which include the influence of two-dimensional (2D) and 3D effects [9-12]. In engineering practice, many other discrepancies can occur, such as: air (free and dissolved) in the liquid, transient cavitation, fluid-structure interaction (particularly if the pipes are not rigidly fixed or the excitation is severe), viscoelastic behavior of the pipe wall (in polymer pipes or if steel pipes deform plastically), as well as leaks and blockages in the pipeline [13]. In this paper, we will discuss the influence of unsteady skin friction indirectly with the help of a frozen-viscosity convolutional 1D model and directly with 3D numerical calculations [14]. It is desirable to solve 1D water hammer equations analytically. Unfortunately an exact solution for the water hammer equations with consideration of unsteady skin friction exists only for transient laminar pipe flow situation [15]. Therefore, a numerical 1D model is used in industry for treatment of both transient laminar and turbulent pipe flows. In our paper the numerical 1D model is based on the method of characteristics (MOC), which is most suitable for solving water hammer equations [1]. However, we will test the analytical solution for transient laminar pipe flow for low-Reynolds number transient turbulent flow situation. Computational fluid dynamics (CFD) is today, with increasingly powerful computers, a very useful tool in the analysis of transient flows and phenomena in hydraulic pipelines and devices [16,17]. So, for example, in [18], Saemi et al. studied the phenomenon of water hammer and analyzed the influence of the outlet boundary condition (flow reduction curves and valve modelling). They also analyzed the 3D effects caused by valve closure and the turbulence structure during the occurrence of water hammer. In [19], Cao et al. presented an alternative 3D CFD model for the systematic investigation of the dynamic characteristics of flow under transient conditions in a hydraulic pipe. The characteristics of the ball valve under static and dynamic conditions were analyzed in detail and the effect of closing time on the head loss coefficient and the discharge coefficient was carefully investigated. In addition, coupled 1D and 3D method can be used for flow situations where 3D effects prevail in a small part of the hydraulic pipeline system [20]. The objective of this paper is to extend the application of the analytical model developed for laminar water hammer [15] to low-Reynolds number turbulent water hammer. The analytical model will be presented in a novel form with defined optimized upper limit of the nominally infinite sum in the solution. Similarly, the effectiveness of Zielke’s weighting function [21] in 1D unsteady friction models will be validated for low-Reynolds number turbulent transient pipe flow. In our 3D investigations the CFD package ANSYS Fluent Release 18.2 [22] will be used. A particular emphasis will be given to selection of the most appropriate turbulence model for a turbulent water hammer. The presented 1D and 3D models will be validated against the results of measurements in a simple pipeline apparatus [8]. Theoretical and experimental investigations in this paper should deepen the understanding of water hammer phenomena in liquid-filled piping systems. 2 METHODS & MATERIALS Water hammer describes the propagation of pressure waves in pipes with liquid. In most engineering applications the cross-sectional dimension of the pipe is negligible compared to the length of the pipe. Assuming uniform flow and neglecting small convective terms yields the following two 1D unsteady pipe flow equations, (i) the continuity equation and (ii) the equation of motion [1,2] ˜Ha2 ˜Q °. 0, (1) ˜t gA ˜x ˜H 1 ˜Q QQ °° f 2 . 0, (2) ˜x gA ˜t 2gDA in which H is piezometric head (head), t time, Q discharge, x distance along the pipe, a pressure wave speed, g gravitational acceleration, A pipe cross-sectional area, f Darcy-Weisbach friction factor and D inner pipe diameter. The symbols are explained as they first appear in the paper. In Equation (2), we traditionally use the steady-state Darcy-Weisbach friction coefficient f. In the case of fast transient phenomena, the steady friction coefficient is modified by introducing a non-stationary friction term. The friction coefficient f is expressed as the sum of a steady (or quasi-steady) part fq and the unsteady part fu. In this paper, we will use the frozen- viscosity convolution-based unsteady friction model [21,23,24] in which the unsteady part is defined by the convolution of a weighting function W0 with past temporal accelerations 32. At°Q * f˜ *Wt td.. . t, (3) ˆ *0 u DQQ 0 °t in which . is kinematic viscosity and W0 weighting function. The unsteady part is a result of inverse Laplace transformation of the quasi-2D water hammer Laplacian-solution in the frequency domain. Weighting functions have been developed for transient laminar [21] and turbulent pipe flow [25,26]. The frozen-viscosity convolution-based unsteady friction model is suitable for water hammer flows where pressure waves travel at sound propagation speed and decelerate/accelerate the flow in short time periods [27]. The importance of unsteady friction in water hammer flows is decreased with increasing pre-transient Reynolds number [28]. The classical 1D water hammer equations will be solved by (i) the analytical method and (ii) the MOC. Finally, (iii) the 3D numerical solution of the Navier-Stokes equations will be presented and applied to water hammer problems. 2.1 Analytical Solution of Water Hammer Equations Analytical solutions of water hammer Eqs. (1) and (2) for frictionless pipe systems and for systems with consideration of quasi-steady friction have been developed by Rich in the middle of the previous century [29]. The time-domain solutions have been obtained by inverse Laplace transforms. Analytical solutions of water hammer equations by consideration of an unsteady state friction term are a difficult task [15]. The time-domain solution for transient laminar pipe flow dependent on the dimensionless water hammer number Wh, defined below Eq. (4), has been recently developed by Urbanowicz et al. [30]. Exact formulas for pressure (pressure head), flow velocity and wall shear stress have been applied for the case of instantaneous valve closure in a reservoir-horizontal pipe-valve system. In this paper we will validate pressure head solution for the transient laminar and low-Reynolds number turbulent cases in a nearly horizontal pipe. The novel analytical form for pressure head in transient laminar pipe flow is [30] 2 . kl°.2 kl, Wh . L.. hxt ˜hx(,)° Q0ˆ() e fWh kl. ,HˆTkl, ,(4) (,) 0 B .1 k˜1.l˜1 . a.. in which .k,l =2(k + l -2) +(x/L)·(-1)k+1 is index, Tk,l =(a/L)t -.k,l time-index dependent function, L pipe length, h(x,0)= h(0,0) +8WhBV0A(x/L) pressure head, B = a/(gA) characteristic impedance, V0 initial average flow velocity, Wh =(Ma0/Re0)(4L/D) water hammer number, Ma0= V0/a initial Mach number, Re0= V0D/. initial Reynolds number, h(0,0) pressure head at the valve (h = H for horizontal pipe), fWh-k,l time-index-water hammer function, and H Heaviside step function. The time-index-water-hammer function is defined by the following expression 2 ˆWh TW , h 4 kl, ˜ kl Tkl, f ˜, °.2 ˜ˆWh. e . Whkl kl, . 2 . ˆWh .ˆ Wh k,lk,l ..1 WhTkl, ˜ˆkl, . erfc . (5) 2 2 T    kl,   The infinite sum in the analytical solution, Eq. (4), does not mean that many terms need to be taken into account to define properly the calculated results. The upper limit of this sum determines the number of amplitudes that will be correctly simulated. In our case study (Section 3.1), during the first four theoretical periods of the pipeline (4Tp; Tp =4L/a), the pressure head history can be properly calculated by just eight terms (l =8). The small number of terms taken into account does not have an impact on the fit of the simulated results, so the proposed analytical solution is free of the typical Gibbs phenomenon present in other series solutions. Zielke’s analytical unsteady friction formulation [21] has been previously validated for low Re-number turbulent pipe flow numerical solutions by Trikha [31], Bergant et al. [8] and Martins et al. [12]. An analytical solution for the transient turbulent pipe flow with consideration of unsteady skin friction is not available yet. It should be noted that semi-analytical methods for transient turbulent pipe flow including unsteady friction term have been developed by Hullender [32], and García García and Alvari [33]. 2.2 1D Method of Characteristics Water Hammer Model The MOC transformation of Eqs. (1) and (2) produces the water-hammer compatibility equations which are valid along the characteristic lines. The compatibility equations in the finite-difference form are numerically stable unless the friction is dominant and the computational grid is coarse and, when written for computational section i, are [1] • along the C+ characteristic line (.x/.t=a) H ˜H °aQ ˜Q , ˜˜ ..... . it i 1,t .t gA , ˜˜ uit di 1,t .t fx . ° Q Qd ˆ0, (6) .... 2 u it, i˜˜,t .t 1 2gDA • along the C– characteristic line (.x/.t=–a) H ˜H ˜aQ ˜Q it, i.˜1,t .t gA °° d.it, °.ui.˜,t .t. 1 fx . ˜°Q. ˆ0, (7) °. 2 d it, °Qu i.˜1,t .t 2gDA in which Dx is space step and Dt time step. The discharge at the upstream side of the computational section ((Qu)i) and the discharge at the downstream side of the section ((Qd)i) are identical for the pure water-hammer case where pressure remains above the liquid vapour pressure (transient liquid pipe flow). At a boundary (reservoir, valve, pump, turbine), a device-specific equation replaces one of the water-hammer compatibility equations. 2.3 3D Computational Fluid Dynamics Water Hammer Model 3D CFD water hammer is simulated numerically as a time-dependent flow of a viscous compressible liquid under isothermal conditions. Flow is governed by the conservation laws for (i) mass ˜ˆ ˜ ° .ˆui .. 0, (8) ˜t ˜xi and (ii) momentum. The corresponding equations are presented in tensor form with Einstein’s summation convention [34]. In the case of transient laminar flow, the momentum conservation equation is ˜˜ ˜p (.ui) ° (.uuij) .. ° ˜t ˜x ˜x ji  ˜ .˜ui ˜uj 2 .˜ul  °ˆ° . ij  , (9) ˜xj  ˜xj ˜xi 3 ˜xl   in which ui is flow velocity, p pressure, . density and µ dynamic viscosity. The Kronecker delta takes the values dij = 0 for i .j and dij = 1 for i = j (i, j, l =1, 2,3). In the case of transient turbulent flow, we use the model of Reynolds-averaged unsteady Navier–Stokes equations (URANS). The momentum conservation equation in this case is ˜˜ ˜p °.ui.. °.uui j..ˆ . ˜t ˜x ˜x ji ˜ ˜ui ˜uj 2 ˜ul ˜ . . . ˆ . ˆ. ii., (10) ij ° uu'' ˜x ˜x ˜x 3 ˜x ˜x j   ji l  j in which ui is time averaged velocity and ui' velocity fluctuations. The Reynolds stresses ˜°uu'' are modeled with an appropriate ii turbulence model. In the case of RANS turbulence models, we used the Boussinesq hypothesis for Reynolds stresses ..ui .uj . 2..ul . ˜uu''°. k. , (11) ii t.. ˜. t  .x .x 3 .x ˆ ji ˆ l  in which µt represents the turbulent viscosity and can be expressed as a function of k and e, or k and .. Additional turbulence transport equations can include the following quantities: turbulent kinetic energy (k), dissipation rate of turbulent kinetic energy (e) or specific dissipation rate of turbulence (.). As the name suggests, Transition SST model predicts laminar-turbulent transition flow regime. Details of the mathematical modelling of turbulence are given in [22,34,35]. In next two sub-sections, we present (i) the valve closure model and, (ii) the mesh model and simulation parameters, two important issues in 3D CFD water hammer modelling. 2.3.1 Valve Closure Model The water hammer is triggered by a sudden stoppage of the steady flow in the pipe by the downstream end valve. In the numerical simulation, this can be implemented as instantaneous closure, which means that at the closure start time t0 the outlet boundary condition is changed to a wall boundary condition. This approach is appropriate for nearly instantaneous flow stoppage (tc << 2L/a) with nearly vertical pressure wave fronts. This does not correspond to the flow conditions, where the flow at the valve is continuously reduced to zero during the longer closing time tc. Here, we constructed a function that prescribes the time dependence of the flow at the valve position. Thus, the axial velocity vx(t) varies continuously from a fixed average value V0 to a value of zero as follows ˆ V : 0 0 tt vt..V1 .cos .0/ °° : tt .tt0 c x˜°.0 ˜ ˜˜tt°tc 0, (12) . : tt t . 00 c  in which t0= 0.2 s and tc = 0.009 s. Other authors have used a similar power function approach for the valve closure [36]. 2.3.2 Mesh Model and Simulation The computational area is represented by a straight round tube of length L = 37.23 m and inner diameter D =22.1 mm. The domain is discretized with a hexahedral mesh with cell size .c ˜2 mm, Fig. 1. The number of cells is Nc =3788484. Martins et al. in [37] suggested ˜150000 as an optimal computational mesh for a geometrically (cross-sectional area of the pipe) and hydrodynamically similar case to ours cells for the laminar regime and ˜270000 cells for the turbulent regime. Since the two numerical models are comparable, we did not perform an analysis of the grid dependence of the results. Fig. 1. Pipeline 3D computational grid Calculations were performed with the following data from the experiment: water temperature Tw = 15.4 °C, reference pressure p0=313.92 kPa, water density at reference pressure .0 = 999.14 kg/m3 and pressure wave speed a = 1319 m/s. The dynamic viscosity of water at this temperature is µ=1.12546×10–3 Pa·s. Because of pressure propagating waves that occur during the water hammer event, we should treat water as a compressible liquid. A compressible liquid model was used, which establishes a non-linear relationship between density and pressure under isothermal conditions. The model is derived from the simplified Tait equation of state [38] ˜.n 1  np p0 ° . .. , (13).˜0 ˆ K0 in which K0= .0 a2=1.7383×109 Pa is a bulk modulus at reference pressure and n =7.15 is density exponent. For the numerical solution, we used a “pressure-based” transient type solver with either laminar or turbulent viscosity model. The turbulent flow regime was treated with the following three turbulence models: (i) Realizable k–e model with standard wall functions, (ii) Shear–Stress Transport (SST) k–. model and (iii) Transition SST model [22]. The boundary conditions were as follows: (i) at the inlet (x =0) the total pressure p0 was prescribed, (ii) at the outlet (x = L) the valve closing function vx(t) was defined by Eq. (12) and (iii) at the pipe wall the velocity was v=0. The initial condition for the pressure was p = p0. The initial velocity values in the laminar regime were prescribed as v=(2V0 [1–4r²/D²], 0) and in the turbulent regime as v=(V0, 0); r is pipe radius. In the solver, the SIMPLE method was selected for velocity– pressure coupling. Spatial discretization was: (i) “least squares cell based” for gradient, (ii) “second order” for pressure and (iii) “second order upwind” for density, velocity and turbulent quantities. The “first order implicit” formulation was used for the transient scheme. The under-relaxation factors were: 0.3 for pressure, 1 for density and turbulent viscosity, 0.7 for velocity and 0.8 for turbulence quantities. In the model, we selected two pipeline sections where we saved the pressure value in each time step, namely at the midpoint of the pipe at x =18.615 m (index mp in the graphs) and at the downstream end valve at x = 37.23 m (index v in graphs). The numerical simulation of the water hammer event took place in three steps: (i) First, a calculation was made for steady-state conditions to develop a steady-state pipe flow velocity profile. (ii) A switch to transient calculation mode was made with a time step.t =10-3 s for a duration of t =1 s. (iii) This was followed by resetting the time to zero, turning on the pressure output at the observed locations and calculating until time t0 =0.196 s, when the time step was reduced to .t =10-5 s. At time t = t0, the valve started to close. At time t = t0+ tc, the valve was completely closed. After the valve was closed, the calculation took another 0.4 s, which means that three full pressure fluctuations between the maximum and minimum values could be observed during this time. During this time period, the water flow was four times in the positive direction and three times in the negative direction. Such a small-time step was necessary, since the pressure rise in rapid downstream-end valve closure event is a very fast phenomenon. 2.4 Experimental Pipeline Apparatus A versatile pipeline apparatus for investigating water hammer and column separation events was constructed at the University of Adelaide, Australia [39]. The apparatus comprises a straight 37.23 m Pressurized reservoir T1 (Ux = ±0.01 m) long sloping copper pipe of 22.1 mm (Ux = ±0.1 mm) internal diameter and of 1.63 mm (Ux = ±0.05 mm) wall thickness, see Fig. 2. The uncertainty in the measurement Ux is expressed as a root-sum-square combination of bias and precision error [40]. A water hammer event in the pipeline was induced by rapidly closing the downstream-end ball valve. The reference [39] describes the experimental pipeline apparatus in detail. 3 RESULTS AND DISCUSSION The results of calculation with 1D analytical and numerical MOC models, and 3D numerical model (3D CFD) are validated against the results of measurements with the Adelaide pipeline apparatus (Fig. 2). We compare the computational and experimental results for the rapid closing of the valve installed at the downstream end of the pipe for two typical cases of water hammer with initial flow velocities V0={0.1;0.3} m/s (initial laminar (Re0=1960) and turbulent pipe flow with low Reynolds number (Re0=5880), respectively) [8]. The static head in the upstream end pressurized tank HT,2=32 m is the same for both cases. The measured wave speed of the propagation of the first steep pressure head rise is a=1319 m/s. The pressure head wave speed is slightly decreasing during the transient event because the liquid has an extra inertia due to the velocity distribution, which is related to the momentum correction factor [41]. The results are compared at the valve (Hv) and at the midpoint of the pipeline (Hmp) (Fig. 2). The goal of the analysis is the validation of the two 1D models by consideration of quasi-2D frozen-viscosity convolutional model of unsteady wall friction, and further to extract the differences in the way of treating the transient flow according to the 1D and 3D methods. 3.1 Comparison of 1D Analytical and Experimental Results 1D analytical calculations were performed for instantaneous valve closure. The results for the run with initial flow velocity V0= 0.1 m/s are compared in Figs. 3a and b. The computational results obtained with the 1D analytical model (Eqs. (4) and (5)) developed for transient laminar pipe flow [30] agree well with the measured results in terms of attenuation, shape and timing of the pressure pulses. The maximum analytically calculated bulk head at the valve Hv–max–A = 45.7 m differs for 1 % from the measured one Hv–max–E = 45.8 m. The head trace at the midpoint has been selected for the assessment of the attenuation rate over the first four positive bulk pressure head pulses (time span of 12L/a). The percentage of the pressure head drop Dhmp1–4 is calculated by the following equation h ° h mp1 mp4 14 (14) ˜hmp ..100 ° , hmp1 L = 37.23 m Fig. 2. Adelaide experimental pipeline apparatus layout 60.0 60.0 50.0 50.0 H (m) H (m) 40.0 30.0 40.0 30.0 mp v 20.0 20.0 10.0 10.0 a) Time (s) b) Time (s) 100.0 100.0 80.0 60.0 80.0 60.0 H (m) H (m) 40.0 40.0 mp v 20.0 20.0 0.0 0.0 -20.0 -20.0 Time (s) Time (s) c) d) Fig. 3. Comparison of measured and 1D analytical computed head histories for laminar (Re0 = 1960) and for turbulent pipe flow case (Re0 = 5880); a) and c) at the valve; b) and d) at the midpoint in which hmpi is pressure head at the midpoint of the pipeline at elevation of 1.02 m (i=1,2,3,4), see Fig. 2. The calculated percentage of the pressure head Dhmp1–4–A =3.8 % is slightly lower than measured one Dhmp1–4–E =5.4 % (1.6 % difference). The actual effective valve closure time is about 0.004 s, well below the wave reflection time 2L/a =2·37.23/1319=0.0565 s. In our case the assumption of the instantaneous valve closure is acceptable as can be seen from the comparison between the measured first pressure head rise and the theoretical one in Fig. 3a. The results for the test with initial pipe flow velocity V0= 0.3 m/s are compared in Figs. 3c and d. The unsteady friction part based on the transient laminar flow does contribute to additional energy loss for the low Reynolds number turbulent flow (Re0< 104) [31]. The maximum analytically calculated bulk head at the valve Hv–max–A = 73.0 m is slightly higher than the experimentally predicted head Hv–max–E = 71.9 m (1.7 % difference). The calculated percentage of the pressure head drop over the observed time span Dhmp1–4–A =7.5 % is slightly lower than measured one Dhmp1–4–E =9.8 % (2.3 % difference). The magnitude of the deviations between calculated and measured results for turbulent flow with low Reynolds number (Figs. 3c and d) is of the same order as the one for laminar flow case (Figs. 3a and b) Finally, it should be noted that in pipeline engineering practice, the acceptable difference bewteen the calculated and measured results is 5 % for the maximum water hammer head and 10 % for the estimated friction factor. 3.2 Comparison of 1D MOC Numerical and Experimental Results 1D MOC calculations were performed for actual valve closing time tc =0.009 s and the number of pipe sections in numerical model N=16. The results for the laminar flow case are compared against the measured results in Figs. 4a and b. The computational results obtained with the frozen-viscosity convolutional model of unsteady friction [21] agree well with the measurement results in terms of attenuation, shape and timing of the pressure pulses. The maximum MOC calculated bulk head at the valve Hv–max–MOC =45.8 m matches the measured head Hv–max–E =45.8 m. There is a slight difference between the calculated Dhmp1–4–MOC =4.3 % and the measured percentage of the pressure head drop at the midpoint of the pipeline Dhmp1–4–E =5.4 % (1.1 % difference). The degree of discrepancies between the 1D MOC and the 1D analytically calculated results when compared to the measured results is small. There are slight differences between the shape and timing of the pressure histories (compare Figs. 3a and b, and Figs. 4a and b). The results for the test with initial pipe flow velocity V0=0.3 m/s are compared in Figs. 4c and d. The maximum MOC calculated bulk head at the valve Hv–max–MOC =73.1 m is slightly higher than the measured head Hv–max–E =71.9 m (1.7 % difference). There is an excellent match between the calculated percentage of the pressure head drop Dhmp1–4–MOC =9.9 % and the measured one Dhmp1–4–E =9.8 %. The inclusion of weights for past velocity changes developed in the frozen-viscosity convolutional model for transient laminar pipe flow [21] contributes to additional energy loss for transient turbulent flow with low Reynolds number (Re0<104) [31]. The magnitude of the deviations between calculated and measured results for turbulent flow with low Reynolds number (Figs. 4c and d) is similar to that for laminar flow case (Figs. 4a and b). 3.3 Comparison of 3D CFD Numerical and Experimental Results Comparisons of 3D numerical and experimental results for transient laminar and low Reynolds number turbulent flow cases are presented in this section. 3.3.1 Laminar Flow Case First, we carried out a numerical simulation of water hammer for the case of laminar pipe flow at initial Reynolds number Re0=1960. The average velocity of the steady water flow in the pipe was V0=0.1 m/s. The parallel calculation took 5.67 days on 23 cores of a Supermicro workstation with four AMD Opteron® 8439 SE processors at a frequency of 2.8 GHz. It should be noted that the calculation with the 1D models only takes a few minutes. Figure 5 compares head histories at the valve and at the midpoint. We may see a good agreement between experiment and numerical calculation. The maximum numerically predicted head at the valve Hv–max–CFDlam =45.2 m is slightly lower than the experimentally predicted one Hv–max–E =45.8 m (1.3 % difference). The calculated 60.0 60.0 50.0 50.0 40.0 (m) H (m) 40.0 H mp 30.0 v 30.0 20.0 20.0 10.0 10.0 a) Time (s) b) Time (s) 100.0 100.0 80.0 60.0 80.0 60.0 (m) H (m) v 40.0 40.0 H mp 20.0 20.0 0.0 0.0 -20.0 -20.0 Time (s) Time (s) c) d) Fig. 4. Comparison of measured and 1D MOC computed head histories for laminar (Re0 = 1960) and for turbulent pipe flow case (Re0 = 5880): a), c) at the valve; b), d) at the midpoint 60.0 60.0 50.0 50.0 (m) 40.0 H mp 30.0 H (m) v 40.0 30.0 20.0 20.0 10.0 10.0 a) Time (s) b) Time (s) Fig. 5. Comparison of measured and computed 3D CFD head histories for laminar pipe flow case (Re0 = 1960): a) at the valve; b) at the midpoint percentage of the pressure head drop over the observed time span Dhmp1–4–CFDlam =2.0 % is lower than the measured one Dhmp1–4–E =5.4 % (3.4 % difference). The 3D CFD compressible liquid model, which establishes a non­linear relationship between density and pressure under isothermal conditions, accurately predicts the water hammer wave mechanics. The head wave history is repeated periodically with a frequency of 8.7 Hz, only the amplitude of the head rise decreases and attenuates during the transient event [8]. It should be pointed out that the damping time of the pressure fluctuation depends on the length of the pipe L at an assumed constant pipe diameter D; the unsteady friction part against the quasi-steady state friction part is decreasing as the pipe length is increasing. In case of unsteady friction dominated pipe flow situations for dimensionless quantities of L/D, quasi-steady friction factor fq0, Ma0 and Re0 one can calculate the steady part (fq0Ma0L/(2D)) and the unsteady part ((2LMa0/(DRe0))0.5) of the friction damping coefficient. For our laminar flow case ( fq0=0.035) are 0.002 and 0.011, respectively, see details in Duan et al. [42]. 3.3.2 Low-Reynolds Number Turbulent Flow Case 3D CFD numerical simulation of water hammer for the case of turbulent flow at Re0=5880 was carried out next. The average velocity of the initial pipe flow velocity was V0=0.3 m/s. The following three commonly used turbulence models were selected [22,34,35]: (i) the two-equation k–e and (ii) k–. models, and (iii) the four-equation Transition SST model. Parallel calculations on the same workstation as for laminar flow case took: 6.86 days for the k–e model, 3.71 days for the k–. model, and 5.53 days for the Transition SST model. Figure 6 compares head histories at the valve and at midpoint for all three turbulence models. Even in this case, the agreement between experiment and numerical calculations is reasonable. The maximum 3D CFD calculated bulk head at the valve Hv–max–CFDturb ={turbulence model k–e: 71.5 m; k–.: 72.6 m; Transition SST: 71.7 m} agree well with the experimentally predicted head Hv–max–E =71.9 m {difference for k–e: 0.6 %; k–.: 1.0 %; Transition SST: 0.3 %}. The calculated percentage of the pressure head drop for the three turbulence models Dhmp1–4–CFDturb ={k–e: 4.1 %; k–.: 5.8 %; Transition SST: 5.1 %} is lower than the measured one Dhmp1-4-E = 9.8% {difference for k–e: 5.7 %; k–.: 4 %; Transition SST: 4.7 %}. The Transition SST model approximates the experiment slightly better (Figs. 6e and f) than the two-equation models (Figs. 6a to d). In the first period, the calculated head differs very little from the measured one, and in the following periods the discrepancy slowly grows. The amplitude of the pressure rise does not decrease as fast as in the experiment, and the frequency of the pressure head fluctuation also decreases compared to the measured values. We may conclude that small discrepancies occur Fig. 6. Comparison of measured and computed 3D CFD head histories for turbulent pipe flow case (Re0 = 5880): a), c), e) at the valve, and b), d), f) at the midpoint because the Re number is very low; therefore, there is no turbulence model, which would have a distinct advantage. 4 CONCLUSIONS It has been shown that both 1D analytical and 1D MOC models with consideration of laminar flow based unsteady state friction term [21] can be successfully used for low Reynolds number turbulent water hammer (Re0<104). The frozen viscosity convolutional model captures the multidimensional shape of the flow velocity profile, which is represented as an average velocity in the 1D model. The multi-dimensional CFD represents a very effective and accurate tool for the analysis of flow and pressure conditions in real engineering problems which was confirmed by 3D CFD simulation of the water hammer. The 3D CFD numerical model can be further extended with two-way interaction between liquid and solid (direct coupling). High pressure loads cause high pipe-wall stress and thus elastic deformation of the pipeline, which causes a reverse change in the flow. This would give us a tool that provides a precise insight into 3D pipe flow conditions. The analytical 1D model is an excellent tool for the verification studies of numerical 1D and 3D models. The results of calculations based on the two 1D and 3D numerical models agree well with the results of the measurements obtained in a simple pipeline apparatus. Quantitative comparison analysis yields up to 2 % difference in maximum head at the valve and up to 5 % relative difference in pressure head drop at the midpoint of the pipeline monitored over the first four positive pressure pulses. Naturally, 1D models are computationally less intensive compared to a 3D CFD model. Further research is sought in theoretical (analytical solution for turbulent water hammer) and experimental investigations of high Reynolds number water hammer (Re0>105), and in the future, in direct numerical simulations of water hammer events [27]. References [1] Wylie, E.B., Streeter, V.L. Fluid Transients in Systems. Prentice Hall (1993) Englewood Cllifs. [2] Chaudhry, M.H. Applied Hydraulic Transients. Springer (2014) New York, DOI:10.1007/978-1-4614-8538-4. [3] Giljen, Z., Nedeljkovic, M., Cheng, Y. The influence of pump-turbine specific speed on hydraulic transient processes. Stroj vestn-J Mech E 70 231-246 (2024) DOI:10.5545/sv-jme.2023.776. [4] Jung, B.S., Karney, B.W. Pressure surge control strategies revised. AWWA Wat Sci 2 e1169 (2020) DOI:10.1002/aws.1169. [5] Li, H., Xu, B., Arzaghi, E., Abbassi, R., Chen, D., Aggidis, G.A., Zhang, J., Patelli, E. Transient safety assessment and risk mitigation of a hydroelectric generation system. Energy 196 117135 (2020) DOI:10.1016/j.energy.2020.117135. [6] Kubrak, M. Modeling of closing functions for gate valves fitted with V-ports. J Pipel Syst Eng Pract 15 04024023 (2024) DOI:10.061/JPSEA2.PSENG-1588. [7] Brunone, B., Golia, U.M., Greco, M. Effects of two-dimensionality on pipe transients modelling. J Hydraul Eng 121 906-912 (1995) DOI:10.1061/ (ASCE)0733-429(1995)121:12(906). [8] Bergant, A., Simpson, A.R., Vítkovský, J. Developments in unsteady pipe flow friction modelling. J Hydraul Res 39 249-257 (2001) DOI:10.1080/00221680109499828. [9] Pezzinga, G. Evaluation of unsteady flow resistances by quasi-2D or 1D models. J Hydraul Eng 126 778-785 (2000) DOI:10.1061/(ASCE)0733­9429(2000)126:10(778). [10] Riedelmeier, S., Becker, S., Schlücker, E. Damping of water hammer oscillations -comparisons of 3D CFD and 1D calculations using two selected models for pipe friction. PAMM Proc Appl Math Mech 14 705-706 (2014) DOI:10.1002/ pamm.201410335. [11] Mandair, S., Magnan, R., Morissette, J.F., Karney, B. Energy-based evaluation of 1D unsteady friction models for classic laminar water hammer with comparison to CFD. J Hydraul Eng 146 04019072 (2020) DOI:10.1061/(ASCE)HY.1943­7900.0001697. [12] Martins, S.C., Covas, D.I.C., Capponi, C., Meniconi, S., Brunone, B. Unified approach for damping rate to transient laminar flow: experiments, computational fluid dynamics, and one dimensional, and global models. J Fluid Eng-T ASME 146 021306 (2024) DOI:10.1115/1.4063697. [13] Plouraboué, F. Review on water-hammer waves mechanical and theoretical foundations. Eur J Mech B-Fluids 108 237-271 (2024) DOI:10.1016/j. euromechflu. 2024.08.001. [14] Kumar, M.R.A., Pu, J.H., Hanmaiahgari, P.R., Lambert, M.F. Insights into CFD modelling of water hammer. Water 15 3988 (2023) DOI:10.3990/w15223988. [15] Urbanowicz, K., Jing, H., Bergant, A., Stosiak, M., Lubecki, M. Progress in analytical modeling of water hammer. J Fluids Eng-T ASME 145 081203 (2023) DOI:10.1115/1.4062290. [16] Neyestanaki, M.K., Dunca, G., Jonsson, P., Cervantes, M.J. A comparison of different methods for water hammer valve closure with CFD. Water 15 1510 (2023) DOI:10.3990/w15081510. [17] Zhang, X. Transient flow characteristics of a pressure differential valve with different valve spool damping orifice structures. Stroj vestn-J Mech E 70 141-158 (2024) DOI:10.5545/sv-jme.2023.691. [18] Saemi, S., Raisee, M., Cervantes, M.J., Nourbakhsh, A. Computation of two- and three-dimensional water hammer flows. J Hydraul Res 57 386-404 (2019) DOI:10 .1080/00221686.2018.1459892. [19] Cao, Y., Zhou, L., Qu, C., Fang, H., Liu, D. 3D CFD simulation and analysis of transient flow in a water pipeline. AQUA - Water Infrastr Ecosyst Soc 71 751-767 (2022) DOI:10.2166/aqua.2022.023. [20] Wang, C., Nilsson, H., Yang, J., Petit, O. 1D-3D coupling for hydraulic system transient simulations. Comput Phys Commun 210 (2017) DOI:10.1016/j. cpc.2016.09.007. [21] Zielke, W. Frequency-dependent friction in transient pipe flow. J Basic Eng-T ASME 90 109-115 (1968) DOI:10.1115/1.3605049. [22] ANSYS® Academic Research CFD, Release 18.2. Ansys (2019) Canonsburg. [23] Vardy, A.E., Brown, J.M.B., He, S., Ariyarante, C., Gorji, S. Applicability of frozen-viscosity models of unsteady wall shear stress. J Hydraul Eng 141 0401064 (2015) DOI:10.1061/(ASCE)HY.1943-7900.0000930. [24] Urbanowicz, K. Fast and accurate modelling of frictional transient pipe flow. J Appl Math Mech 98 802-823 (2018) DOI:10.1002/zamm.201600246. [25] Vardy, A.E., Brown, J.M.B. Transient turbulent friction in smooth pipe flows. J Sound Vib 259 1011-1036 (2003) DOI:10.1016/jsvi.2002.5160. [26] Vardy, A.E., Brown, J.M.B. Transient turbulent friction in fully-rough pipes flow. J Sound Vib 270 233-257 (2004) DOI:10.1016/S0022-460X(03)00492-9. [27] Guerrero, B., Lambert, M.F., Chin, R.C. Extension of the 1D unsteady model for rapidly accelerating and decelerating turbulent pipe flows. J Hydraul Eng 148 04022014 (2022) DOI:10.1061/(ASCE)HY.1943-7900.0001998. [28] Meniconi, S., Duan, H.F., Brunone, B., Ghidaoui, M.S., Lee, P.J., Ferrante, M. Further developments in rapidly decelerating turbulent flow modelling. J Hydraul Eng 140 04014028 (2014) DOI:10.1061/(ASCE)HY.1943-7900.0000880. [29] Rich, G.R. Water hammer analysis by the Laplace-Mellin transformation. J Fluid Eng-T ASME 67 361-376 (1945) DOI:10.1115/1.4018265. [30] Urbanowicz, K., Bergant, A., Stosiak, M., Karpenko, M., Bogdevicius, M. Developments in analytical wall shear stress modelling for water hammer phenomena. J Sound Vib 562 117848 (2023) DOI:10.1016/j.jsv.2023.117848. [31] Trikha, A.K. An efficient method for simulating frequency-dependent friction in transient liquid flow. J Fluid Eng-T ASME 97 97-105 (1975) DOI:10.1115/1.3447224. [32] Hullender, D.A. Water hammer peak pressures and decay rates of transients in smooth lines with turbulent flow. J Fluid Eng-T ASME 140 061204 (2018) DOI:10.1115/1.4039120. [33] García García, F.J., Alvariño, P.F. On an analytic solution for general unsteady/ transient turbulent pipe flow and starting turbulent flow. Eur J Mech B-Fluids 74 200-210 (2019) DOI:10.1016/j.euromechflu.2018.11.014. [34] Tucker, P.G. Advanced Computational Fluid and Aerodynamics. Cambridge University Press (2016) New York, DOI:10.1017/CBO9781139872010. [35] Wilcox, D.C. Turbulence Modeling for CFD. DCW Industries (2010) La Cañada. [36] Martins, N.M.C., Wahba, E.M. On the hierarchy of models for pipe transients: From quasi-two-dimensional water hammer models to full three-dimensional computational fluid dynamics models. J Press Vess-T ASME 144 021402 (2022) DOI:10.1115/1.4051930. [37] Martins, N.M.C., Carriço, N.J.G., Ramos, H.M., Covas, D.I.C. Velocity-distribution in pressurized pipe flow using CFD: Accuracy and mesh analysis. Comput Fluids 105 218-230 (2014) DOI:10.1016/j.compfluid.2014.09.031. [38] Tait, P.G. Report on Some of the Physical Properties of Fresh Water and of Sea Water. Challenger Expedition, Johnson Reprint Corporation (1965) New York. [39] Bergant, A., Simpson, A.R. Water Hammer and Column Separation Measurements in an Experimental Apparatus. Research Report no. R128, Department of Civil and Environmental Engineering, University of Adelaide (1995) Adelaide. [40] Coleman, H.W., Steele, W.G. Experimentation and Uncertainty Analysis for Engineers. John Wiley and Sons (1989) New York. [41] Bergant, A., Karadžic, U., Vítkovský, J., Vušanovic, I., Simpson, A.R. A discrete gas-cavity model that considers the frictional effects of unsteady pipe flow. Stroj vestn-J Mech E 51 692-710 (2005). [42] Duan, H.-F., Ghidaoui, M.S., Lee, P.J., Tung, Y.K. Relevance of unsteady friction to pipe size and length in pipe fluid transients. J Hydraul Eng 138 154-166 (2012) DOI:10.1061/(ASCE)HY.1943-7900.0000497. Acknowledgements Anton Bergant gratefully acknowledges the support of the Slovenian Research and Innovation Agency (ARIS) conducted through the research programme P2-0126. Kamil Urbanowicz declares that his contribution is financially supported by the Minister of Science under the ‘Regional Initiative of Excellence’ (RID) programme. In addition, this paper has been written to mark the 30th anniversary of the Slovenian Association of Hydraulic Research (SDHR). Received: 2024-10-01, revised: 2025-05-12, accepted: 2025-05-27 as Original Scientific Paper Data availability Data supporting the findings of this study are available from the corresponding author upon reasonable request. Authors Contributions Anton Bergant: Conceptualization, Metodology, Computational Analysis, Measurements, Validation, Writing – Original Draft, Writing – Review and Editing, Supervision; Zlatko Rek: Conceptualization, Metodology, Computational Analysis, Validation, Writing – Original Draft, Writing – Review and Editing; Kamil Urbanowicz: Conceptualization, Metodology, Computational Analysis, Validation, Writing – Original Draft, Writing – Review and Editing. Analiticne, numericne 1D in 3D raziskave vodnega udara v enostavnem preizkusnem cevovodu Povzetek Prispevek obravnava analiticno in numericno simulacijo vodnega udara v sistemu rezervoar-cevovod-ventil (RCV). Razvita je bila analiticna rešitev enacb vodnega udara z nestalnim clenom stenskega trenja za prehodni laminarni tok v RCV sistemu. Za simulacijo poljubnega pretocnega stanja je bilo razvitih vec 1D numericnih metod. Fizikalno zasnovana metoda karakteristik se je izkazala za racunsko ucinkovito in prikladno za obravnavanje kompleksnih robnih pogojev. Natancnost 1D numericnega modela se poveca z vpeljavo clenov, ki zajemajo vpliv vecrazsežnega prostora (primer nestalnega stenskega trenja). 3D model neposredno napoveduje te vplive in predstavlja odlicno orodje za raziskovanje vecrazsežnih lastnosti tekocin (numericni laboratorij). Rezultati izracunov, ki temeljijo na 1D in 3D numericnih modelih, se dobro ujemajo z rezultati meritev, ob upoštevanju ustrezne napovedi in modeliranja vplivnih fizikalnih kolicin v stanjih laminarnega in turbulentnega vodnega udara z nizkim Reynoldsovim številom. Kvantitativna primerjalna analiza pokaže do 2 % razliko v maksimalni tlacni višini pri ventilu in do 5 % relativno razliko v padcu tlacne višine na sredini cevovoda, ki se spremlja v prvih štirih tlacnih pulzih. Kljucne besede cevovod, vodni udar, analiticna rešitev, metoda karakteristik, racunalniška dinamika tekocin, nestalno stensko trenje Study on the Influence of the Splitter Blade Length on Radial and Axial Force of a Centrifugal Pump Qijiang Ma1,2,3,4 – Zhenbo Liu4 – Sen Jiang5 1 School of Mechanical and Electrical Engineering, Chuzhou University, China 2 School of Mechanical Engineering, Hefei University of Technology, China 3 Anhui Liuxiang Special Ship Co., Ltd., China 4 National Research Center of Pumps, Jiangsu University, China 5 Ningbo Yi Ji Valve Manufacturing Co., Ltd., China qijma@chzu.edu.cn Abstract To enhance the operational stability of centrifugal pumps, this study investigates the influence of splitter blade length on the axial and radial forces of centrifugal pumps. Using the SST k-. turbulence model and experimental research, the external characteristics, axial force, radial force, and time-frequency characteristics of pressure pulsation were compared among impellers without splitter blades and those with splitter blades of two different lengths. The results show that impellers with a conventional structure achieve higher efficiency near the design operating point. However, under low-flow conditions, the rectifying effect of splitter blades allows impellers with splitter blades to achieve higher efficiency. For impellers with splitter blades, the axial force shows a periodic behavior, presenting two peak values — one large and one small — within each cycle. The addition of splitter blades shifts the radial force acting on the impeller to one side, requiring an increase in the support strength and stiffness of the rotor system. Furthermore, splitter blades influence the pressure pulsation at the impeller's inlet and outlet, the original impeller has an additional pulsation energy at 145 Hz but the splitter blades are different. Thus, an impeller with splitter blades can reduce the frequency pulsation and optimizing flow conditions. This study provides valuable theoretical insights and data support for the hydraulic structural optimization of centrifugal pumps. Keywords centrifugal pump, axial force, radial force, splitter blade length, pressure fluctuation Highlights . Influence of splitter blade length on the axial and radial forces of centrifugal pump is studied. . Impellers with a conventional structure achieve higher efficiency near the design operating point. . Splitter blades shifts the radial force acting on the impeller to one side. . Splitter blades influence the pressure pulsation at the impeller's inlet and outlet. 1 INTRODUCTION Centrifugal pumps, as one of the most commonly used fluid transportation devices, are widely applied in various fields such as water treatment [1], chemical production [2], energy conversion [3], and agricultural irrigation [4]. Their operational efficiency and stability directly impact the economic viability and safety of production processes. The fundamental working principle of centrifugal pumps is based on the action of centrifugal force. By rotating, the impeller draws fluid from the pump inlet into the blades and accelerates it to a high specific velocity under centrifugal force, simultaneously converting kinetic energy into pressure energy [5]. As such, the impeller serves as the core component of a centrifugal pump. In the design and operation of centrifugal pumps, dynamic characteristics are a critical factor that cannot be overlooked. Al-Obaidi used numerical calculation and experiment to study the complex hydrodynamic characteristics of axial flow pumps [6], including the influence of flow conditions [7], the influence of the number of blades [8], guide vane configurations [9], and the influence of cavitation characteristics [10]. Moreover, during operation, the blades are subjected to fluid dynamics forces, resulting in radial and axial forces [11]. Radial force primarily acts perpendicular to the rotational axis on the impeller and is caused by the centrifugal force of the fluid and the interaction between the impeller and the fluid. Axial force, on the other hand, is generated along the pump shaft and is mainly influenced by the flow direction of the fluid within the pump. Notably, the inlet and outlet angles of the pump play a crucial role in the distribution of these forces [12]. Radial and axial forces are the primary factors affecting the stability of centrifugal pumps. Excessive radial force can increase friction between the pump casing and the impeller, leading to heat generation and localized wear [13]. Similarly, excessive axial force may cause bending of the pump shaft, compromising operational stability and reducing the pump's lifespan. Therefore, effectively controlling these forces is a key task in centrifugal pump design. In recent years, advancements in computational fluid dynamics (CFD) have allowed researchers to delve deeper into the factors influencing radial and axial forces [14]. For instance, Zhu et al. [15] investigated the axial force of a guide vane mixed-flow pump using experimental and numerical simulation methods. They optimized the axial force by incorporating balance holes and balance disc structures, employing orthogonal experiments and BP neural network algorithms to achieve optimal efficiency and axial force. Liu et al. [16] studied the axial force characteristics of a centrifugal pump with a floating impeller structure, developing mathematical expressions for fluid leakage, fluid pressure, and axial force acting on a stainless steel disc in the axial gap. Their findings revealed that the floating impeller significantly reduces axial force. Dong et al. [17] examined the impact of particle flow on the axial force of a centrifugal pump, using a semi-open impeller pump as the study object. Their research showed that the total axial force in clear water conditions is greater than in particle flow, but the total axial force increases as the solid-phase volume fraction rises. Cao et al. [18] investigated the influence of the backflow hole structure on the pump's axial force and energy characteristics, finding that when the ratio of backflow balance hole area to sealing ring area exceeds six, better axial force balance can be achieved. Similarly, significant progress has been made in the study of impeller radial forces. Cui et al. [19] explored the effect of radial force on the vibration characteristics of centrifugal pumps using numerical simulations to examine the relationship between flow rate, radial force, and vibration characteristics. Their results showed that impeller radial displacement aligns with radial force variations but lags behind in phase. Moreover, radial force-induced vibrations increase flow losses within the pump, reducing efficiency. González et al. [20] studied the impact of the clearance between the impeller and volute on radial forces in centrifugal pumps through experimental and numerical analysis of both steady and unsteady radial forces. They found that, under certain flow conditions, larger diameters generate higher radial forces, although this does not apply universally to all operating points. Ou et al. [21] conducted numerical studies on the hydrodynamic radial force in mixed-flow pumps, discovering that the transient hydrodynamic radial force of the impeller shows periodic variation over time. Under uniform inflow conditions, the average transient radial force approaches zero, while in certain flow conditions, recirculation flow structures significantly affect pressure pulsation and hydrodynamic radial force. The aforementioned studies highlight that the hydraulic structure of the impeller has a substantial impact on the axial and radial forces within the entire rotor system, further influencing the vibration and noise characteristics of the pump's operational system. Consequently, reducing axial and radial force pulsations has remained a key focus for both academia and industry. In recent years, with advancements in the hydraulic design methods of centrifugal pumps, the splitter blade design approach has been introduced to enhance pump performance [22]. Splitter blades are secondary blades added behind the main flow path of the pump. By altering the flow path on the blades, they optimize the fluid flow state, reduce flow losses, and thus improve the performance and stability of the pump [23]. In centrifugal pumps, the design parameters for splitter blades are particularly critical. Different configurations of blade shape, length, and angle can lead to variations in fluid flow characteristics within the pump, thereby affecting the operational state and lifespan of the entire system [24]. For instance, Yan et al. [25] employed unsteady CFD analysis to enhance the efficiency of a double-volute centrifugal pump and reduce the unsteady radial forces on the impeller. By introducing splitter blades, they improved the pump's hydraulic performance and mitigated unstable radial forces. Similarly, Zeng et al. [26], aiming to reduce pressure and radial force pulsations at the design point of a centrifugal pump, considered the circumferential position, leading-edge location, and deflection angle of splitter blades. They demonstrated experimentally that splitter blades are effective in suppressing secondary flows and reducing fluid-induced vibrations. Xie et al. [27] studied the resonant response of impellers to multi-frequency fluid excitations using an acoustic­fluid-structure coupling method. They investigated the influence of the number and circumferential position of splitter blades on the natural frequency of impellers. Their findings revealed that adding splitter blades between the main impeller blades reduced the natural frequency of the impeller in both air and water, while the circumferential position of the splitter blades had minimal impact on the natural frequency. These studies collectively indicate that splitter blades significantly enhance the internal flow structure of pumps. However, the presence of splitter blades also influences the distribution characteristics of radial and axial forces on the impeller, particularly with respect to blade length [28]. Existing research suggests that different blade lengths can alter the magnitude and distribution of radial and axial forces [29]. Appropriately designed blade lengths can not only improve pump efficiency and flow rate but also effectively reduce vortex flows within the pump, minimize noise, and enhance stability [30]. Conversely, improperly designed blade lengths may increase these forces, leading to greater vibrations, increased noise, and premature equipment wear [31]. While these findings primarily focus on the main blades of the impeller, they also provide valuable insights for the design of splitter blades. Currently, research on splitter blades remains limited, and studies specifically addressing the effects of splitter blade length on impeller axial and radial forces have yet to be published. Therefore, an in-depth investigation of the hydrodynamic characteristics of splitter blades and their impact on centrifugal pump performance not only offers guidance for pump optimization and improved energy efficiency but also provides critical theoretical and practical references for the development of related industries. This study will employ a combination of experimental and numerical simulation methods. By constructing centrifugal pump models with varying blade lengths, the research will investigate the variations in radial and axial forces under different operating conditions. Numerical simulations will be conducted using CFD software to capture the actual fluid flow characteristics within the pump, and these results will be compared with experimental data. Simultaneously, the experimental data will be used to validate the accuracy of the simulation model and further analyze the impact of different design parameters on pump performance. This systematic research approach not only allows for a comprehensive analysis of the effects of blade length on the pump’s dynamic characteristics but also provides practical guidance for future studies and applications. 2 METHODS & MATERIALS 2.1 Pump Model and Splitter Blade Structure In this paper, we focus on the ISG25-220 model, a single-stage single-suction horizontal centrifugal pump, as our primary subject of research, illustrated in Fig. 1. The centrifugal pump is designed to have a flow rate (Q) of 88 m3/h , a head (H) of 12 m, and an operational speed (n) of 1450 rpm. For the design of the pump's hydraulic structure, we utilize a high-efficiency centrifugal pump model with a comparable specific speed, which is calculated as follows: 365 . n= = 21. (1) s 34/ H We primarily design two schemes: one with shunt blades and the other without, for comparative purposes. Additionally, to evaluate the impact of diverter blade length on the hydraulic performance of the centrifugal pump, we implement two different lengths of diverter blades for further analysis. Consequently, three distinct hydraulic design schemes for the impeller are selected, as depicted in Fig. 2. Fig. 1. Physical centrifugal pump model Fig. 2. Splitter blades structure; a) impeller #A, b) impeller #B, c) impeller #C The rotary bias method is mainly used in the design of the shunt blade. Since the original non-shunt impeller has 6 blades, the angle between each blade is 60°. Therefore, when designing the splitter blade, the original blade is rotated by 30° and its inlet edge is cut, leaving only the tail of the blade. In this paper, the length of the splitter blade is set to 0.7 times and 0.3 times of the original blade length respectively. The main parameters of the centrifugal pump model and the splitter blade are shown in Table 1. Table 1. Main parameters of the centrifugal pump model Parameter Value Designed flow rate, Qd [m3/h] 88 Designed rated head, H [m] 12 Rated speed, n [rev/min] 1450 Impeller inlet diameter, D1 [mm] 112.5 Impeller outer diameter, D2 [mm] 219 Number of blades, Z 6 Impeller outlet width, b2 [mm] 22.02 Blade inlet angle, ß1 [°] Blade outlet angle, ß2 [°] 16.5 Blade wrap angle, . [°] 90 Guide vane angle, a [°] 15 Volute base circle diameter, D3 [mm] 176.5 Volute outlet diameter, D4 [mm] 100 2.2 Pre-Processing of Numerical Calculation Some preparatory work is required before numerical calculations are performed on a centrifugal pump model. The first is to model the computational domain through 3D design software. The hydraulic power of the model pump was modeled by UG NX software in three dimensions, including inlet section, outlet section, impeller, volute and pump gap chamber. In the three schemes in this paper, the clearance of the impeller mouth ring is equal and has little influence on the flow field in the shunt impeller, so it can be ignored in the modeling process. In addition, considering that inlet reflux and outlet reflux will lead to divergence of the calculation process, the inlet and outlet areas of the impeller are extended respectively [32]. The three-dimensional water structure of the centrifugal pump is shown in Fig. 3. Grid is the carrier of simulation and analysis, and the quality of grid has an important impact on the calculation accuracy and efficiency [33]. In this paper, ANSYS ICEM is used to divide the hexahedral structure of the whole basin. In the process of discretization of the calculation domain, the mixed mesh method is mainly adopted, that is, hexahedral structured mesh is mainly used for the inlet segment, outlet segment, volute and pump cavity, which can not only accurately control the streamline distribution and the orthogonal direction of the boundary layer, but also flexibly adjust the distance between nodes to adjust the density of the boundary layer mesh [34]. As for the impeller, due to its high hydraulic structure complexity, the generated hexahedral mesh quality is poor, so the tetrahedral unstructured mesh is used to generate. For the three impeller schemes that need to be compared in this paper, the same mesh generation standard is adopted, and the final mesh quality is higher than 0.5. Figure 3 shows the meshing scheme and encryption details for different computing domains. Fig. 4. Mesh details of different simulation domain The presence of boundary layer mesh can more accurately capture the flow in the near-wall region and reflect the flow characteristics in the near-wall region [35]. Thus, it is necessary to encrypt the boundary layer of blade surface and volute tongue to ensure calculation accuracy. In order to meet the requirements of different flow patterns on the number of grids in the near wall area, in order to ensure that there are enough nodes in the area, y+ value can be used to test the location of the node closest to the wall. Here, y+ value represents the distance from the node closest to the wall to the wall, and is a dimensionless variable, whose definition is as follows [36]: ˜ uy y ° , (2) v where, u represents the wall friction velocity, [m/s]; y represents the distance between the two nodes closest to the wall, [m]; v represents kinematic viscosity, [Pa·s]. The y+ values of blade surface and volute wall are shown in Fig. 5. Fig. 5. Distribution of y+ on the blade 2.2 Grid Independence Test In the numerical calculation, it is necessary to analyze the influence of the number of grids on the calculation results to avoid the solution error caused by the size of the grid, that is, the number of grids [37]. Thus, it is necessary to use different mesh sizes for grid division. The grid independence analysis is carried out based on the change of the head with the number of grids under the design condition, and an appropriate number of grids is selected for further study. To evaluate the sensitivity of the computational mesh, this study adopts the grid convergence index (GCI), a method initially introduced by Li et al. [38]. The GCI approach necessitates that the obtained results adhere to the condition of monotonic convergence to ensure reliability. In this work, three distinct grids with varying levels of refinement were constructed and analyzed under the rated operating conditions of the original impeller design. This process enables a systematic assessment of how mesh resolution influences the accuracy and stability of the simulation results. The formula for calculating GCI is as follow [39]: . b.1,b GCI ˜ | FSp , (3) r °1 b.1,b where Fs is the safety factor, taking values in the range 1.25 to 3.00, and 1.25 when three or more grids are used. Symbol r denotes the grid refinement ratio and is defined as follows: 1 . N ˆD b˜1rb˜1,b °. . , (4) N . b . where N is the number of control cells of the grid; D is the computational dimension; the subscript b represents different grid schemes, with larger b indicating a denser grid. Symbol f is the relative error of the results obtained from numerical calculations using two sets of grids with adjacent quantities: f . f bb°1˜b°1,b . , (5) fb°1 where fb is the numerical discrete solution of the selected convergence parameter. Symbol . is the convergence accuracy. When three sets of grids are used for GCI analysis, . can be calculated iteratively based on Eq. (6) as follows: .f3 .f2 .˜° 1 ln ˆ g ˜ , (6a)   ln() r f .f 21, . 21 . ˆr21 ˜ , .s  g°.˜.ln..˜ , (6b) r .s .32,  . f ° f ˆ s ˜ sgn . 32 ., (6c) f ° f . 21 . where g(.) is the .-order error term coefficient that does not vary with the grid. It should be noted that when ( f3 – f2) / ( f2 – f1) < 0, i.e., s = -1, the oscillatory iterative convergence solution is prone to occur, and then the grid needs to be replaced and recreated. As can be seen from Table 2, with the increase of the number of grids, the value of GCI decreases from 1.21 % to 0.62 %, meeting the requirement that the value required by the GCI convergence standard is less than 1. Therefore, in the final calculation scheme, the grid density scheme is selected for grid generation in all calculation domains. Table 2. Calculated results of GCI with different sizes of grids Mesh number×106 Relative error GCI [%] 20.3 0.000097 1.21 23.6 30.5 0.000034 0.62 2.2 Governing Equations and Boundary Conditions The numerical calculation process of the flow field in centrifugal pump is actually a process of solving the governing equation, and in order to solve the governing equation, it is necessary to select a suitable turbulence model closed equations. Reynolds time average method uses time mean value and pulsation value to represent instantaneous value, which has high reliability and small calculation amount, and can basically meet the engineering application requirements, so it is the most widely used in practical engineering [40]. The turbulence model based on Reynolds-averaged Navier– Stokes (RANS) equation can be divided into vortex viscosity model and Reynolds stress model. The vortex viscosity model is mainly used to solve the turbulence viscosity coefficient, which can be divided into zero equation model, single equation model and two-sided equation model, among which the two-sided equation model is the most widely used. The standard Reynolds stress model is mainly used to solve the Reynolds stress transport equation [41] Because the SST k-. model uses a mixed function combining the standard k-. and the transformed k-e model, the standard k-. model is used in the near wall region and the transformed k-e model is used away from the near wall region. At the same time, considering the transmission of turbulent shear force in the internal flow field, the flow separation amount of fluid under negative pressure gradient can be predicted very accurately, without excessive prediction of eddy viscosity, which is especially suitable for the simulation of boundary layer requiring high precision. Therefore, the SST k-. model closed control equations are used in this paper to numerically calculate the internal flow field of low specific speed centrifugal pump under different working conditions. In the SST k-. model, the eddy viscosity coefficient and k 2.3 Experimental Equipment equation and k equation are expressed as ˜(k) ˜(ku ) ˜. .˜k  it * ° . ˆ° ° Pk  k , (7) ˜t ˜x ˜x ˜x ii . k3 . j  ˜() ˜(u) ˜ . .˜  ° i .° t ° P ˆ  3 k ˜t ˜x ˜x ˜xk ij. 3 . j 21 ˜k˜  3 ° 21(  F1) , (8)  ˜xx˜  2 jj ak1 ˜ .° , (9) t max(a1., SF2) .°ui °uj .°ui 2 °uk .°uk . Pk ˜t .... .3t .k , (10) °x °x °x 3 °x °x ˆ jj  jk ˆ k  where, F1 and F2 serve as blending functions within the SST k–. turbulence model, while S denotes the shear force constant. The model utilizes an empirical constant ß, specifically set to 0.09, alongside a1, another empirical constant. The turbulent eddy frequency, . [s.¹], while the turbulent kinetic energy, k [m²/s²]. These parameters work together to accurately represent and predict flow behavior within the model's two-equation framework. The setting of boundary conditions will essentially affect the accuracy of the calculation results. The steady conditions of the low specific speed centrifugal pump model are set as follows: the inlet boundary condition is set as the total pressure, the outlet boundary condition is set as the mass flow, the impeller is set as the frozen rotor when it contacts with the inlet section and the impeller and the volute, and the rotation Angle is specified as 360°. Set the number of steady calculation steps of the model to 2500. The unsteady setting of the outlet and inlet of the model is the same as that of the steady setting. The difference is that the impeller wall is set to instantaneous rotation, the inlet of the impeller and the contact surface between the volute and the impeller are set to instantaneous freezing rotor total time 0.103448 s, and the time step is t =1.72414×10 s, that is, the impeller of the centrifugal pump model rotates 5 times. The impeller is turned 3° at each time step. To ensure simulation accuracy, a convergencecriterionof10.4wasselectedforboththecontinuityand momentum equations. Based on the open test bench of Zhenjiang Machinery, China industry testing institution, the performance test of the pump was carried out. Figure 4 shows a schematic diagram of the test bench. First let the motor no load, torque calibration at rated speed, and then connect the motor to the test pump shaft. Connect each test device and instrument according to the test standard, then fill the test pipe with normal temperature water, turn on the machine to check whether the instrument and equipment can work normally, open the pipe vent, and drain the air in the pipe. Before starting the test mixed-flow pump, the inlet and outlet valves on the test pipe are adjusted to the maximum opening, and then the speed is adjusted by the frequency conversion regulator to make the speed reach and maintain at the rated speed of 1450 r/min. Keep the inlet valve open to the maximum, slowly reduce the outlet valve open, parameter measurement under different flow conditions. Observe the indicator number on the pump product parameter measuring instrument and system analysis program. When the indicator number is stable at a certain value or fluctuates back and forth within a small range of a certain value, it is considered as the measured value and the corresponding data is recorded. After the data measurement of the minimum flow condition is completed, the speed is still adjusted by the variable frequency regulator to slowly reduce to zero. After the shutdown is complete, wait for the water flow in the pipeline to stop, and then turn the outlet valve open to the maximum again. In order to minimize the accidental error and improve the reliability of the test data, the energy characteristic data recorded in the three repeated tests were arithmetically averaged and used as the final test result. 2.4 Verification of Simulation The pump data acquisition terminal and processing program can calculate the head of the test pump according to the inlet and outlet pressure of the test pipeline, the shaft power can be calculated according to the torque and speed data measured by the torque and speed measuring instrument, and the efficiency can be calculated according to the head and shaft power. The three expressions are as follows: . p˜ pˆ. v2 ˜ v2 ˆ 21 21 HZ=2 ˜ Z1 °. .°.., (11) g 2g . ... nM M. 0 . °. P˜ , (12) 9552 Fig. 6. Schematic diagram of experiment system; 1 water tank, 2 valve, 3 inlet pressure transmitter, 4 pump, 5 torque sensor, 6 motor, 7 outlet pressure transmitter, 8 torque speed transmitter, 9 flowmeter, 10 valve, 11 data acquisition terminal, and 12 computer °gQH ˜. .100 %. (13) 1000P In Equation (11), p1 and p2 are respectively pump inlet and outlet pressure, [Pa]; Z1 and Z2 are the installation heights of pump inlet and outlet respectively, [m]. In this test, the installation heights of pump inlet and outlet are the same, so Z1= Z2; v1 and v2 are respectively the average flow rate of the inlet and outlet of the pump, [m/s], which can be calculated from the inner diameter of the inlet and outlet pipeline; M is the measured torque of the pump set, M0 is the no-load torque, [N·m]; n is the speed, [r/min]; Q is the flow rate, [m3/h]; H is the head, [m]; P is the shaft power of the mixed flow pump, [kW]. Fig. 7. Comparison of simulations and experimental results Figure 6 shows the performance curves obtained from numerical calculations and tests of the model pump. From the trend of the head and efficiency, the numerical and experimental results are in perfect agreement. At about Q = 95 m3/h, the efficiency of the centrifugal pump has the maximum value and is the optimal working condition point. With the decrease in flow condition, the head curve gradually increases. At about Q =45 m3/h, the local lowest point of the lift curve appears, which is the point of stall condition, but the drop of the lift is not much, indicating that the flow loss is not serious. From the overall change of the head and efficiency, the head curve changes gently in the low flow condition, indicating that the flow structure in the pump shows a certain trend change. By comparing the experimental and numerical results, it can be found that the error between the numerical results and the experimental values is less than 10 % under each flow condition, and further the error is under 0.2 for each test, which meets the accuracy required by the numerical calculation. Therefore, the numerical method adopted in this paper is feasible. 3 RESULTS AND DISCUSSION 3.1 Comparison of Energy Performance Figure 8 compares the head and efficiency curves of centrifugal pumps with and without splitter blades. The results reveal that pumps with splitter blades exhibit distinct performance characteristics compared to conventional impellers, especially under low-flow operating conditions. Unlike the consistent negative slope observed in conventional impellers, the head curve of pumps with splitter blades transitions from a negative slope to a positive slope, indicating the presence of a noticeable extremum point. In terms of specific head values, centrifugal pumps equipped with splitter blades achieve significantly higher head performance than conventional impellers. It is almost 8 % higher than the original impeller design under the condition of 70 m3/h. For longer splitter blade structures, the head is slightly higher, although the difference is marginal. Regarding efficiency curves, conventional impellers demonstrate higher the efficiency near the design operating point, whereas impellers with splitter blades achieve better efficiency under low-flow conditions due to the rectifying effect of the splitter blades. Consequently, while splitter blades positively influence the head curve of centrifugal pumps, this enhancement comes at the cost of some efficiency loss. 3.2 Analysis of Flow Field of Centrifugal Pump Blades with Different Splitter Blade Length Figures 9 and 10 illustrate the pressure and velocity distributions within the impeller for different impeller designs. As shown in the figures, both pressure and velocity are relatively low at the impeller inlet but gradually increase as the fluid moves toward the impeller outlet, reaching their peak values at the exit. Comparing different impeller designs, it is evident that the addition of splitter blades results in a more uniform pressure distribution along the impeller Fig. 8. Impact on the pump performance of splitter blades; a) head, and b) efficiency Fig. 9. Pressure distribution in different kinds of impeller periphery. Notably, near the volute tongue, the high-pressure region is significantly reduced, and the circumferential pressure uniformity is improved. This explains why the splitter blades contribute to a higher head in the centrifugal pump. However, the introduction of splitter blades also leads to the formation of a distinct low-velocity region near the tongue, which plays a crucial role in the development of vortex structures. The presence of these vortices increases hydraulic losses within the pump, ultimately leading to a reduction in efficiency under the design operating conditions. 3.3 Analysis of Axial Force of Centrifugal Pump Blades with Different Splitter Blade Length Figure 10 illustrates the schematic diagram of the axial force acting on the impeller of the centrifugal pump under design operating condition. According to the pump fluid dynamics theory, the axial force of the impeller primarily originates from the fluid within the impeller flow passage and the fluid in the gaps between the front and rear shrouds of the impeller. Zhu et al. [42], developed a mathematical model of the axial force exerted on the impeller by the fluid in the front and rear gap regions under the rotational action of the impeller. The main expression of the model is as follows: 0 * ** 22 F˜°F.F˜°.r 2ˆrP..1 .r. dr 12 0 r 2 1 1 22 . 2ˆrP...r0 . dr. (14) .r0 r 2 2 Here, F1 represents the axial force acting on the front shroud of the impeller, while F2 denotes the axial force acting on the rear shroud. However, this model only accounts for the axial force caused by the gap fluid and needs to include the contribution of the fluid within the impeller flow passage to represent the total axial force acting on the impeller. Although some researchers have proposed empirical formulas to calculate the axial force caused by the fluid within the impeller flow passage, in this study, more accurate axial force values can be directly obtained through post-processing in simulation software. Therefore, the final axial force is the sum of three components: the axial force on the front shroud, the axial force on the rear shroud, and the axial force caused by the fluid within the impeller flow passage. Figure 12 presents th e time-dependent variations in the axial force acting on the impeller under design conditions for different impeller structures. To minimize random errors, the axial force over the last three revolutions at the same impeller phase was averaged arithmetically. Additionally, the axial forces for different flow field phases of the two splitter blade configurations were compared, with specific time points marked in the Fig. 12. Since the splitter blades are only located within the impeller flow passage, they do not affect the fluid in the gaps between the front and rear shrouds. Consequently, the axial force exerted by the splitter blades has only the flow passage component. From Fig. 10, it is observed that for the original impeller, the axial force shows a periodic variation over time, with a time interval of approximately 0.0069 s between two peak values, equivalent to 20-time steps. Moreover, the axial force is negative, indicating that its direction is toward the impeller inlet. The maximum axial force is approximately 114 N, while the minimum is around 62 N, nearly doubling within a single cycle. This suggests that the impeller is susceptible to fluid excitation forces during pump operation, potentially leading to vibrations. In contrast, for the impellers with splitter blades, the magnitude of the SV-JME . VOL 71 . NO 5-6 . Y 2025 . 163 axial force is significantly increased, but periodic behaviors similar to the original impeller are observed, with the same periodic interval. However, in this case, two peaks – one large and one small – are present within a single cycle, which is distinct from the original impeller. Additionally, the axial force direction for both splitter blade configurations remains the same as that of the original impeller. As the length of the splitter blades decreases, the axial force slightly increases, rising from approximately 149 N to 170 N. However, in terms of relative fluctuation amplitude, the increase in axial force is significantly reduced. This indicates that splitter blades play a role in mitigating axial vibrations of the impeller. 3.4 Analysis of Radial Force of Centrifugal Pump Impeller with Different Splitter Blade Lengths Figure 13 shows the curves of the radial force acting on the impeller varying with time under different design conditions. Meanwhile, the radial force fluctuation curve within one impeller rotation period was selected and its schematic distribution along the circumference of the impeller was plotted, as shown in Fig. 14. It can be seen from the figures that the impellers under the three types of structures all exhibit certain radial force fluctuation characteristics, but the fluctuation characteristics of the original impeller and the impeller with splitter blades are significantly different. For the original impeller, the radial force pulsation curve has a more pronounced periodicity, with only one peak value within one period. However, for the impeller with splitter blades, the periodicity of the radial force is weakened. Comparing impeller #B and impeller #C, it can be seen that the difference in the extremes of the radial force over multiple rotation periods is large, exceeding 40 N, which is much greater than the radial force pulsation value under the design condition. This indicates that the existence of the splitter blades leads to a serious imbalance of the radial forces. From Figure 12, it can also be seen that the radial force on the original impeller shows a central symmetric structure along the rotation axis, but for impeller #B and impeller #C impellers, this pattern is broken, and the radial force acting on the impeller is biased to one side. This requires the rotor bearings of the centrifugal pump to have the ability to withstand larger support forces in order to reduce the radial excursion of the impeller during rotation and maintain relative stability. Therefore, in the actual pump design process, when splitter blades are added to the impeller, the support stiffness and strength of the rotor system bearings need to be further improved. 3.5 Analysis of Pressure Pulsation Characteristics in Impeller with Different Blade Structures To further clarify the infl uence of different impeller structures on the pressure pulsation characteristics within the pump, this section obtained the time-domain and frequency-domain distribution curves a) b) c) Fig. 12. Variation trend of axial force on impeller with time under different impeller structures; a) impeller #A, b) impeller #B, and c) impeller #C of the three pressure monitoring points are shown in Fig. 15. The pulsation coefficient for the original impeller is still slightly higher time-domain diagrams of the pressure pulsation are shown in Fig. than the two diffuser blade impeller structures. This suggests that the 16. To reduce the random error, this paper selected the dimensionless flow stability at the outlet of impeller #B and impeller #C remains pressure coefficient Cp for comparison, which is defined as: better. In the middle of the volute outlet pipe (P3), the pressure pulsation characteristics for the three impeller structures do not show ° N . i 1 U 2 where p is the instantaneous pressure at the monitoring point, [Pa]; U is the peripheral velocity at the impeller outlet, [m/s].  .˜ p / .. .ˆ Cp ˜ , (15) 1 2 From the analysis of Fig. 14, it can be seen that at the impeller inlet (P1), the pressure pulsation curves for the three impeller structures all exhibit a certain periodicity, but within the same period, the original impeller has fewer oscillations than impeller #B and impeller #C, which is related to the number of blades. Although the splitter blades are added downstream of the impeller flow passage, they still have an influence on the pressure fluctuations at the impeller inlet. Additionally, the pulsation amplitude of the original impeller is significantly higher than that of impeller #B and impeller #C. The pulsation coefficient Cp of the original impeller reaches 0.3, while the pressure coefficients of the two diffuser blade impellers are both less than 0.2. This indicates that the flow at the impeller inlet becomes more stable after adding the splitter blades. At the impeller outlet (P2), the periodicity of the pressure pulsation curves still exists for the three impeller structures, but the extreme value of the N significant differences, as this location is far away from the impeller and less influenced by the impeller. To further extract the character istics of pressure pulsation, the time-domain signals at each monitoring point were subjected to fast Fourier transform (FFT) to obtain the frequency-domain distribution shown in Fig. 17. From the figure, it can be seen that at the impeller inlet, the dominant frequency of the original impeller is 145 Hz, which is the blade passing frequency of the original impeller. Furthermore, the pulsation amplitude at the dominant frequency of the original impeller is the highest, followed by the second harmonic, and the third harmonic is further reduced. This phenomenon indicates that the energy at the inlet and outlet of the original impeller gradually attenuates with increasing frequency. For the two impeller structures with splitter blades, their dominant frequency is 300 Hz, which is twice the blade passing frequency of the original impeller. This is due to the addition of the splitter blades. In other words, the blade passing frequency of the splitter blades is twice that of the original impeller blades. In this case, the energy at the dominant frequency of the splitter blades is predominant, while the energy at other harmonics is relatively small. Ignoring the other harmonic energies and only comparing the blade passing frequency and the energy at twice the frequency among the three impellers, it can be found that the amplitudes at 300 Hz are not significantly different, but the original impeller has an additional pulsation energy at 145 Hz. This further indicates that the impeller with splitter blades can reduce the frequency pulsation at the impeller inlet and promote more stable flow. Comparing impeller #B and impeller #C, it can be found that the amplitude at the dominant frequency is slightly lower for impeller #C, which also suggests that the length of the splitter blades is an important factor affecting the pulsation at the impeller inlet. At the impeller outlet, the frequency amplitudes for the three impellers have all decreased. The main frequency characteristics of the impellers with splitter blades are the same as at the impeller inlet, with the blade passing frequency still dominating, but the energy at twice the blade passing frequency is significantly increased. This indicates that the flow structure at the impeller outlet becomes more unstable, which is consistent with the physical reality. At this point, the dominant frequency of the original impeller is still the SV-JME . VOL 71 . NO 5-6 . Y 2025 . 165 a) b) c) Fig. 16. Time domain of monitoring points at designed point; a) point P1, b) P2, and c) P3 blade passing frequency, but the energy at the second, third, and fourth harmonics gradually decreases, which is different from the rapid attenuation of the harmonic energy at the impeller inlet. Under different diffuser blade lengths, the shorter splitter blades have lower dominant frequency energy, and the flow condition is slightly better than the longer splitter blades. At the volute outlet, the frequency characteristics of the three impellers do not differ significantly, with the dominant frequency energy concentrated at low frequencies, which is due to the unsteady flow structure at the volute outlet. In this case, the length of the splitter blades does not have a noticeable impact on the flow field at this location. Therefore, from the pressure pulsation frequency response at different locations, it can be further understood that the splitter blades can reduce the amplitude of the dominant frequency at the impeller inlet, optimize the flow structure at the impeller inlet and outlet, but the influence of the splitter blades on the flow structure within the volute is not significant. a) b) c) Fig. 17. Frequency domain of monitoring points at designed point; a) point P1, b) P2, and c) P3 4 CONCLUSIONS The performance of pumps with splitter blades is significantly different from ordinary impellers, especially in the low flow rate operating range. The head curve of a centrifugal pump with splitter blades has a distinct peak, and there is a positively sloped region at low flow rates. The head is almost 8 % higher under low flow rate while the efficiency of the ordinary impeller structure is higher near the design operating point, but at low flow rates, the straightening effect of the splitter blades gives the impeller with splitter blades a higher efficiency. For the impeller structure with splitter blades, the axial force on the impeller shows a certain periodicity, with the same time interval as the impeller structure without splitter blades. However, within one period, the axial force on the impeller with splitter blades has two peak values – one large and one small. The length of the splitter blades can also affect the magnitude of the axial force, and can help reduce the axial vibration of the impeller. After adding the splitter blades, the original radially symmetric force distribution on the impeller is changed, and the radial force acting on the impeller becomes biased to one side. This requires an increase in the support strength and stiffness of the rotor system. Additionally, the splitter blades can influence the pressure pulsation at the impeller inlet and outlet, optimizing the flow conditions, but their impact on the flow structure within the volute is not significant. This article provides an in-depth analysis of the impact of splitter blade length on the axial and radial forces of a centrifugal pump. The findings are highly significant for optimizing pump design, enhancing operational stability, and extending service life under various working conditions. Future research should further explore the effects of the placement angle, degree of distortion, and number of splitter blades on pump performance to provide more comprehensive design guidance for centrifugal pumps. References [1] Pu, W., Ji, L., Li, W., Shi, W., Tian, F., Xiao, C.et al. Study on the particle dynamic characteristics in a centrifugal pump based on an improved computational fluid dynamics-discrete element model. Phys Fluids 36 123331 (2024) DOI:10.1063/5.0242078. [2] Pu, W., Ji, L., L,i W., Yang, Q., Liu, Z., Yang., Y. et al. Experimental study on the unsteady evolution mechanism of centrifugal pump impeller wake under solid-liquid two-phase conditions: Impact of particle concentration. Phys Fluids 36 113327 (2024) DOI:10.1063/5.0239240. [3] Ji, L., Pu, W., Li, W., Shi, W., Yang, Y. Xiao, C. et al. PIV experimental study on dynamic and static interference flow field of multi-operating centrifugal pump under the influence of impeller wake. Exp Therm Fluid Sci 161 111355 (2024) DOI:10.1016/j.expthermflusci.2024.111355. [4] Li, W., Li, H., Liu, M., Ji, L., Agarwal, R.K., Jin, S. Energy dissipation mechanism of tip-leakage cavitation in mixed-flow pump blades. Phy Fluids 36 015115 (2024) DOI:10.1063/5.0183540. [5] Ji, L., Pu, W., Li, W., Shi, W., Tian, F, Yang, Y. et al. Flow instability in mixed-flow/axial-flow pump: A review of relationship between tip leakage flow distortion and rotating stall. P I Mech Eng A-Pow 239, 225-253 (2024) DOI:10.1177/09576509241287838. [6] Al-Obaidi, A.R., Alhamid, J. Experimental and simulation analyses of the hydraulic complex internal flow characteristics in an axial pump based on varying frequency vibration ranges technique. Int J Interact Des Manuf 19 3661-3681 (2025) DOI:10.1007/s12008-024-02012-9. [7] Al-Obaidi, A.R., Alhamid, J. Investigation of the main flow characteristics mechanism and flow dynamics within an axial flow pump based on different transient load conditions. Iran Sci Technol Trans Mech Eng 47 1397-1415 (2023) DOI:10.1007/s40997-022-00586-x. [8] Al-Obaidi, A.R., Alhamid, J. Analysis of unsteady internal flow characteristics in axial pump with varying number of blades using computational modelling and vibration techniques. Flow Meas Instrum 99 102654 (2024) DOI:10.1016/j. flowmeasinst.2024.102654. [9] Al-Obaidi, A.R. Effect of different guide vane configurations on flow field investigation and performances of an axial pump based on CFD analysis and vibration investigation. Exp Tech 48 69-88 (2024) DOI:10.1007/s40799-023­00641-5. [10] Al-Obaidi, A.R. Experimental diagnostic of cavitation flow in the centrifugal pump under various impeller speeds based on acoustic analysis method. Arch Acoust 48 159-170 (2023) DOI:10.24425/aoa.2023.145234. [11] Giljen, Z., Nedeljkovic, M., Cheng, Y. The Influence of pump-turbine specific speed on hydraulic transient processes. Stroj vestn-J Mech E 70 231-246 (2024) DOI:10.5545/sv-jme.2023.776. [12] Kilavuz, F., Kiral, B.G. Design optimization of mechanical valves in dishwashers based on the minimization of pressure losses. Stroj vestn-J Mech E 70 194-208 (2024) DOI:10.5545/sv-jme.2023.768. [13] Li, W,, Yan,g Q., Yang, Y., Yi, Y., Ji, L., Shi, W. et al. Optimization of pump transient energy characteristics based on response surface optimization model and computational fluid dynamics. Appl Ener 362 123038 (2024) DOI:10.1016/j. apenergy.2024.123038. [14] Namazizadeh, M., Gevari, M., Mojaddam, M., Vajdi, M. optimization of the splitter blade configuration and geometry of a centrifugal pump impeller using design of experiment. J Appl Fluid Mech 1, 89-101 (2020) DOI:10.29252/jafm.13.01.29856. [15] Zhu, D., Xiao, R., Yao, Z., Yang, W., Liu, W. Optimization design for reducing the axial force of a vaned mixed-flow pump. Eng Appl Comp Fluid 14 882-896 (2020) DOI:10.1080/19942060.2020.1749933. [16] Liu, Z., Xu, L., Jia, X., Wu, J., Wang, D. Analysis of liquid flow and axial force calculation in axial clearance for floating impeller of centrifugal pump. Trans Chinese Soc Agri Eng 29 79-85 (2013) DOI:10.3969/j.issn.1002­6819.2013.12.011. [17] Dong, W., Liu, Z., Dong, Y., Yang, Z., Lu, Q., Li, Q. Study on pressure distribution of pump chamber and axial force in particle-laden flow of centrifugal pump. P I Mech Eng A Pow 236 831-839 (2022) DOI:10.1177/09576509211047655. [18] Cao, W., Dai, X., Hu, Q. Effect of impeller reflux balance holes on pressure and axial force of centrifugal pump. J Cent South Univ 22 1695-1706 (2015) DOI:10.1007/s11771-015-2688-2. [19] Cui, B., Li, J., Zhang, C., Zhang, Y. Analysis of radial force and vibration energy in a centrifugal pump. Math Probl Eng 2020 6080942 (2020) DOI:10.1155/2020/6080942. [20] González, J., Parrondo, J., Santolaria, C., Blanco, E. Steady and unsteady radial forces for a centrifugal pump with impeller to tongue gap variation. J Fluid Eng-T ASME 128 454-462 (2006) DOI:10.1115/1.2173294. [21] Ou, M., Shi, W., Jia, W., Fu, Q. Numerical simulation and experimental validation on hydrodynamic radial force of mixed-flow pump impeller. T Chin Soc Agric Eng 31 71-76 (2015) DOI:10.11975/j.issn.1002-6819.2015.09.012. [22] Zhu, L., Jin, Y., Li, Y., Jin, Y., Wang, Y., Zhang, L. Numerical and experimental study on aerodynamic performance of small axial flow fan with splitter blades. J Therm Sci 22 333-339 (2013) DOI:10.1007/s11630-013-0632-z. [23] Liang, Z., Wang, J., Jiang, B., Zhou, B., Yang, W., Ling., J. Large-eddy simulation of flow separation control in low-speed diffuser cascade with splitter blades. Processes 11 3249 (2023) DOI:10.3390/pr11113249. [24] Clark, C., Pullan, G., Curtis, E., Goenaga, F. Secondary flow control in low aspect ratio vanes using splitters. J Turbomach 139 091003 (2017) DOI:10.1115/1.4036190. [25] Yan, P., Chu, N., Wu, D., Cao, L, Yang, S., Wu, P. Computational fluid dynamics-based pump redesign to improve efficiency and decrease unsteady radial forces. J Fluid Eng 139 011101 (2017) DOI:10.1115/1.4034365. [26] Zeng, G., Li, Q., Wu, P., Qian, B, Huanh, B., Li, S. et al. Investigation of the impact of splitter blades on a low specific speed pump for fluid-induced vibration. J Mech Sci Tech 34 2883-2893 (2020) DOI:10.1007/s12206-020-0620-7. [27] Xie, S., Dai, K. Investigation on the added mass effect of centrifugal impellers with splitter blades. 3rd International Conference on Mechanical Design and Simulation 12639 278-285 (2023) DOI:10.1117/12.2682032. [28] Ji, L., Li, W., Shi, W., Agarwal, R.K. Application of Wray-Agarwal turbulence model in flow simulation of a centrifugal pump with semi-spiral suction chamber. J Fluid Eng 143 031203 (2021) DOI:10.1115/1.4049050. [29] Li, W., Ji, L., Li, E., Shi, W., Agarwal, R., Zhou, L.. Numerical investigation of energy loss mechanism of mixed-flow pump under stall condition. Renew Energ 167 740­760 (2021) DOI:10.1016/j.renene.2020.11.146. [30] Ji, L., Li, W., Shi, W., Tian, F., Agarwal, R. Diagnosis of internal energy characteristics of mixed-flow pump within stall region based on entropy production analysis model. Int Commun Heat Mass 117 104784 (2020) DOI:10.1016/j. icheatmasstransfer.2020.104784. [31] Ji, L,, Li, W., Shi, W., Chang, H., Yang, Z. Energy characteristics of mixed-flow pump under different tip clearances based on entropy production analysis. Energy 199 117447 (2020) DOI:10.1016/j.energy.2020.117447. [32] Yu, H., Wang, T., Dong, Y., Gou, Q., Lei, L., Liu, Y. Numerical investigation of splitter blades on the performance of a forward-curved impeller used in a pump as turbine. Ocean Eng 281 114721 (2023) DOI:10.1016/j.oceaneng.2023.114721. [33] Song, H., Zhang, J., Huang, P., Cai, P., Cai, H., Cao, P. Analysis of rotor-stator interaction of a pump-turbine with splitter blades in a pump mode. Mathematics 8 1465 (2020) DOI:10.3390/math8091465. [34] Siddique, M., Samad, A., Hossain, S. Centrifugal pump performance enhancement: Effect of splitter blade and optimization. P I Mech E A-J Pow 236 391-402 (2022) DOI:10.1177/09576509211037407. [35] Khalafallah, M.G., Saleh, H.S., Ali, S.M., Abdelkhalek, H.M. CFD investigation of flow through a centrifugal compressor diffuser with splitter blades. J Eng Appl Sci 68 43 (2021) DOI:10.1186/s44147-021-00040-w. [36] Ji, L., Li, W., Shi, W., Tian, F., Agarwal, R. Effect of blade thickness on rotating stall of mixed-flow pump using entropy generation analysis. Energy 236 121381 (2021) DOI:10.1016/j.energy.2021.121381. [37] Ji, L., He, S., Li, Y., Li, W., Shi, W., Li, S. et al. Investigation of energy loss mechanism of shroud region in A mixed-flow pump under stall conditions. P I Mech E-A Pow 237 1235-1250 (2023) DOI:10.1177/09576509231162165. [38] Li,W., Huang, Y., Ji, L., Ma, L, Agarwal, R.K., Awais, M. Prediction model for energy conversion characteristics during transient processes in a mixed-flow pump. Energy 271 127082 (2023) DOI:10.1016/j.energy.2023.127082. [39] Li, W., Li, S., Ji, L., Li, E., Shi, W., Agarwal, R. et al. Effect of inlet elbow on rotation stall in waterjet propulsion pump. Fund Res 4 898-906 (2024) DOI:10.1016/j. fmre.2022.05.029. [40] Li, W., Liu, M., Ji, L., Wang, Y., Awais, M., Hu, J. et al. Research on the matching characteristics of the impellers and guide vanes of seawater desalination pumps with high capacity and pressure. J Marine Sci Eng 10 115 (2022) DOI:10.3390/ jmse10010115. [41] Ji, L., Li, W., Shi, W., Tian, F., Li, S., Agarwal, R. Influence of different blade numbers on the performance of “saddle zone” in a mixed flow pump. P I Mech Eng A-J Pow 236 477-489 (2022) DOI:10.1177/09576509211046074. [42] Zhu, Z., Lin, Y., Li, X., Zhai, L, Lin, T. Axial thrust instability analysis and estimation theory of high speed centrifugal pump. Phys Fluids 34 075118 (2022) DOI:10.1063/5.0098194. Acknowledgments The work was sponsored by the Science and Technology Plan Project of Chuzhou (2024CX007) and Research Initiation Fund Project of Chuzhou University (2023qd21). Received: 2024-12-28, revised: 2025-03-27, 2025-04-02, accepted: 2025-04-09 as Original Scientific Paper. Conflict of interest The authors declared that they have no conflicts of interest to this work. Data Availability The data supporting the findings of this study are included in the article. Author Contribution Qijiang Ma: Writing– original draft, Writing – review & editing , Software, Conceptualization. Zhenbo Liu: Writing– original draft, Writing – review & editing, Supervision. Sen Jiang: Writing– original draft, Writing – review & editing. Raziskava vpliva dolžine vmesnih lopatic na radialne in aksialne sile v centrifugalni crpalki Povzetek Za izboljšanje obratovalne stabilnosti centrifugalnih crpalk ta študija preucuje vpliv dolžine vmesnih lopatic na aksialne in radialne sile v centrifugalnih crpalkah. S pomocjo SST k-. turbulencnega modela in eksperimentalne analize smo primerjali karakteristiko, aksialno silo, radialno silo ter casovno-frekvencne znacilnosti pulziranja tlaka med rotorji brez vmesnih lopatic in rotorji z vmesnimi lopaticami dveh razlicnih dolžin. Rezultati kažejo, da prvotni rotorji dosegajo višjo ucinkovitost v bližini tocke najvecjega izkoristka. Vendar pa pod pogoji nizkega pretoka usmerjanje toka z vmesnimi lopaticami omogoca, da rotorji z vmesnimi lopaticami dosegajo višji izkoristek. Pri rotorjih z vmesnimi lopaticami aksialna sila kaže periodicnost z dvema vrhovoma – enim vecjim in enim manjšim – znotraj vsakega cikla. Dodatek vmesnih lopatic povzroci spremembo radialne sile na rotor, kar zahteva povecanje nosilnosti in togosti rotorskega sistema. Poleg tega vmesne lopatice vplivajo na tlacne pulzacije na vstopu in izstopu iz rotorja; prvotni rotor izkazuje tlacne pulzacije pri frekvenci 145 Hz, kar pa je pri rotorjih z vmesnimi lopaticami drugace. Tako lahko rotorji z vmesnimi lopaticami zmanjšajo frekvencne pulzacije in izboljšajo pretocne pogoje. Študija daje pomemben teoreticni vpogled in podatkovno podporo za hidravlicno konstrukcijsko optimizacijo centrifugalnih crpalk. Kljucne besede centrifugalna crpalka, aksialna sila, radialna sila, dolžina vmesnih lopatic, tlacne pulzacije Multi-Objective Optimization Design of the Ejector Plate for Rear-Loader Garbage Trucks Fu-sheng Ding1 – Hong-ming Lyu1 – Jun Chen2 – Hao-ran Cao1 – Lan-xiang Zhang1 1 Yancheng Institute of Technology, School of Automotive Engineering, China 2 Jiangsu Yueda Special Vehicle Co., Ltd., China hongming_lv@ycit.edu.cn Abstract This work presents a multi-objective optimization design approach for the ejector plate, a critical component of rear-loader garbage trucks, with the goal of ensuring structural integrity and optimizing lightweight performance. A parametric finite element model of the ejector plate is developed with optimization objectives focused on minimizing mass, controlling deformation, and reducing the maximum von Mises stress. Through sensitivity analysis, seven key variables are selected for optimization. A Box-Behnken design (BBD) is used to systematically explore these parameters, and a Kriging surrogate model is constructed to approximate the objective function, with accuracy benchmarked against response surface methodology (RSM). The Non-dominated Sorting Genetic Algorithm II (NSGA-II) is applied to derive the optimal solution, achieving a lightweight design meeting all structural requirements. The results show that the mass of the ejector plate of the rear-loading waste compactor can be reduced by 6.06 % through structural optimization, while meeting the strength and deformation criteria. This improvement not only enhances waste transportation efficiency, but also lowers production costs and enhances material utilization. Keywords garbage truck, ejector plate, multi-objective optimization, NSGA-II, Kriging Highlights . Optimizing the ejector plate considering the mass, maximum displacement and maximum von Mises stress of the ejector plate. . Using sensitivity analysis to identify key design variables of the ejector plate. . Using Kriging method to construct a highly accurate surrogate model. . Applying the non-dominated sorting genetic algorithm II for multi-objective optimization of the ejector plate. 1 INTRODUCTION As the global economy grows and urbanization accelerates, the generation of household waste has increased significantly. This trend is putting immense pressure on urban environmental health and poses challenges for sustainable urban development [1]. Given the substantial volume of urban waste, efficient transportation and disposal methods have become critical concerns for municipalities and related departments [2]. Rear loader garbage trucks are specialized vehicles designed to operate with waste compaction transfer stations. Their widespread adoption by sanitation and municipal departments can be attributed to their large loading capacity and effective sealing. The ejector plate, a key component of these trucks, plays a vital role in the loading and unloading of waste. Different types of waste impose varying demands on the strength of the ejector plate. During the compaction and loading/unloading processes, the ejector plate is subjected to forces from both the hydraulic cylinder and the waste itself, resulting in varying degrees of deformation. Such deformation can reduce the compaction ratio within the truck's box body, ultimately diminishing loading and unloading efficiency and affecting the truck's overall waste capacity. Therefore, research into the ejector plate is crucial for the effective design of garbage trucks. The three dimensional (3D) model of the ejector plate was developed in SolidWorks software, which was then subjected to finite element analysis (FEA) to evaluate its stress and displacement distributions under operational conditions [3]. The garbage truck loading mechanism served as the research focus, with parametric analysis conducted to establish the functional relationships among compression-filling force, pusher stroke, and installation angle [4]. The garbage truck manipulator was digitally prototyped in SolidWorks for performance enhancement, followed by a multi-domain simulation workflow: multi-body dynamic (MBD) analysis in ADAMS determined operational load spectra under various working conditions, FEA in HyperWorks assessed structural integrity, and topology optimization integrated with genetic algorithms achieved 16.73 % mass reduction while maintaining mechanical performance [5]. Lightweighting vehicles is an essential strategy for saving energy and reducing emissions [6,7]. Reducing the weight of vehicles can reduce energy consumption while improving dynamics and braking performance [8,9]. The main approaches to lightweight design include structural optimization, process optimization and material lightweighting [10-12]. Structural optimization can be further categorized into size optimization, shape optimization, and topology optimization [13]. Currently, size and shape optimization are widely used in engineering applications [14,15]. Size optimization can be divided into discrete and continuous categories based on the continuity of design variables. Continuous size optimization results typically require rounding to fit available size parameters, making discrete size optimization more suitable for practical project needs [16]. Furthermore, discrete optimization allows for the simultaneous optimization of multiple variables, potentially yielding superior results. In the realm of optimal design algorithms, Goldberg was the first to apply genetic algorithms (GA) to multi-objective optimization in 1989 [17]. Subsequent research has built upon this foundation. For instance, Paz et al. [18] utilized multivariate adaptive regression spline techniques combined with multi-objective genetic algorithms to enhance the energy absorption properties of automotive components while reducing mass. Velea et al. [19] conducted multi-objective optimization of a composite body for electric vehicles, considering weight, cost, and stiffness. Duan et al. [20] applied a multi-objective particle swarm optimization algorithm for lightweight design of body-in-white structures to meet reliability requirements. Jiang et al [21] employed the Kriging surrogate model alongside the non-dominated sorting genetic algorithm II (NSGA-II) for the multi-objective optimal design of suspension arms and torsion beams, achieving weight reduction without compromising structural or vehicle performance. Xie et al. [22] implemented a multi-objective optimization method incorporating the TOPSIS method for the lightweight design of commercial vehicle cabs, successfully meeting design requirements while achieving weight reduction. The lightweight design of the ejector plate is directly related to the load capacity and reliability of the garbage truck. Therefore, pursuing a lightweight design for the ejector plate while considering both structural performance and load capacity is essential. This paper focuses on the ejector plate of a specific garbage truck, utilizing 3D design software for parametric modeling and conducting finite element analysis based on actual working conditions. A sensitivity analysis of mass, displacement, and equivalent force concerning fourteen key parameters is performed, identifying 7 parameters with greater sensitivity as design variables. These variables are then optimized using a combination of the Kriging surrogate model and the NSGA-II algorithm to achieve the lightweight design of the ejector plate. 2 METHODS AND MATERIALS 2.1 Working Principle of the Ejector Plate As shown in Fig. 1, the rear loader garbage truck consists primarily of several components, including the vehicle chassis, box body, hydraulic cylinder, packing mechanism, and ejector plate, as illustrated in Fig. 2. The box body is mounted on and rigidly connected to the chassis frame. At the rear of the truck, the packing mechanism is responsible for loading and initially compacting the waste. The waste is then pushed into the box body by the action of a scraper on the packing mechanism. After each compaction, both the compacted waste and the newly loaded waste are gradually pushed toward the rear of the box body, with the ejector plate continuously moving back as new waste is added. Fig. 1. Construction of the rear-loader garbage truck; 1 chassis, 2 cylinder 3 ejector plate, 4 box body, and 5 precompressor 2.2 Design Requirements for the Ejector Plate The ejector plate consists of the front plate, frame skeleton, guide rail skeleton, cylinder support, and other components. As a critical load-bearing element of rear-loader garbage trucks, the mechanical properties of the ejector plate significantly impact the overall performance of the vehicle. During operation, the ejector plate experiences the combined effects of hydraulic cylinder thrust, friction between the ejector plate and the guide rail, and the extrusion pressure from the refuse. The deformation of the ejector plate skeleton directly affects the gap between the ejector plate and the cargo area, which, in turn, influences the complete discharge of the waste. Therefore, the design of the ejector plate must ensure that the deformation of the frame skeleton remains below 10 mm. In addition, under maximum external load conditions, the stress on all components should not exceed 355 MPa. To meet the performance requirements, the design of the ejector plate should also prioritize minimizing weight and maintaining high quality. Fig. 2. 3D model of the ejector plate 2.3 Optimization Flow for the Ejector plate This study employs a combination of the finite element method (FEM), Kriging Surrogate Model, and GA for optimization. The finite element method is utilized to calculate the maximum working load conditions of the ejector plate. Based on the results of the finite element analysis, a lightweight multi-objective structural optimization program is proposed which focuses on minimizing the total deformation, mass and von Mises stress. The Kriging surrogate model is known for its high accuracy, adaptability and scalability, making it widely applicable in the structural optimization design of various mechanical components [23,24]. Genetic algorithms are favored for multi-objective optimization due to their broad adaptability, global search capability, parallelism, and independence from derivative information [25]. The optimization process is illustrated in Fig. 3. The specific steps are as follows: • Parametric modeling: The ejector plate is modeled parametrically in 3D modeling software, and its finite element model is created in HyperWorks software, followed by static analysis. • Sensitivity analysis: 14 parameters within the ejector plate skel­eton are analyzed for sensitivity, leading to the selection of 7 pa­rameters that significantly influence the mass, total displacement, and von Mises stress of the ejector plate skeleton for further anal­ysis. • Experimental design: The 7 highly sensitive parameters are sam­pled and tested using the Box-Behnken experimental design meth­od. • Surrogate model generation: A Kriging surrogate model is devel­oped based on the experimental design parameters. • Optimization: A multi-objective genetic algorithm (NSGA-II) is employed to obtain the optimal set of solutions and to validate the appropriateness of the selected optimization parameters. 2.4 Finite Element Analysis of the Ejector plate 2.4.1 Finite Element Model When modeling, capturing the essence of the research object and appropriately simplifying the original model are essential for improving simulation efficiency and quality. The rubber buffer blocks on the ejector plate, and the sealing plates on both sides serve auxiliary roles that have minimal impact on the forces in the model, allowing them to be omitted from the modeling process. The focus is on the overall structural characteristics of the ejector plate rather than localized welding issues. Rigid units are employed to simulate the weld relationships between different components, facilitating the transfer of forces. This simplified approach enhances the fidelity of the finite element results while significantly improving solution speed, enabling a more efficient design process that can accommodate a large number of experimental simulations. In the ejector plate assembly, the cylinder support is a casting and is simulated using 8 nodes hexahedron elements. The other components are sheet metal parts, characterized by smaller dimensions in the thickness direction and larger dimensions in length and width, and are simulated using 4 nodes shell elements. Considering the enterprise's specifications for mesh size, as well as the need to maintain computational accuracy while minimizing computational cost, an 8 mm mesh was adopted for the model. Upon completing the meshing of the ejector plate, the model consists of 117,823 nodes and 117,117 elements, as shown in Fig. 4. The material of the ejector plate is high-quality structural carbon steel (Q355B), with properties including a Young's modulus of 2.1 × 105 MPa, a density of 7850 kg/m³, and a Poisson's ratio of 0.3, and a minimum yield strength of 355 MPa. As one of the primary load-bearing components of the truck, the ejector plate is subjected to the thrust from the hydraulic cylinder, the force exerted by the waste on its front plate during loading, and the support force from the guide rail. Due to the considerable stroke length of the ejector plate, a multi-stage hydraulic cylinder is employed; the shorter the stroke of the hydraulic cylinder, the greater the thrust produced. According to the hydraulic cylinder's operating behavior, the ejector plate experiences the greatest load and deformation when the truck is nearly full of waste and the ejector plate is close to its innermost position within the box body. The hydraulic system operates at a pressure of 17.6 MPa. The ejector plate cylinder is a 3-stages rod cylinder, with the diameters of the rods measuring 90 mm (1st stage), 70 mm (2nd stage), and 50 mm (3rd stage). The maximum thrusts exerted by these rod sections are calculated to be 111.9 kN, 67.7 kN, and 34.5 kN, respectively. When waste is loaded into the box body, the force exerted by the waste on the ejector plate must overcome both the hydraulic cylinder's force and the friction between the ejector plate and the guide rail before the ejector plate can move inward. Thus, the ejector plate experiences maximum load when the hydraulic cylinder's thrust is 111.9 kN. The ejector plate moves slowly within the box body, neglecting the effects of dynamic loading. The friction force, which cannot be directly measured, is assumed to be approximately 12 kN (about 10 % of the maximum thrust). During slow movement, the force from the waste on the front plate is balanced by the hydraulic cylinder's force and N The NSGA-II algorithms is used to obtain the Pareto optimal solution set Y Fig. 3. Flow chart for optimization of the ejector plate structure the friction between the ejector plate and the guide rail, resulting in a Taking into account the working conditions of the ejector plate, pressure of 0.06 MPa acting on the front plate. overload protection, personnel protection, material properties, etc., a Fig. 6. Distribution of the load on ejector plate In accordance with actual working conditions, fixed displacement constraints in the y-direction (forward and backward) are applied at the cylinder support of the ejector plate. Displacement constraints in the x-direction (left and right) are applied to the lateral friction blocks of the ejector plate guide, and z-direction (vertical) displacement constraints are applied to the upper and lower friction blocks of the guide, as shown in Fig. 5. A uniform compressive pressure load of 0.06 MPa is applied to the side of the front plate in contact with the waste. Considering that the waste cannot be loaded to the top of the box body, this load is applied only up to a height of 960 mm from the bottom plate, as illustrated in Fig. 6. safety factor of 1 is selected, and the allowable stress of the material is 355 MPa. 2.4.2 FEA Result Analysis The finite element analysis yields to the static deformation and von Mises stress contours for the ejector plate and its skeleton, as illustrated in Fig. 7 through Fig. 10. Fig. 8. Static deformation of the ejector plate skeleton As shown in Fig. 7, the maximum deformation of the ejector plate is approximately 21.4 mm, which occurs in the central and upper areas of the front plate. This deformation corresponds to the actual working conditions. The ejector plate has a width of 1800 mm and the front plate is made of a 2.3mm thick steel plate. Although this design results in a relatively low stiffness, the significant deformation of the front plate does not affect its performance, provided that the deformation of the ram skeleton remains within acceptable limits. Figure 8 indicates that the maximum deformation of the ejector plate skeleton is approximately 7.36 mm, concentrated in the middle of the upper tube, with progressively smaller deformations further down. This deformation is within the allowable limit of 10 mm specified in the design, demonstrating that the current ejector plate meets the stiffness requirements and possesses a degree of stiffness redundancy. Fig. 9. von Mises stress distribution in the ejector plate From Figure 9, it can be seen that the maximum von Mises stress of the ejector plate is 304.7 MPa, which is located in the middle and upper sections of the front plate. In Figure 10, the von Mises stress contour for the ejector plate skeleton shows a maximum stress of 292 MPa, which is located in the lower sections of the side tube and below the design allowable stress limit of 355MPa. To address the potential influence of mesh size on the simulation results, a systematic mesh sensitivity study was carried out specifically on the side tube (identified as the critical region with the highest stress concentration). Three different mesh configurations were tested: 8 mm, 6 mm, and 4 mm average element size. The corresponding results are summarized in Table 1 below. Table 1. Mesh sensitivity analysis for the side tube Mesh size [mm] 8 6 4 Max stress [MPa] 291.7 295.5 296.9 Stress increment [MPa] - +3.8 +5.2 Maximum stress difference - +1.3 % +1.8 % Number of side tube elements 2760 4450 9905 From Table 1, it can be seen that as the mesh is refined from 8 mm to 4 mm, the maximum stress exhibits a monotonic but diminishing increase, with a total variation of 1.8 % (5.2 MPa absolute difference), and the number of side tube elements shows an increase in geometric multiples, with a total variation of 258.9 % (7145 elements difference). The relative stress change between successive refinements (6 mm vs. 8 mm: 1.3 %; 4 mm vs. 6 mm: 0.5 %) demonstrates convergence behavior, indicating that further mesh refinement beyond 6mm yields marginal improvements in accuracy. The sub-2% discrepancy between the 8 mm and 4 mm meshes falls within typical engineering tolerance thresholds for such analyses [26], confirming that the 8 mm mesh provides sufficient accuracy while maintaining computational efficiency. The results of finite element analysis indicates that the current ejector plate skeleton meets the strength requirements and has a significant material surplus, allowing for potential lightweighting studies of the ejector plate structure. 2.5 Sensitivity analysis 2.5.1 Determination of Design Variables and objective functions According to the structural characteristics and the stress distribution of the ejector plate’s skeleton , the upper tube thickness P1, reinforced tube . thickness P2, side tube thickness P3, side steel tube thickness P4, angle iron thickness P5, bracket thickness P6, bottom skeleton thickness P7, reinforced tube . thickness P8, reinforced steel tube thickness P9, angle brace thickness P10, bottom plate thickness P11, rail angle brace plate thickness P12, rail skeleton thickness P13, and rail support plate thickness P14 are used as the candidate design variables. The maximum static deformation D, maximum von Mises stress S, and mass M are used as objective functions. The position of each design variable of the ejector plate’s skeleton is shown in Fig. 11. 2.5.2 Sensitivity analysis of design variables Sensitivity analysis is a widely used tool in structural optimization, aimed at identifying design variables that significantly impact structural performance. The sensitivity value associated with changes in a design variable intuitively indicates both the magnitude and direction of its effect on performance, allowing for the rapid screening of critical design variables for optimization [27]. This analysis typically examines how the output response of a system is affected by its input parameters. For complex mathematical models, determining higher-order sensitivities can be challenging. Therefore, conventional sensitivity analysis focuses on the first-order partial derivatives of the design response with respect to the design variables. The mathematical expression is given by: ˜ gX .gX (). (). ˆ. (1) sen x .x , ° i . i where g(X) is the system performance index, X is the system design parameter vector, and xi is the ith parameter in the system design parameter vector. The sensitivity of the 14 design variables in the ejector plate structure corresponding to the three performance indicators of mass, maximum von Mises stress and total deformation are calculated, as shown in Fig. 12, respectively. From Figure 12, we can see that: 1. The significance of all variables for total mass is positively correlated, with P7, P9, and P12 having the greatest significance for total mass. 2. P3, P12, P13, and P14 exhibit the most significant negative correlation with maximum von Mises stress, while the other variables show less significant correlations. 3. All variables show a negative correlation with total displacement. P1, P2, P3, P4, P7, P11 and P12 are most significantly negatively correlated with maximum deformation, while the remaining variables have a lesser impact. Considering the sensitivity of each design parameter to the response value, 7 parameters are selected as design variables for the final size optimization. The initial values of each parameter and their allowable range of variation are presented in Table 2. Table 2. Initial values and range of the ejector plate design variables [mm] Design variables Parameters Original Value range Upper tube P1 3 2~4 Reinforced tube I P2 3.2 2~4.4 Side tube P3 3.2 2~4.4 Side steel tube P4 3.2 2~4.4 Bottom skeleton P7 5 4~6 Guide rail angle brace plate P12 3.2 2~4.4 Guide rail skeleton P13 4 3~5 2.6 Surrogate Model A surrogate model is a mathematical representation constructed from a finite set of data obtained through computational or physical experiments. Complex models that require significant computational resources can be approximated using finite element-based models or less computationally intensive statistical models as surrogate models [28]. In optimization contexts, where evaluating objective and/or constraint functions demands considerable computational effort, a surrogate model effectively replace these functions [29]. The original model can be accurately predicted only if a highly precise surrogate model is developed. The design of experiments (DOE) is a crucial step in creating a surrogate model, utilizing mathematical analyses such as probability theory and linear generation to identify reasonable discrete sample points. The distribution of these sample points within the design space significantly influences the construction of the Kriging surrogate model, making the selection of an appropriate experimental design method critical for subsequent optimization. Commonly used experimental design methods include central composite design (CCD), Box-Behnken design (BBD), Latin hypercube design (LHS), and optimal space filling design (OSF). For this study, the Box-Behnken design method was chosen due to its advantages of requiring fewer trials and not necessitating the measurement of vertex points compared to other design methods. 2.6.1 Box-Behnken Design (BBD) Using the upper tube thickness P1, the reinforced tube . thickness P2, the side tube thickness P3, the side steel tube thickness P4, the bottom skeleton thickness P7, the rail angle brace thickness P12, and the rail skeleton thickness P13 as independent variables, and the total mass (M) of the ejector plate, the maximum displacement (D), and the maximum von Mises stress (S) as the evaluation indices, a total of 62 experimental samples and their corresponding response values were extracted, as shown in Table 3. Table 3. Test samples and response values for BBD tests P1 P2 P3 P4 P7 P12 P13 MD S Tests [mm] [mm] [mm] [mm] [mm] [mm] [mm] [kg] [mm] [MPa] 1 3 3.2 3.2 2 4 2 4 169.4 7.856 414.3 2 3 3.2 3.2 4.4 4 2 4 174.8 7.509 418.4 3 3 3.2 3.2 2 6 2 4 178.9 7.685 415.3 4 3 3.2 3.2 4.4 6 2 4 184.2 7.341 419.4 5 3 3.2 3.2 2 4 4.4 4 175.7 7.530 224.5 6 3 3.2 3.2 4.4 4 4.4 4 181.0 7.195 226.6 7 3 3.2 3.2 2 6 4.4 4 185.1 7.386 224.4 8 3 3.2 3.2 4.4 6 4.4 4 190.5 7.052 227.0 … …………………… … … 60 3 3.2 3.2 3.2 5 3.2 4 180.0 7.36 291.7 61 3 3.2 3.2 3.2 5 3.2 4 180.0 7.36 291.7 62 3 3.2 3.2 3.2 5 3.2 4 180.0 7.36 291.7 2.6.2 Kriging Surrogate Model Kriging refers to a surrogate model based on Gaussian process modeling, which first originated in a geostatistical paper by Krige [30], and is now one of the most widely used surrogate modeling methods. The Kriging surrogate model can be expressed as follows: yx()˜ fx()° Zx(), (2) where y(x) is an unknown function of the optimization objective, f(x) is a polynomial function on the variable x, and Z(x) is a stochastic function mainly used to correct the error of the global model. The expression for f(x) is: kk k.1 k 2 fx()˜.0 °..iix°..ii ix°...ijxxij, (3) i˜1 i˜1 i˜1ji˜°1 where k is the number of design variables, in this paper k=7,ß is the coefficient to be determined for each equation. Z(x) is a random function that follows a Gaussian normal distribution with zero mean. EZx() , ˜ °. 0 (4) VZx˜ ()°.. 2. (5) Z(x) creates a ‘local’ deviation that satisfies Eq. (4) with minimum variance to obtain the best valuation of the response surface of the kriging model interpolated to the N sample data points. The covariance matrix of is given by the following equation: CovZxZx˜ () .. .. 2 ˆ[(, )]., (6) ,( ) R rxx ° ij ij where R is the N×N symmetric positive definite correlation matrix; r(xi, xj) is a Gaussian correlation function, which can represent the spatial correlation between any two sample points xi and xj. The expression is as follows: . ° M ° 2 ˆ The fitting accuracy of the Kriging surrogate model and response surface methodology (RSM) surrogate model are obtained by solving Eq. (10) and Eq. (11), as shown in Fig. 13, and the specific values are shown in Table 4. It can be seen that the fitting accuracy values between the three objective functions of the ejector plate and the 7 key design variables are all greater than 0.90, while the coefficient of determination R2 of the general engineering requirements for the k 1 where .k is the unknown parameter used for fitting; M is the number ˜  (7) rx(,)e ˜ x xx .. .. response surface model should not be less than 0.9 [31]. xp , j k ki kj Table 4 shows the specific index values of the accuracy test results of each surrogate model. Table 4. Accuracy of objective function of design variables; xki and xkj are the components of the kth sampling point xi, xi. Based on the unbiasedness of Z(x) and the minimum variance of the estimate, the relevant parameter .k given by the maximum possible estimate, i.e. when .> 0, is derived to maximize the following equation: .2 n ln()..ln R()x A ˜° l , (8) 2 .2 where nl is the number of response values, s is the variance estimate, and |R(x)| is the correlation value between the point to be measured and the sampling point. In other words, turn to find the minimum value. min( x) R( x) m . 2. (9) .° 1 x˜0 2.6.3 Surrogate Model Accuracy Validation The surrogate model may have some errors in the process of establishment, in order to verify the accuracy of the established surrogate model, the deterministic coefficient R2, root mean square error (RMSE) are introduced. The value of R2 ranges between [0,1], when the value of R2 tends to be closer to 1 and the value of RMSE tends to be closer to 1, they indicate that the accuracy of the fitted surrogate model is higher, and vice versa is lower. The formulas for calculating are as follows:  n ..yi °y  i ˆ. 2 2 i˜1 .. ˜° (10) 2 R 1  n yi °yi  , i˜1 Kriging surrogate model RSM surrogate model Performance R2 RMSE R2 RMSE mass 0.9505 0.0538 0.9468 0.0583 displacement 0.9659 0.0411 0.9503 0.0623 stress 0.9327 0.0601 0.9031 0.0736 The results show that both Kriging model and RSM model exhibit high accuracy, with R² values exceeding 0.9 and RMSE values falling below 0.08. Overall, the Kriging model slightly exceeds the RSM model in terms of accuracy. Therefore, the Kriging surrogate model will be used for the subsequent multi-objective optimization. 2.7 Multi-Objective Optimization of the Ejector Plate 2.7.1 Mathematical Modeling In order to achieve the purpose of lightweighting of the ejector plate, according to the Kriging surrogate model established by fitting, the mass of the ejector plate M, the maximum static deformation of the ejector plate skeleton D and the maximum von Mises stress S are taken as the optimization objective function, the maximum static deformation of the ejector plate skeleton D and the maximum von Mises stress S are taken as the constraints, and a total of 7 key dimensions of the ejector plate obtained from screening P1~P4, P7, P12 and P13 are taken as the independent variables of the objective function, so that the following multi-objective optimization mathematical model can be established. The multi-objective optimization mathematical model of ejector plate is as follows: Find xi (i = 1~4,7,12,13) 2 N min: {M(x)} max: {D(x), S(x)} , (12)  ... ˆ.. 1 RMSE ˜ ° (11) y y i , i N subject to D(x)=10mm, where yi are true values of the test point, .yi predicted values of the S(x)=355 MPa, surrogate model for the test point, yi average of true values, and N number of sample points. ximin =xi =ximax Fig.13. Surrogate model accuracy validation; a) mass b) displacement c) Von Mises stress where x is the set of design variables; M(x) is the mass of the ejector plate [kg]; D(x) is the displacement of the ejector plate skeleton [mm]; S(x) is the von Mises stress of the ejector plate skeleton [MPa]; ximin and ximax are the lower and upper bounds of each design variable. 2.7.2 Multi-Objective Optimization Using NSGA-II NSGA-II is one of the most widely used and effective multi-objective evolutionary algorithms, originally proposed by Deb [25]. Unlike traditional genetic algorithms, NSGA-II introduces a fast non-dominated sorting method, an elite maintenance strategy, and an efficient congestion distance estimation process. These enhancements significantly accelerate iterative convergence, reduce computational complexity, and ensure population diversity [32]. In this paper, the NSGA-II algorithm is employed to address the multi-objective optimization problem. The population size is set to 120, with a maximum of 200 iterations. The algorithm utilizes a crossover probability of 0.9, a crossover distribution exponent of 10.0, and a mutation distribution exponent of 20.0. 3 RESULTS AND DISCUSSION After 240,001 iterations, 742 non-dominated solutions were obtained. The Pareto front for the multi-objective optimization of the lightweight design of the ejector plate, produced by the NSGA­II algorithm, is shown in Fig. 14. In multi-objective optimization problems, the global optimal solution is typically not unique; instead, it consists of multiple optimal solutions, collectively known as the Pareto optimal solution set. After solving the developed response surface model, three sets of candidate solutions were identified, as presented in Table 5. Table 5 displays the optimized design parameters suggested by the software platform after problem resolution. Considering its superior performance in minimizing mass, displacement, and stress, Scheme 3 is identified as the most favorable design for the ejector plate. The Table 5. Candidate combinations for design parameter optimization optimal design variables obtained after optimization were rounded for correction, and the finite element model of the ejector plate was re-analyzed for static characteristics based on this set of parameters. The performance indices corresponding to the design sizes of the ejector plate before and after optimization are summarized in Table 6. Table 6. Comparison of the design variables before and after optimization Design variables Original Optimal Variation Error P1 [mm] 3.0 2.3 –0.7 P2 [mm] 3.2 2.5 –0.7 P3 [mm] 3.2 3.0 –0.2 P4 [mm] 3.2 2.5 –0.7 P7 [mm] 5.0 4.0 –1.0 P12 [mm] 3.2 3.2 0 P13 [mm] 4 3.5 –0.5 M [kg] 180.0 169.1 –10.9 6.06 % D [mm] 7.36 8.20 +0.84 11.4 % S [MPa] 291.7 345.8 +54.1 18.5 % As seen from Table 6, it demonstrates that the maximum displacement of the ejector plate skeleton has increased by 0.84 mm from 7.36 mm before optimization to 8.20 mm after optimization and the maximum equivalent stress of the ejector plate skeleton has increased by 54.1 MPa from 291.7 MPa before optimization to 345.8 MPa after optimization. Although both indicators have increased, they meet the design requirements. After optimization, the ejector plate's total mass decreased from 180 kg to 169.1 kg, achieving a reduction of 10.9 kg (6.06 %) compared to the original design. The optimization effect is clear, achieving the intended optimization goals. The effectiveness of the proposed scheme has been validated in actual production, as illustrated in Fig. 15. The mass of the ejector plate has been reduced, material utilization has improved, and safety is maintained. No P1 [mm] P2 [mm] Optimized design parameters P3 [mm] P4 [mm] P7 [mm] P12 [mm] P13 [mm] M [kg] Optimal results D [mm] S [MPa] 1 2.4278 2.5407 2.6980 2.5852 4.1747 3.1341 3.5848 167.1 8.901 355.0 2 2.4408 2.5404 2.7123 2.5228 4.1790 3.1284 3.4678 167.0 8.862 355.0 3 2.4231 2.4693 2.7431 2.5534 4.1815 3.1013 3.5779 167.0 8.853 355.0 4 CONCLUSIONS This work presents an innovative lightweight design framework for rear-loader garbage truck ejector plates, which integrates multi-objective optimization algorithms with surrogate modeling techniques to achieve balanced improvements in structural performance and economic efficiency. Through systematic optimization, the final design demonstrates a total mass reduction to 169.1 kg (6.06 % lighter than the original design) while satisfying all operational constraints. The main findings can be summarized as follows: • Comparative analysis of surrogate modeling approaches shows that the Kriging model exhibits superior prediction accuracy com­pared to the RSM model, particularly in capturing the response characteristics of the ejector plate system. • The implementation of global sensitivity analysis allows effective identification of critical design parameters that have dominant in­fluences on key performance indicators, thereby reducing compu­tational costs by reducing experimental design parameters. • The NSGA-II evolutionary algorithm is shown to be effective in generating Pareto optimal solutions for push plate parameters, achieving convergence within 200 generations while maintaining solution diversity through tournament selection and simulated bi­nary crossover mechanisms. In particular, the developed methodology demonstrates strong scalability and provides a generalized framework for the optimization of various load-bearing components in heavy-duty vehicles, including but not limited to chassis structures and actuator systems. Future work should focus on experimental validation of optimized designs under real-world operating conditions and extension. References [1] Mostafayi Darmian, S., Moazzeni, S., Hvattum, L.M. Multi-objective sustainable location-districting for the collection of municipal solid waste: Two case studies. Comput Ind Eng 150 106965 (2020) DOI:10.1016/j.cie.2020.106965. [2] Wan, H., Ma, J., Yu, Q., Sun, G., He, H., Li, H. Modeling and optimization of multi-model waste vehicle routing problem based on the time window. J Database Manag 34 1-16 (2023) DOI:10.4018/jdm.321543. [3] Negreanu, G.P., Voicu, G., Lazea, M., Constantin, G.-A., Stefan, E.-M., Munteanu, M.-G., et al. Finite element analysis of the compaction plate from a garbage truck. E3S Web Conf 180 04006 (2020) DOI:10.1051/e3sconf/202018004006. [4] Zheng, X., Guo, J., Li, S. Design and mechanical analysis of the loading mechanism of the back-loading garbage truck. IOP Conf Ser Mater Sci Eng 612 032043 (2019) DOI:10.1088/1757-899x/612/3/032043. [5] Luo, L., Li, X., Yuan, M. (2023). Lightweight optimization of garbage collection truck manipulator based on parametric modelling. J Phys Conf Ser 2587 012032 DOI:10.1088/1742-6596/2587/1/012032. [6] Cavazzuti, M., Baldini, A., Bertocchi, E., Costi, D., Torricelli, E., Moruzzi, P. High performance automotive chassis design: A topology optimization based approach. Struct Multidiscip Optim 44 45-56 (2011) DOI:10.1007/s00158-010-0578-7. [7] Li, S., Zhou, D., Pan, A. Integrated lightweight optimization design of wall thickness, material, and performance of automobile body side structure. Struct Multidiscip Optim 67 95 (2024) DOI:10.1007/s00158-024-03810-1. [8] Yu, L., Gu, X., Qian, L., Jiang, P., Wang, W., Yu, M. Application of tailor rolled blanks in optimum design of pure electric vehicle crashworthiness and lightweight. Thin-Walled Struct 161 107410 (2021) DOI:10.1016/j.tws.2020.107410. [9] Del Pero, F., Berzi, L., Antonacci, A., Delogu, M. Automotive lightweight design: Simulation modeling of mass-related consumption for electric vehicles. Machines (Basel) 851 (2020) DOI:10.3390/machines8030051. [10] Xu, X., Chen, X., Liu, Z., Xu, Y., Zhang, Y. Reliability-based design for lightweight vehicle structures with uncertain manufacturing accuracy. Appl Math Model 95 22-37 (2021) DOI:10.1016/j.apm.2021.01.047. [11] Shangguan, Y., Wang, W., He, A., Ma, W., Tian, H. Lightweight design for the aluminum alloy-carbon fiber hybrid structure of the EMU car body. J Mech Sci Technol 37 6441-6452 (2023) DOI:10.1007/s12206-023-1116-z. [12] Fu, L., Zhang, S., Qiu, P., Xiao, H., Deng, B., Lu, X. Optimization of clinching joint process with preforming between ultra-high-strength steel and aluminum alloy sheets. Metals (Basel) 14 767 (2024) DOI:10.3390/met14070767. [13] Dogan, O., Karpat, F., Yuce, C., Kaya, N., Yavuz, N., Sen, H. A novel design procedure for tractor clutch fingers by using optimization and response surface methods. J Mech Sci Technol 30 2615-2625 (2016) DOI:10.1007/s12206-016­0522-x. [14] Zhang, Y., Shan, Y.C., Liu, X.D., He, T. Integrated shape-morphing and surrogate model-assisted multi-objective structural optimization of an automobile wheel made of lightweight materials to improve its impact resistance. Eng Optim 56, 1972-1998 (2023) DOI:10.1080/0305215x.2023.2291481. [15] Wang, Z.H., Zhang, T.Q., Zhang, Y., Jiang, Z.P., Wu, F.H. Shape optimization method for wheel rim of automobile wheels based on load path analysis. P I Mech EnC-J Mech 237 267-280 (2023) DOI:10.1177/09544062221119360. [16] Cheng, M.-Y., Prayogo, D., Wu, Y.-W., Lukito, M.M. A hybrid harmony search algorithm for discrete sizing optimization of truss structure. Autom Constr 69 21­33 (2016) DOI:10.1016/j.autcon.2016.05.023. [17] Goldberg, D.J.O. Genetic Algorithms in Search, Optimization and Machine Learning (1989). Addison-Wesley Professional, Boston. [18] Paz, J., Díaz, J., Romera, L., Costas, M. Crushing analysis and multi-objective crashworthiness optimization of gfrp honeycomb-filled energy absorption devices. Finite Elem Anal Des 91 30-39 (2014) DOI:10.1016/j.finel.2014.07.006. [19] Velea, M.N., Wennhage, P., Zenkert, D. Multi-objective optimisation of vehicle bodies made of frp sandwich structures. Comp Struct 111 75-84 (2014). DOI:10.1016/j.compstruct.2013.12.030. [20] Duan, L., Li, G., Cheng, A., Sun, G., Song, K. Multi-objective system reliability-based optimization method for design of a fully parametric concept car body. Eng Optim 49 1247-1263 (2017) DOI:10.1080/0305215x.2016.1241780. [21] Jiang, R., Jin, Z., Liu, D., Wang, D. (2021). Multi-objective lightweight optimization of parameterized suspension components based on NSGA-II algorithm coupling with surrogate model. Machines (Basel) 9 107 DOI:10.3390/machines9060107. [22] Xie, C., Wang, D., Wang, S., Kong, D., Du, C. A performance-driven lightweight optimization method for material-structure coupling design of cab. Struct Multidiscip Optim 67 12 (2024) DOI:10.1007/s00158-023-03730-6. [23] Qian, J., Yi, J., Cheng, Y., Liu, J., Zhou, Q. A sequential constraints updating approach for kriging surrogate model-assisted engineering optimization design problem. Eng Comput 36 993-1009 (2020) DOI:10.1007/s00366-019-00745-w. [24] Kaintura, A., Spina, D., Couckuyt, I., Knockaert, L., Bogaerts, W., Dhaene, T. A Kriging and stochastic collocation ensemble for uncertainty quantification in engineering applications. Eng Comput 33 935-949 (2017) DOI:10.1007/s00366­017-0507-0. [25] Deb, K., Pratap, A., Agarwal, S., Meyarivan, T. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans Evol Comput 6 182-197 (2002) DOI:10.1109/4235.996017. [26] Zienkiewicz, O.C., Zhu, Z.J. The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity. Int J Numerical Meth Eng 33 1365-1382 (1992) DOI:10.1002/nme.1620330703. [27] Kim, S.Y. Constrained sensitivity filtering technique for topology optimization with lower computational expense. Int J Numer Meth Eng 124 4075-4096 (2023). DOI:10.1002/nme.7301. [28] Kudela, J., Matousek, R. Recent advances and applications of surrogate models for finite element method computations: A review. Soft Comput 26 13709-13733 (2022) DOI:10.1007/s00500-022-07362-8. [29] Juliani, M.A., Gomes, W.J.S. An efficient Kriging-based framework for computationally demanding constrained structural optimization problems. Struct Multidiscip Optim 65 4 (2022) DOI:10.1007/s00158-021-03095-8. [30] Krige, D.G. A statistical approach to some basic mine valuation problems on the witwatersrand. J Chem Metall Min Soc S Afr 52 119-139 (1951). [31] Van, A.-L., Nguyen, T.-C., Bui, H.-T., Dang, X.-B., Nguyen, T.-T. Multi-response optimization of gtaw process parameters in terms of energy efficiency and quality. Stroj vestn-J Mech E 70 259-269 (2024) DOI:10.5545/sv-jme.2023.890. [32] Huang, W., Huang, J., Yin, C. Optimal design and control of a two-speed planetary gear automatic transmission for electric vehicle. Appl Sci (Basel) 10 6612 (2020) DOI:10.3390/app10186612. Acknowledgements This work was supported by the National Natural Science Foundation of China under Grant (51875494). Received 2024-10-05, revised 2025-04-07, 2025-04-30, accepted 2025-05-13 as Original Scientific Paper. Data availability The data that support the findings of this study are available from the corresponding author upon reasonable request. Author contributions Fu-sheng Ding: methodology, formal analysis, validation, writing – original draft; Hong-ming Lyu: supervision, project administration, methodology, review & editing; Jun Chen: conceptualization, methodology, review & editing. Hao-ran Cao: methodology, formal analysis, data curation. Lan-xiang Zhang: formal analysis, data curation. All authors have contributed significantly to this work and thoroughly reviewed and approved the final version of the manuscript. AI-Assisted writing AI-based tools were employed solely to enhance the readability and grammatical accuracy of the manuscript. Veckriterijska optimizacija zasnove iztisne plošce za smetarska vozila z zadenjskim nakladanjem Povzetek V delu je predstavljena veckriterijska optimizacija zasnove iztisne plošce, kljucne komponente smetarskih vozil z zadenjskim nakladanjem, z namenom ohranjanja trdnosti ob hkratni optimizaciji mase. Zasnovan je bil parametricni model s koncnimi elementi iztisne plošce, pri cemer so bili cilji optimizacije usmerjeni v zmanjšanje mase, maksimiranje omejitev deformacij ter zmanjšanje maksimalne von Misesove napetosti. Z analizo obcutljivosti je bilo identificiranih sedem kljucnih spremenljivk. Sistematicna zasnova parametrov je bila izvedena s pomocjo Box-Behnkenovega nacrtovanja eksperimentov (BBD), ciljna funkcija pa je bila približana s Krigingovim nadomestnim modelom, katerega ucinkovitost je bila primerjana z metodo odzivne površine (RSM). Za dolocitev optimalne rešitve je bil uporabljen veckriterijski genetski algoritem NSGA-II, ki omogoca lahkotnejšo zasnovo ob izpolnjevanju vseh strukturnih zahtev. Rezultati kažejo, da je možno maso iztisne plošce smetarskega vozila z zadenjskim nakladanjem z optimizacijo zmanjšati za 6,06 %, pri cemer so izpolnjeni kriteriji trdnosti in deformacij. Ta izboljšava ne poveca le ucinkovitosti transporta odpadkov, ampak tudi zmanjša proizvodne stroške in izboljša izkorišcenost gradiv. Kljucne besede smetarsko vozilo, iztisna plošca, veckriterijska optimizacija, NSGA-II, Kriging Identification Method of Tire-Road Adhesion Coefficient Based on Tire Physical Model and Strain Signal for Pure Longitudinal Slip Jintao Zhang – Zhecheng Jing – Haichao Zhou – Yu Zhang – Guolin Wang Jiangsu University, School of Automotive and Traffic Engineering, China hczhou@ujs.edu.cn Abstract To precisely calculate the tire-road adhesion coefficient of rolling tires at various slip rates, and enhance the safety and stability of vehicle operation, an approach for estimating the tire-road adhesion coefficient based on strain sensors and brush models was proposed. First, a finite element model of 205/55R16 radial tire was established, and the effectiveness of the model was verified through static ground contact and radial stiffness experiments. Then, the circumferential strain signal of the inner liner centerline of the tire during braking was extracted, utilizing the average peak angle spacing of the first-order and second-order circumferential strain curves, and the contact area length was estimated using the arc length formula. Subsequently, the braking simulation of rolling tires confirmed the asymmetry of pressure distribution within the ground contact area, estimating the position of slip points within the contact area based on arbitrary pressure distribution function and brush model, while nonlinear regression was utilized to fit the estimation function of slip point under various slip rates. Finally, a functional relationship was developed between tire-road adhesion coefficient and slip rate, considering the friction characteristics between tire rubber and road surface, while the friction model is based on exponential decay. The results suggest that the methods described above enable estimation of the tire-road adhesion coefficient under different slip rates, providing valuable insights for intelligent tire applications in vehicle dynamics control. Keywords intelligent tire, tire-road adhesion coefficient estimation, slip point, slip rate, nonlinear regression Highlights . Estimated ground contact mark length using strain signals from the tire’s inner liner. . Analyzed pressure distribution in tire contact areaes under pure longitudinal slip. . Estimated slip points in contact areas and derived a function for various slip rates. . Created a function to estimate road adhesion under slip with successful validation. 1 INTRODUCTION The tire is the only part of the vehicle in contact with the road surface, transmitting the forces and torques required to move the vehicle. Therefore, it is important to accurately obtain the adhesion coefficient between the tire and the road surface during vehicle driving to improve the braking performance, driving smoothness, and handling stability of the vehicle [1,2]. At the same time, as one of the key components of the vehicle dynamics control system, the tire plays a "passive" role in vehicle dynamics control. The adhesion coefficient between the tire and the road surface cannot be directly obtained. With the continuous development of vehicle intelligence, the concept of the intelligent tire was proposed by installing sensors inside the tire to make the tire an ‘active’ component in vehicle dynamics control and to obtain contact information between the tire and the road surface directly [3]. At present, most researchers' methods of estimating tire-road adhesion coefficient through intelligent tires are mainly divided into experiment-based methods and model-based methods. In estimating the tire-road adhesion coefficient through experiment-based methods, most researchers try to measure the parameters related to friction between tire and road, such as tire noise, longitudinal and lateral tire deformation, and find the correlation between the sensor signal characteristics and the tire-road adhesion coefficient. Alonso et al. [4] utilized the different noise levels of cars driving on roads with different adhesion conditions to build a system architecture and signal processing algorithm for measuring noise, which was used to measure road adhesion conditions. Some other researchers estimate road conditions by using microphones to record distinct noise data vehicles generate as they traverse different roads [5,6]. Kuno et al. [7] studied a real-time monitoring method for asphalt pavement conditions based on charge-coupled device (CCD) cameras, using average absolute deviation and reference brightness signals to detect the glossiness of the pavement, determine the pavement condition, and improve the accuracy of vehicle avoidance control systems. Xu et al. [8] designed an intelligent tire system by installing a three-axis acceleration sensor in the center of the tire lining layer to extract the three-direction acceleration signal of the tire and analyze the characteristic parameters of the circumferential acceleration. Using the support vector machine algorithm, this approach can classify and recognize different road surfaces. Estimation of tire-road adhesion coefficients through experiment-based method is relatively demanding. Since the sensor needs to be installed inside the tire, the signal is easily interfered by external noise during vehicle driving and faces problems such as sensor power supply and high price [9]. At the same time, a large amount of signal data needs to be tested for signal characterization and feature extraction as training samples for artificial intelligence algorithms such as neural networks. Although the accuracy of the estimation is high, once it deviates from the training conditions, the accuracy of the estimation is greatly reduced. Usually, the estimated tire-road adhesion coefficient is a constant value and does not describe the process of tire-road adhesion coefficient variation with tire slip rate [10]. In model-based research, most researchers attempted to use simplified mathematical models to determine the adhesion coefficient Fig. 1. Flowchart for identifying tire-road adhesion coefficient between tire and road. There are three main types of mathematical methods used for this purpose: the vehicle dynamics model, the glide-based method, and the tire model. Matsuzaki et al [11] installed an acceleration sensor on the tire lining to extract the acceleration signal. They obtained deformation information and tire force in the tire connection area by integrating the signal. Finally, the obtained data parameters are plugged into the adhesion coefficient calculation formula, which is derived from the tire brush theory model, to estimate the tire-road adhesion coefficient. Rajamani et al. [12] delve into developing an independent friction coefficient estimation algorithm for each vehicle wheel. Three distinct parameters are formulated to estimate slip ratio and longitudinal tire force. Subsequently, a recursive least squares parametric formula is employed to estimate the friction coefficient. Nishihara et al. [13] defined the concept of grip margin using the brush model and friction circle theory, re-derived and analyzed the brush model theory, and completed the estimation of the tire-road adhesion coefficient. However, the pressure distribution of the brush model used was a quadratic parabolic distribution. The model-based method avoids the need for extensive experiments or simulations for data collection and is more generalizable than the training models of intelligent algorithms such as neural networks. However, in the model-based method, the estimation of the tire-road adhesion coefficient often requires an accurate mathematical model to describe the mechanical characteristics of the tire-road contact, such as the contact pressure distribution. At the same time, in the process of constructing the algorithm, the issue of necessary grounding parameter inputs is sometimes overlooked because of the limitations of the on-board sensors, which do not allow real-time access to certain feature parameters. In summary, the goal is to improve the vehicle's drive braking performance under purely longitudinal slip conditions and promote the development of intelligent chassis technology. This paper designs a method for estimating the tire-road adhesion coefficient at different slip rates by combining the sensor signals with a mathematical model of the tire, used to calculate the optimum slip rate and peak adhesion coefficient of a tire, as shown in Fig. 1. First, the finite element model of a passenger car tire 205/55R16 was established, and it is validated using the contact sizes and radial stiffness in S1. Second, based on the circumferential strain signal characteristics of the centerline of the inner liner of the rolling tire, the estimation methods for contact angle and contact area length were designed in S2. In S3, synthetic application of a tire arbitrary pressure distribution function, pure longitudinal slip model, and a computing method for the tire-road adhesion coefficient are proposed. To complete the computing method, the S3 includes three substeps: an arbitrary pressure distribution function based on the UniTire model is introduced to reflect the rolling tire contact pressure in the substep of S3.1. Then, with the help of the theoretical analysis of the pure longitudinal slip brush model, the calculation methods for tire longitudinal force and slip point position are pit forward in the substep S3.2. Finally, in S3.3, according to the friction theorem, and taken tire physical parameters and vehicle speed into consideration, the estimation function of the tire-road adhesion coefficient under different slip rates is built. Estimations of the tire-road adhesion coefficients under different slip rates of rolling tire are also achieved. 2 METHODS & MATERIALS 2.1 Tire Model Simplifications and assumptions were employed in the development of the finite element model of the passenger car tire 205/55R16. These are deemed less influential in terms of simulation precision and are outlined as follows: 1. To facilitate computations, define a two-dimensional plane as a rigid surface body for simulating the road surface. 2. Ignore the influence of temperature during tire rolling. 3. Ignore the effects of road surface roughness. A tire is not made of a single rubber material but various materialswith relatively complex structure. Therefore, a finite element model is established, which is divided into structural and tread pattern models [14]. First, the 205/55R16 tire was cut in the laboratory to obtain the cross-sectional structure of the tire. A two dimensional (2D) axisymmetric structure model was established based on the actual tire structure profile and material distribution, in which the rubber element types are CGAX3H and CGAX4H, and the reinforcing materials, such as the carcass and belts, are established by the SFMGA1 plane element. The rubber solid element is embedded with rebar elements to simulate the features of the tire cord rubber composite material. The geometric parameters of the reference tread pattern blocks are used to import a 2D cross-section tire model into ABAQUS software. The 2D cross-section tire model is rotated to create a 3D solid mesh model as shown in Fig. 2. Fig. 2. Finite element model for the passenger car tire 205/55R16; a) 2D cross-section tire model, b) tread pattern block, c) one tire patch, and d) 3D tire model The single-pitch 3D tread pattern solid model is established using the CATIA software and imported into the HYPERMESH software to generate the 3D mesh. Considering the incompressibility of rubber materials, the corresponding grid type is the hybrid element C3D8RH [15]. The 3D tire model, including the complex tread pattern, is shown in Fig. 2d with 123842 elements and 155703 nodes. We specified the contact properties between the tire tread and road as hard contact and model the variation of the friction coefficient between the tire tread and road with slip speed using the Coulomb friction model. The tire's rubber material exhibits typical viscoelastic, incompressible, and doubly nonlinear characteristics of both material and geometry. The accuracy of the material model significantly impacts the analysis's convergence during tire dynamic performance. Building upon our team's prior research [16,17], the Yeoh model represents the rubber material and is defined as: W = c1(I1–3)+c2(I1–3)2+c3(I1–3)3, (1) where W is strain energy density, in which c1, c2, c3 are material constants and I1 is the first strain invariant. 222 I ˜.°.°. ˜ tr()C , (2) 11 23 where .i (i =1,2,3) are stretch ratios and the square root of the right Cauchy-Green strain tensor (C). The selection of intelligent tire sensors mainly includes three-axis accelerometers, strain sensors, optical sensors, and others. Strain sensors offer advantages such as high precision, wide measurement range, fast response, durability, reliability, and ease of installation and integration. Compared to accelerometers, the signals from strain sensors are less susceptible to interference from tire noise. Based on the previous research conducted by the research team on the force-sensitive response regions of tires, it was found that under braking conditions, the deformation caused by the contact between the tire and the ground primarily occurs at the tread. Therefore, circumferential strain signals were extracted and analyzed by comparing the strain at the tread and the inner liner centerline, as shown in Fig. 3. The tread pattern significantly influences the circumferential strain signal at the tread centerline, resulting in large fluctuations in the signal, which are particularly evident at the signal's troughs. These fluctuations hinder the effective extraction of sensitive signal features. Fig. 3. Strain signal extraction locations Compared to the circumferential strain signal at the tread centerline, the strain signal at the inner liner centerline is less affected by the tread pattern as seen in Fig. 4. The signal curve is smoother and more stable, with its trend opposing the tread centerline strain signal. This discrepancy is primarily related to the tire's structure, material properties, and load distribution. When the tire is subjected to external loads and in contact with the ground, the tread experiences compression and elongation due to frictional forces, resulting in significant localized deformation. The tire structure is multi-layered and composite, and due to the differing elastic moduli and response characteristics between the tread and the inner liner, the inner liner undergoes deformation that is opposite to that of the tread. Therefore, by extracting the circumferential strain signal of the point in the inner liner centerline under the rolling tire braking condition, the signal is selected to represent the real-time strain sensor signal of the rolling tire. The specific location of this point is shown in Fig. 2d. Fig. 4. Comparison of circumferential strain signals between tread and inner liner 2.2 Experimental Verification To verify the accuracy of the tire finite element model established, the results were verified using the tire static ground contact distribution characteristics test and tire radial stiffness test. A static ground contact test was performed using a tire drum testing machine; during the test, the tire was loaded to a rated load of 3800 N, and the inflation pressure was set to 210 kPa. The ground contact area of the tire under static loading is obtained through red ink printing, and the ground contact geometric parameters are extracted, as shown in Figs. 5 and 6. Fig. 5. Tire stiffness test benchP Table 2. Comparison of test and simulation results of tire sinkage and static loading radius Fig. 6. Comparison of the sizes of tire contact ar ea Table 1. The tire ground contact characteristics Title Experiment values [mm] Simulate values [mm] Error [%] Contact width 147 148 0.67 Contact length 161 163.8 1.74 Table 1 shows that the error between the tire contact geometric properties obtained by finite element simulation and the test result is less than 2 %, indicating that the tire finite element model established in this paper can accurately describe the tire ground contact characteristics. To further verify the accuracy of the tire finite element model, a comprehensive strength testing machine was used to conduct tire radial stiffness tests. The test and simulation results of tire sinkage and static loading radius under different test conditions are shown in Table 2, and the test and simulation results of radial stiffness are shown in Fig. 7. It can be seen from the Table 2 and Fig. 7 that the tire sinkage, static loading radius, and radial stiffness of the test and simulation are consistent. 2.3 Estimate Contact Patch Length According to the simulation results of the circumferential strain signals of the liner centerline and tread centerline of rolling tire [18], this paper introduces a virtual strain sensor that utilizes strain information from a node on the longitudinal centerline of the tire inner liner, as illustrated in Fig. 2d. When the tire is subject to the 8 % slip rate, a 3800 N load, and 0.21 MPa inflation pressure, the circumferential strain signal of the tire inner liner centerline during braking is obtained through the virtual strain sensor. The ground pressure signal at the centerline of the tire's inner liner is obtained and analyzed. The results are shown in Fig. 8. .˜QP ° . 24 355. ˆ ˆ˜CD ° . 27 227. . , (3) ˜EF ° 21 176. . ˆ ˜° 33 277 ˆ . . . GH where .QP is actual rolling tire contact angle, .CD is the peak angle spacing of the second-order, .EF is the peak angle spacing of the first-order, .GH is the peak spacing of the circumferential strain. Experiment values Simulation values Test conditions 210 kPa 250 kPa 210 kPa 250 kPa / Sinkage [mm] Static radius [mm] Sinkage [mm] Static radius [mm] Sinkage [mm] Static radius [mm] Sinkage [mm] Static radius [mm] 123 kg 6.01 308.14 5.49 309.01 6.48 309.47 5.98 309.97 246 kg 12.28 301.87 11.13 303.37 12.17 303.78 11.62 304.33 369 kg 18.06 296.09 16.29 298.21 18.28 297.67 16.99 298.96 .QP is the actual rolling tire contact angle. It can be obtained that the difference between the rolling tire contact angle and the peak angle spacing of the first-order and second-order derivative curves of the circumferential strain is about 3°, and the difference with the peak angle spacing of the circumferential strain signal is about 10°. Therefore, the rolling tire contact angle can be estimated by averaging the circumferential strain signal's first- and second-order peak pitch angles, as shown in Eq. (4). ˜°.˜ .˜ . /,2 (4) CD EF where . is estimated contact angle. Fig. 8. Comparison between peak angle difference of circumferential strain of tire inner liner and contact angle To investigate the influence of rolling tires on the tire contact angle during braking, we conducted simulation analysis of pure longitudinal slip conditions under free rolling and slip rates of 2 %, 4 %, 6 %, 8 %, and 10 %, with load of 3800 N and inflation pressure of 0.21 MPa. Figure 9 shows the simulation results. During free-rolling, the front and rear contact angles are equal. When the tire begins to slip, the position of the midpoint of the tread in the finite element simulation starts to move toward the front. As the slip ratio increases, the front contact angle increases while the rear contact angle decreases, increasing the asymmetry of the contact angles. When the tire is in a high slip ratio, the trend of increasing contact angle asymmetry slows down, consistent with the results in [19]. Therefore, when estimating the tire's contact area length during braking based on the contact angles, the assumption of equal front and rear contact angles should not be used. Instead of using trigonometric functions, the arc length formula does not require distinguishing between the magnitudes of the front and rear contact angles. This makes it more suitable for estimating the contact area length of a rolling tire during the braking process. The formula for estimating the contact area length is given by Eq. (5) L ˜ (/ 360)(° 2. R), . (5) where L is estimated contact area length. Figure 10 shows a rolling tire’s estimated contact area length under different loads, compared with the contact area length obtained from finite element simulation. It can be observed that since the arc length is greater than the straight-line length, the estimated contact area length is larger than the simulation result, with a relative error exceeding 6 %. Under low and high load conditions, the error between the estimated and simulated values is minor, and it becomes larger under moderate load. To reduce the error in the estimated contact area length, a correction factor ß accounts for the length difference between the arc and the straight line, as shown in Eq. (6). L ˜ (/. 360)(° 2. R) .ˆ. (6) To determine the value of the correction factor ß, values of ß were chosen as 1 mm, 2 mm, 3 mm, and 4 mm to calculate the relative error in the contact area length for each correction factor, as shown in Fig. 11. It can be observed that as the value of ß increases, the estimation accuracy of the contact area length improves. After considering the overall performance for different values of ß, a correction factor ß of 3 mm was selected. Simulation results were further used to obtain .1 of the grounding angle and L1 of the contact area length under different loads. Estimated values for the contact angle .2 and the contact area length L2 are determined using Eqs. (4) and (6). Comparison shown in Fig. 12 reveals that the disparity between simulated value .1 and estimated value .2 of the contact angle is less than 1°, and the maximum error between simulated value L1 and estimated value L2 of the contact area length does not exceed 3 mm, verifying the accuracy of this estimation method. Fig. 12. Comparison of contact parameter estimation results under different loads 3 RESULTS AND DISCUSSION 3.1 Contact Pressure Distribution of Rolling Tires Accurately describing the pressure distribution within the tire's contact area is crucial for ensuring the accuracy of estimates concerning longitudinal forces, slip point, and tire-road adhesion coefficients. Typically, pressure distribution within the tire's contact area is assumed to follow a symmetrical parabolic distribution. However, this symmetrical parabolic distribution only applies when the tire is in a state of free rolling within small loads. Under varying working conditions, a symmetrical parabolic distribution cannot accurately represent actual tires. Hence, according to the tire contact pressure distribution function presented in UniTire model [20], this article utilizes arbitrary pressure distribution within the UniTire model to uniformly express the pressure distribution of rolling tires during braking. The pressure distribution of the tire is defined as follows: Fz 2n qu()˜ A° (1. u )° (1. Bu), (7) z 2a where Fz is the vertical load, a is the half-length of the contact area, and A, B are undetermined coefficients. The vertical load offset of the tire is set to .. According to the mechanical characteristics of the tire, the relationship between contact pressure distribution and tire vertical load Fz needs to be satisfied: . 1 aq udu.F °z() z .˜1 . 12 ˆ aqz()uudu.Fz. , (8) ° .˜12n ()..˜ .˜ ) uA(1 u )(1 Bu .  while the following boundary conditions should be met: ˜()°.°˜10 1 ()  ˜()u .0,u.[.1,]1  ˜()u °0,uˆ[.1,]1  1. (9) .˜()udu°2 .1 1  ˜().. °. uudu 2 ..1 a  According to the boundary conditions the expression of A and B can be defined as: . 2n °1 A ˜ . 2n ˆ . (10) 32 3) .(n+ .B ˜. .. . 2n °1 a The influence of n and ./ a on the contact pressure distribution function is shown in Fig. 13. The parameter n significantly influences the uniformity of pressure distribution, while parameter ./ a mainly affects the peak distribution on either side of the tire. Table 3. Simulation schemes Number Load [N] Inflation pressure [MPa] Friction coefficient 1 2800 0.21 0.5 2 3800 0.21 0.5 3 4800 0.21 0.5 4 3800 0.26 0.5 5 3800 0.31 0.5 6 3800 0.21 0.1 7 3800 0.21 0.9 a) To verify the contact pressure distribution characteristics of a rolling tire during braking, it is necessary to determine the values of the parameters for the arbitrary pressure distribution function and assess whether this function can accurately describe the contact pressure distribution between the tire and the ground during braking. We designed seven different finite element simulation schemes using the control variable method to determine the effect of load, inflation pressure, and longitudinal force on pressure distribution under free rolling and braking conditions. The details of each simulation scheme are shown in Table 3. Figure 14 illustrates the pressure distribution within the contact area of rolling tires under varying loads. The tire inflation pressure is 0.21 MPa, the speed is 70 km/h, the friction coefficient is 0.5, and the slip rate of rolling tires during braking is 8 %. For the convenience of research, the pressure distribution simulation result is obtained by extracting the longitudinal centerline of the tire contact area. As shown in Fig. 14a, as the weight on a tire increases or decreases, and the length of its contact area with the ground also changes. When the tire bears more weight, the contact area becomes longer due to radial deformation of the tire. This happens because the tire deforms radially when subjected to greater loads. Additionally, as the weight increases, the highest-pressure point in the center of the b) a) b) Fig. 14. Comparison of contact pressure distribution under the different loads: a) free rolling; b) braking tire's contact area decreases while the highest pressure on both sides increases. Similar simulation results are obtained under braking conditions, as shown in Fig. 14b. The comparison of contact pressure between the free rolling and the braking conditions shows that pressure distribution is approximately symmetric in free rolling but asymmetric in braking conditions. Figure 15 illustrates the pressure distribution within the contact area of rolling tires under varying inflation pressure; the tire load is 3800 N, the speed is 70 km/h, the friction coefficient is 0.5, and the slip rate of rolling tires during braking is 8 %. Figure 15a indicates that when the inflation pressure is 0.31 MPa, the contact area length of the rolling tires is less than that of 0.21 MPa. The tire's overall stiffness changes with the inflation pressure when the load remains unchanged. The higher the tire's inflation pressure, the less it deforms against the ground and the smaller the ground contact area becomes. Similar results are obtained under braking conditions. Figure 16 illustrates the pressure distribution within the contact area of rolling tires under varying friction coefficients; the tire inflation pressure is 0.21 MPa, the speed is 70 km/h, the load is 3800 N, and the slip rate of rolling tires during barking is 8 %. Figure 16a shows that under free rolling conditions, the peak value of tire contact pressure increases with the increase of friction coefficients. However, the increasing trend is relatively small, and the length of the tire contact area remains unchanged. In Figure 16b, it can be observed that under braking conditions with an 8 % slip ratio, the pressure distribution between the tire and the ground shifts toward the front. Furthermore, as the road surface friction coefficient increases, the higher friction coefficient intensifies the frictional force between the tire and the ground, making the changes in the front end of the contact pressure distribution more pronounced. Based on the finite element simulation results of contact pressure distribution under different operating conditions shown in Figs. 14­16, it is demonstrated that the contact pressure distribution during tire braking is not parabolic symmetrical. Instead, the pressure distribution shifts toward the front end of the contact area depending on changes in operating conditions. In Figure 13, the arbitrary pressure distribution function solved in MATLAB can be adjusted by changing the values of parameters, thereby modifying the shape of the pressure distribution. Therefore, to describe the forward shift in contact pressure distribution during braking, parameters are set to –0.04 and 2, respectively, to approximate the characteristics of the contact pressure distribution between the rolling tire and the ground under braking conditions. 3.2 Estimate the Slip Point Based on Tire Brush Model Under the condition of pure longitudinal slip, the tire only slips in the rolling direction, and the slip angle and lateral speed are zero, therefore, the linear speed of the tire is equal to the vehicle speed. During a vehicle's braking process, the tires undergo three phases: free-rolling, sliding while rolling, and pure slip, which is a dynamic and progressive process. Under braking conditions, a tire's slip rate a) a) is calculated from the linear and rolling tire angular velocities at a known effective rolling radius, which is defined as shown in Eq. (11). v °.r Sx ˜ , (11) where v and . represent the vehicle speed and wheel angular speed respectively; r is the free-rolling radius of the wheel. Accurate calculation of a tire's slip rate during vehicle braking requires accurate determination of the tire's linear and angular velocities as well as the size of the effective rolling radius. Wheel speed sensors on current vehicles can accurately measure the rolling tire's angular velocity while the vehicle is in motion, and the tire's linear velocity can be obtained from a body viewer or global positioning system (GPS). Several researchers have implemented the observation of longitudinal and transverse vehicle velocities using whole-vehicle dynamics models. However, the observed velocities represent the traveling speed at the vehicle's center. In order to study the relationship between the linear velocities of the four tires and the longitudinal and transverse velocities at the center of mass of the vehicle when the vehicle is moving, the relationship between the linear velocities of the four tires and the longitudinal and transverse velocities at the center of mass of the vehicle can be obtained by analyzing the whole vehicle's seven-degree-of-freedom dynamics model, which can be expressed as: . v. fl . .ˆvfr . ˜ (vx ˜ (vx ° . B r) cos. 2 B r) cos. 2 . v( y . (vy . lr)sin .f . lr)sin . ,f (12) . .. vrl ˜ vx ° B rv, rl 2 ˜ vx . B r 2 where vfl, vfr, vrl, and vrl denote the tire speeds of the left front wheel, right front wheel, left rear wheel, and right rear wheel, respectively; d is the tire angle of rotation; B is the wheelbase between the tires of the vehicle; lf and lr denote the axle distances between the front and rear axles of the vehicle. A uniform definition of the slip rate of a tire based on the update of the imprint coordinates is used in the Unitire tire model. Based on this idea, this paper expresses the bristle slip within the contact area by determining the slip rate at a point in the center of the contact area to simplify the description of the bristle slip at different parts. The brush model is a basic physical representation of a tire. It simplifies the tire as a rigid ring connected to a set of bristles. This model assumes that the bristles undergo elastic deformation when the tire touches the ground and carries the vertical, longitudinal, and lateral loads. During pure longitudinal slip of the rolling tire, the contact area is partitioned into a slip zone and an adhesion zone [21]. The location of the slip point is critical for delineating the slip zone and the adhesion zone. Moreover, the slip points vary with different slip rates, influencing the longitudinal force's magnitude within the contact area, as shown in Fig. 18. When there is a relative slip between the tire and the road contact area, the friction between the tire and the ground is dynamic. In such cases, the longitudinal force can be calculated by taking the integral of the vertical load of the bristles and the friction coefficient. This is shown in Eq. (13): 1 F ˜ a uq()udu, (13) sl ° sz u c where a is half the length of the tire contact area; uc is the slip point; us is the friction coefficient; qz(u) is a function for arbitrary pressure distribution in the tire contact area. When there is no relative slip between the tire and the road surface, the longitudinal force at each point of the tread in the adhesion zone is the product of the longitudinal stiffness and the longitudinal deformation of the tread, expressed as: kS qu()˜°kx˜tx xau. (14) x tx 1.Sx where ktx is the longitudinal stiffness of the element bristles. Our research group has identified the stiffness of 205/55R16 radial tires [22], and the longitudinal stiffness of the tire selected in this article is 3680000 N/m2. The normalized coordinate of the tire contact area is defined wit symbol u; through integration, the longitudinal force expression of the grounding imprinting adhesion area can be obtained as: uc kS Fad ˜ a tx x audu. (15) . . 1 ° Sx 1 The slip point is the demarcation point between the adhesion and slip zones. At the slip point, the longitudinal forces in the adhesion and slip zones are equal. Therefore, the joint Eqs. (13) and (15) will lead to Eq. (16). kS tx xau ° uq (). (16) u c szc 1˜ Sx Upon inputting known parameters into Eq. (16) and utilizing MATLAB software to solve the calculation, the expression is further articulated as: SV-JME . VOL 71 . NO 5-6 . Y 2025 . 187 54 kS u ˜ AFu (Bz ° z ° Bz °1) ° 2 tx xz , (17) c zsii i i 1 . Sx where zi (i = 1, 2, 3, 4, 5) results in the equation having a total of 5 roots. By analyzing Eq. (17) for the slip point uc, it can be concluded that the slip rate, load, and longitudinal stiffness directly affect the calculation results of the slip point uc. Meanwhile, since Eq. (17) for the slip point uc is a complex equation with 5 roots, the process of solving is too cumbersome, and sometimes judgment failure may occur under the condition of small slip rates. Therefore, it is meaningful to simplify the solution process. As shown in Fig. 19, the tire-road adhesion coefficient estimation function designed in this paper for different slip rates is mainly for the tire’s braking condition, and the longitudinal stiffness of the tire tread does not change during the vehicle's driving process. Based on the slip point solution in Eq. (17), the conjecture assumes the formula defined in Eq. (18) for calculating the slip point at different slip rates. The process of assumption of the fitting formula Eq. (18) is as follows: 1. This paper mainly estimates the tire-road adhesion coefficient under different slip rates, so the slip rate is set as a main variable. 2. According to the analysis of the slip point solution formula in Eq. (17), it is found that tire load and longitudinal stiffness of tread also affect the calculation accuracy of slip point. 3. Tire load and longitudinal stiffness of tread are introduced into the assumption of slip point fitting formula into Eq. (18). 4. Finally, through adjusting the order of the fitting equation function defined in Eq. (18), the slip point fitting function was established and solved. u˜a°.S 3 aS 2 .akS . (18) .°. .aF c 1 x 2 x 3 txx 4 z Figure 20 demonstrates the slip point fitting results for different slip rates at a load of 4000 N with constant longitudinal tire tread stiffness. Equation (18) can be fitted fairly good, which proves that the assumptions made in Eq. (18) in this paper are correct. This article exploits the calculation method to obtain the split point uc under six load conditions of 2000 N, 2500 N, 3000 N, 3500 N, 4000 N, and 4500 N for data fitting, the slip rate ranges from 0.02 to 0.6 with intervals of 0.02, the friction coefficient is set at 0.5, and the tire’s rolling speed is maintained at 19.44 km/h. Table 4 shows the values of a1, a2, a3 and a4 under six groups of load conditions. From Table 4, it is evident that with the increase in load, the values of fitting coefficients a1, a2 and a3 exhibit a gradual decrease, and the rate of decrease diminishes progressively. The fitting coefficient a4 undergoes a fluctuating pattern with increasing load, showing a nonlinear relationship. To make the formula of split point uc respond to the changing trend of tire load, the relationship between a1, a2, a3 and a4 coefficients and load is analyzed respectively, nonlinear fitting is carried out in the MATLAB fitting toolbox. The fitting relationship between a1, a2, a3 and a4 and the load is shown in the following equation: . a ˜ pF3 ° pF2 ° pF ° p 11 z 2 z 3 z 4 . 32 a ˜ pF ° pF ° pF ° p . 21 z 2 z 3 z 4 . 5 43 2 . (19) a ˜ pF ° pF ° pF ° pF ° pF ° p . 31 z 2 z 3 z 4 z 5 z 6 . 5 43 2 a ˜ pF ° pF ° pF ° pF ° pF ° p ˆ 41 z 2 z 3 z 4 z 5 z 6 Table 4. Fitted values of a1, a2, a3, a4 under six sets of loads Load [N] a1 a2 a3 a4 RMSE 2000 –177.6 110 –0.000006653 0.0005979 0.02931 2500 84.13 66.7 –0.000005148 0.0005864 0.02986 3000 –55.04 49.68 –0.000004399 0.0005921 0.02685 3500 –30.6 33.89 –0.00000363 0.0005819 0.03029 4000 –19.68 25.36 0.000003129 0.000289 0.0324 4500 –18.63 24.41 0.00000308 0.0006173 0.02453 Table 5 shows the fitting data of a1, a2, a3, a4. To further validate the accuracy of the fitted slip point estimation formula, we substituted the estimated slip point into the longitudinal force formula and compared it with the longitudinal force obtained through finite element simulation. Three simulation scenarios were established to simulate the transition of the rolling tire from free rolling to lock up under load conditions of 2800 N, 3800 N, and 4800 N. The inflation pressure was set at 0.21 MPa, and the coefficient of friction was set to 0.5. Figure 21 presents the tire’s estimated and simulated longitudinal forces under three different load conditions. It can be observed that, with the increase in slip rate, the simulated longitudinal force of the tire increases almost linearly, reaching a maximum value around a slip rate of 10 %, after which it remains constant. The load influences the estimated longitudinal force. When the load is 2800 N, the trends of the simulated and estimated longitudinal forces are generally consistent. As the load increases, the estimated longitudinal force rises more slowly in the low-slip region, and the peak shifts to a higher slip rate. In the high-slip region, the estimated longitudinal force aligns closely with the simulated values, indicating that the Table 5. Fitted values of p1, p2, p3, p4, p5, p6 p1 p2 p3 p4 p5 p6 R2 1.822×10–8 2.158×10–4 0.8604 –1179 / / 0.9925 –5.586×10–9 7.201×10–5 –0.3165 498.5 / / 0.9947 –6.392×10–21 9.935×10–17 –6.044×10–13 1.8 ×10–9 –2.622 ×10–6 0.001489 0.9999 3.744×10–19 –5.771×10–15 3.489×10–11 –1.035 ×10–7 1.505 ×10–4 –0.08527 0.9998 tire’s longitudinal force estimated at the slip point under Coulomb on surface roughness [23]. In [23], the authors compared the friction friction is accurate, with virtually no error in the high-slip region. forces under Coulomb friction and exponential decay friction models of the tire contact area, which is expressed as: uc kS 1 a tx x audu . auszqudu () .. F . F 1. S ad sl .1 x uc ˜° ° . (20) Fz Fz Given the tire longitudinal stiffness and load one needs to decide upon the appropriate fricrion model. The friction mechanism between tire rubber and the road surface is complex, involving factors such as the viscoelasticity of the rubber material, road surface roughness, and deformation in the tire contact area. It can be divided into adhesive friction and hysteresis friction [23]. The friction between tire rubber and the road surface is nonlinear, meaning the friction coefficient is not constant. When the slip ratio is low, adhesive friction plays the primary role during braking. As the slip ratio increases, hysteresis friction gradually becomes more prominent. When the slip ratio approaches full slip, the friction force tends to saturate and slightly decrease. Therefore, this paper adopts an exponential decay friction model to describe the nonlinear friction characteristics between tire rubber and the ground, aiming to improve the accuracy of friction force prediction under dynamic conditions, as shown in Eq. (21). . u ˜u °(u .u )e..s , (21) sk hk where uk is the dynamic friction coefficient at the highest slip speed, uh is the static friction coefficient at zero slip velocity, a is a self-defined attenuation coefficient according to road roughness, and s is the slip speed. The static friction coefficient is closely related to the roughness of the road surface, while the decay coefficient is highly dependent To validate the Eq. (20), which incorporates an exponential decay friction model and can accurately estimate the longitudinal force and friction characteristics between tire rubber and road surface, we obtained the optimal slip ratio and peak adhesion coefficient during braking, and provided accurate inputs for active safety control of vehicles. This study used finite element simulation software (Abaqus) to perform simulation analysis of the tire transitioning from free rolling to full slip under three load conditions: 2800 N, 3800 N, and 4800 N. A comparison between the finite element simulation results and the estimates from Eq. (20) is shown in Fig. 23, which illustrates the comparison between the simulated longitudinal force and the estimated longitudinal force obtained using the exponential decay friction model. Using the exponential friction model, the peak longitudinal force changes with the load increase; when the load is 2800 N, 3800 N, and 4800 N, the peak longitudinal force appears near the slip rate of 14 %, 18 %, and 22 %. After that, it began to decline, and the simulation results were consistent with the estimated results. Based on the simulation results, the feasibility of Eq. (20) in estimating the tire-road adhesion coefficient at different slip rates was verified. The vehicle speed sensor and wheel speed sensor measure the tire’s linear and angular velocities, and the tire’s slip rate is solved according to Eq. (7). Finally, the estimation of the tire-road adhesion coefficient at different slip rates is completed based on the method of estimating the contact area length and the slip point solution method established in Section 2.3 and 3.2. Finally, the tire-road adhesion coefficient under different slip rates can be estimated by inputting vehicle speed, tire angular velocity, longitudinal stiffness, load, and slip point. Figure 24 shows the estimated values of tire-road adhesion coefficient under different slip rates when the decay coefficient of the exponential decay friction model is 0.05, 0.1, 0.2, 0.3, 0.4, 0.5. It can be obtained that the estimated tire-road adhesion coefficient rises rapidly in an approximately linear relationship before reaching the peak adhesion coefficient when the slip rate is 15 % to 20 %, and then begins to decline slowly. This is due to the local increase of relative slip in the tire contact area before reaching the peak adhesion coefficient, resulting in slower tire-road adhesion coefficient with the increase of slip rate. Because the dynamic friction factor between friction pairs is smaller than the static friction factor, the estimated tire-road adhesion coefficient gradually decreases after reaching the peak adhesion coefficient, consistent with the resulting trend in the references [24,25]. As noted in the previous discussion and based on [26], the decay rate is related to the surface roughness of the road. It can be observed that a higher decay rate results in a lower peak adhesion coefficient between the tire and the road surface, as well as a lower optimal slip rate, which is consistent with [27,28]. Therefore, the tire-road adhesion coefficient estimation under different slip rates designed in this paper can effectively identify the optimal tire slip rates and the vehicle’s peak tire-road adhesion coefficients during the driving process. It provides accurate inputs for the whole vehicle’s drive anti-skid control and, simultaneously, can improve the control accuracy of the vehicle’s anti-lock breaking system (ABS) and other active safety control systems [29]. 4 CONCLUSION This study combines strain sensors with a brush model to develop a method for estimating the tire-road adhesion coefficient of rolling tires at various slip rates during braking. The following conclusions were drawn from this study: 1. During the braking process of rolling tires, the average peak angle spacing between the first-order and second-order curves of the tire's inner layer centerline circumferential strain can effectively represent the tire contact angle, and an estimation method for the contact angle of rolling tires is put forward by using the circumferential strain signal characteristics. Considering the asymmetry of the front and rear contact angles of rolling tires in a braking state, an arc length formula is employed to estimate the rolling tire contact area length; the arc length formula calculates the arc length based on the central angle and radius, with a defined factor aimed at minimizing the discrepancy between the arc length and the straight line. The maximum relative error of the estimated values using the arc length formula compared with the simulation results is less than 3 mm. This method offers an approach for utilizing innovative tires to estimate the contact area length. 2. Considering that the rolling tire contact pressure distribution is affected by conditions such as load, inflation pressure, friction coefficient, etc., an arbitrary pressure distribution function is selected to describe the pressure distribution of the rolling tire during braking. The slip point is the pivotal factor for distinguishing between the adhesion and slip zones. Given the influence of load, tire longitudinal stiffness, and slip rate on the rolling tire slip point, as well as the high-order function of the pressure distribution function, a nonlinear regression equation is used to establish the estimation function for the slip point. The comparison between simulation values and estimated values of longitudinal force confirms the accuracy of the estimation method using the nonlinear regression equation. 3. An estimation method based on the designed contact area length and slip point, employing the brush model, developed a method for estimating tire-road adhesion coefficient. This method establishes the functional relationship between the tire-road adhesion coefficient and the slip rate, determining the required input parameters: slip rate, load, tire longitudinal stiffness, slip point, and contact area length. Assuming a known load, the tire-adhesion coefficient estimation method presented in this paper enables the estimation of the road adhesion coefficient of rolling tires at various slip rates during braking. By integrating sensors and mathematical models, this method offers a technical solution for applying intelligent tires in vehicle control, contributing positively to the future development of intelligent tire technology. 4. Due to the difficulty of controlling the slip rate and measuring the slip point location, this study did not conduct actual vehicle experiments on various road surfaces. The method designed in this paper for estimating the tire-road adhesion coefficient under different slip rates mainly considers the pure longitudinal slip condition of the tire. It does not take into consideration the vehicle's lateral slip, composite conditions during driving, and estimation of parameters such as wheel steering, slip angle, wheel alignment angle, etc., which has some limitations on the prospect of the vehicle's application. In future research work, contact pressure distribution tests of rolling tires, identification and calibration tests of the parameters in the slip point fitting equations, and tests of tire-road adhesion coefficients at different slip rates will be conducted to optimize further the method of tire-road adhesion coefficient estimation established in this paper. References [1] Li, G., Xie, R., Wei, S., Zong C. Estimation of vehicle status and road adhesion coefficient based on dual-volume Kalman filter. Scie Sinic Technol 45 403-414 (2015) DOI:10.1360/N092014-00207. [2] Wang, Y., Liang, G., Wei Y. Intelligent tire pavement identification algorithm based on support vector machine. Autom Eng 42 1671-1678 (2021) DOI:10.19562/j. chinasae.qcgc.2020.12.009. [3] Khaleghian, S., Emami, A., Taheri, S. A technical survey on tire-road friction estimation. Friction 5 123-146 (2017) DOI:10.1007/s40544-017-0151-0. [4] Alonso, J., López, J. M., Pavón, I., Recuero, M., Asensio, C., Arcas, G., et al. On-board wet road surface identification using tyre/road noise and support vector machines. Appl Acoust 76 407-415 (2014) DOI:10.1016/j.apacoust.2013.09.011. [5] Kongrattanaprasert, W., Nomura, H., Kamakura, T., Ueda, K. Detection of road surface states from tire noise using neural network analysis. IEEJ Trans Indust Applic 130 920-925 (2010) DOI:10.1541/ieejias.130.920. [6] Erdogan, G., Alexander, L., Rajamani, R. Estimation of tire-road friction coefficient using a novel wireless piezoelectric tire sensor. IEEE Sens J 11 267-279 (2010) DOI:10.1109/JSEN.2010.2053198. [7] Kuno, T., Sugiura, H. Detection of road conditions with CCD cameras mounted on a vehicle. Syst Comput Jpn, 30 88-99 (2015) DOI:101002/(SICI)1520684X(19991 2)30:14<88::AID-SCJ9>3.0.CO;2-8. [8] Xu, N., Askari, H., Huang, Y., Zhou, J., Khajepour, A. Tire force estimation in intelligent tires using machine learning. IEEE T Intell Transp 23 3565-3574 (2020) DOI:10.1109/TITS.2020.3038155. [9] Xiong, Y., Yang, X. A review on in-tire sensor systems for tire-road interaction studies. Sensor Rev 38 231-238 (2018) DOI:10.1108/SR-07-2017-0132. [10] Rajamani, R., Phanomchoeng, G., Piyabongkarn, D., Lew, J.Y. Algorithms for real-time estimation of individual wheel tire-road friction coefficients. IEEE-ASME T Mech 17 1183-1195 (2011) DOI:10.1109/TMECH.2011.2159240. [11] Matsuzaki, R., Kamai, K., Seki, R. Intelligent tires for identifying coefficient of friction of tire/road contact surfaces using three-axis accelerometer. Smart Mater Struct 24 025010 (2014) DOI:10.1088/0964-1726/24/2/025010. [12] Rajamani, R., Piyabongkarn, N., Lew, J., Yi, K., Phanomchoeng, G. Tire-road friction-coefficient estimation. IEEE Contr Syst Mag 30 54-69 (2010) DOI:10.1109/MCS.2010.937006. [13] Nishihara, O., Masahiko, K. Estimation of road friction coefficient based on the brush model. T Jpn Soc Mech Eng 75 041006 (2011) DOI:10.1115/1.4003266. [14] Zhou, H., Li, H., Yang, J., Chen, Q., Wang, G., Han, T., et al. A strain-based method to estimate longitudinal force for intelligent tires by using a physics-based model. Stroj vestn-J Mech E 67 153-166 (2021) DOI:10.5545/sv-jme.2020.7068. [15] Li, B., Quan, Z., Bei, S., Zhang, L., Mao, H. An estimation algorithm for tire wear using intelligent tire concept. P I Mech Eng D J Aut 235 2712-2725 (2021) DOI:10.1177/0954407021999483. [16] Zhou, H., Li, H., Liang, C., Zhang, L., Wang, G. Relationship between tire ground characteristics and vibration noise. Stroj Vestn-J Mech E 67 11-27 (2021) DOI:10.5545/SV-JME.2020.6946. [17] Mei, F. Research on Finite Element Automatic Modeling Technology of Radial Tire Tread Pattern. PhD Thesis, Shandong University, Shandong (2020). [18] Gu, T., Li, B., Quan, Z., Bei, S., Yin, G., Guo, J., et al. The vertical force estimation algorithm based on smart tire technology. World Electr Vehic J 13 104 (2022) DOI:10.3390/wevj13060104. [19] Romano, L., Sakhnevych, A., Strano, S., Timpone, F.. A novel brush-model with flexible carcass for transient interactions. Meccanica, 54 1663-1679 (2019) DOI:10.1007/s11012-019-01040-0. [20] Guo, K., Lu, D., Chen, S.K., Lin, W.C., Lu, X. P. The UniTire model: a nonlinear and non-steady-state tyre model for vehicle dynamics simulation. Vehicle Syst Dyn 43 341-358 (2005) DOI:10.1080/00423110500140690. [21] Mavros, G., Rahnejat, H., King, P.D. Transient analysis of tyre friction generation using a brush model with interconnected viscoelastic bristles. P I Mech Eng K-J Mul 219 275-283 (2005) DOI:10.1243/146441905X9908. [22] Riehm, P., Unrau, H.J., Gauterin, F., Torbrügge, S., Wies, B. 3D brush model to predict longitudinal tyre characteristics. Vehicle Syst Dyn 57 17-43 (2019) DOI:10.1080/00423114.2018.1447135. [23] Wang, H., Al-Qadi, I.L., Stanciulescu, I. Effect of surface friction on tire-pavement contact stresses during vehicle maneuvering. J Eng Mech 140 04014001 (2014) DOI:10.1061/(ASCE)EM.1943-7889.0000691. [24] Rajendran, S., Spurgeon, S.K., Tsampardoukas, G., Hampson, R. Estimation of road frictional force and wheel slip for effective antilock braking system (ABS) control. Int J Robust Nonlin 29 736-765 (2019) DOI:10.1002/rnc.4366. [25] Han, T. Research on Intelligent Tire Vertical/Longitudinal Force Estimation Algorithm Based on Flexible Ring Model. PhD Thesis, Jiangsu University, Zhenjiang (2020) [26] Gustafsson, F. Slip-based tire-road friction estimation. Automatica 33 1087-1099 (1997) DOI:10.1016/S0005-1098(97)00003-4. [27] Yu, M., Wu, G., Kong, L., Tang, Y. Tire-pavement friction characteristics with elastic properties of asphalt pavements. Appl Sci-Basel 7 1123 (2017) DOI:10.3390/ app7111123. [28] Gerthoffert, J., Cerezo, V., Thiery, M., Bouteldja, M., Do, M.T. A brush-based approach for modelling runway friction assessment device. Int J Pavement En 21 1694-1702 (2020) DOI:10.1080/10298436.2018.1563786. [29] Sheng, S.G. Research on Stability Control of Distributed Drive Pure Electric Vehicle Based on Tyre Force Observation, MSc Thesis Jilin University, Jilin (2023) Acknowledgements This study was financially supported by the National Natural Science Foundation of China (52072156,52272366) and the Postdoctoral Foundation of China (2020M682269). Received 2024-05-12, revised: 2024-11-13, 2025-01-06, 2025-02-07, accepted: 2025-04-01. Author Contribution Jintao Zhang: writing – original draft, visualization, validation, methodology. Zhecheng Jing: investigation, resources, conceptualization, data curation. Haichao Zhou: writing – original draft, methodology, investigation, conceptualization. Yu Zhang: conceptualization. Guolin Wang: validation, supervision, conceptualization. Data Availability The data supporting the findings of this study are included in the article. Metoda dolocanja koeficienta oprijema med pnevmatiko in cestišcem na podlagi fizikalnega modela pnevmatike in deformacijskega signala pri cistem vzdolžnem zdrsu Povzetek Za natancen izracun koeficienta oprijema med kotalno pnevmatiko in cestišcem pri razlicnih stopnjah zdrsa ter izboljšanje varnosti in stabilnosti vozila je bila predlagana metoda ocene koeficienta oprijema med pnevmatiko in cestišcem, ki temelji na deformacijskih senzorjih in šcetkastih (brush) modelih. Najprej je bil izdelan FEM model radialne pnevmatike dimenzij 205/55R16, katerega ucinkovitost je bila potrjena s poskusi staticnega stika s podlago in merjenjem radialne togosti. Nato je bil med zaviranjem pridobljen signal obodne deformacije vzdolž sredinske linije notranje površine pnevmatike, dolžina kontaktne površine pa je bila ocenjena z uporabo formule za dolžino loka na podlagi povprecne razdalje med vrhovi prvega in drugega reda krivulj obodne deformacije. Simulacija zaviranja kotalnih pnevmatik je potrdila nesimetricnost porazdelitve tlaka na kontaktni površini. Položaj tocke zdrsa znotraj kontaktne površine je bil ocenjen na podlagi poljubne funkcije porazdelitve tlaka in šcetkastega modela. Za dolocitev ocenjevalne funkcije tocke zdrsa pri razlicnih stopnjah zdrsa je bila uporabljena nelinearna regresija. Vzpostavljena je bila funkcionalna zveza med koeficientom oprijema pnevmatika–cestišce in stopnjo zdrsa, ki upošteva znacilnosti trenja med gumijasto površino pnevmatike in cestišcem, pri cemenr je bil uporabljen model trenja na osnovi eksponentnega upadanja. Rezultati kažejo, da opisana metoda omogoca oceno koeficienta oprijema med pnevmatiko in cestišcem pri razlicnih stopnjah zdrsa, kar zagotavlja dragocena spoznanja za aplikacije inteligentnih pnevmatik na podrocju nadzora dinamike vozil. Kljucne besede inteligentna pnevmatika, ocena koeficienta oprijema med pnevmatiko in cestišcem, tocka zdrsa, stopnja zdrsa, nelinearna regresija Tooth Contact Analysis of Involute Beveloid Gear Based on Higher-Order Curve Axial Modification Yongping Liu1,2 – Qi Chen1 – Changbin Dong1 1 Lanzhou University of Technology, School of Mechanical and Electrical Engineering, China 2 Lanzhou University of Technology, State Wenzhou Engineering Institute of Pump & Valve, China cameliu@163.com Abstract This study investigates the tooth flank contact characteristics of a beveloid gear pair through the lens of higher-order curve tooth modification of the involute beveloid gear. The machining coordinate system of the modified gear pair is established, and its tooth surface equations are derived based on the principle of gear meshing and coordinate transformation. In this context, a contact analysis of the modified gear is conducted, examining the impact of varying parameters on the contact trace and contact ellipses, as well as the implications for meshing characteristics in the presence of assembly errors. The findings indicate that the contact form of the high-order curve axial modification of the beveloid gear pair is point contact. Furthermore, the maximum modification magnitude and the order of the modification curve influence the meshing performance of the beveloid gear pair. Additionally, the beveloid gear pair demonstrates enhanced tolerance to the center distance and the axis crossed error, while exhibiting reduced tolerance to the axis intersected error. Keywords involute beveloid gear, higher order curve axial modification, tooth contact analysis, transmission error, assembly error Highlights . The model of involute beveloid gear pair with higher-order curve axial modification is constructed. . The contact paths and contact ellipses of the modified beveloid gear pair are derived. . The transmission error of the modified gear pair is calculated. . The transmission error difference of the modified gear pair with different assembly errors is analyzed. 1 INTRODUCTION Involute beveloid gears offer a number of advantages, including ease of manufacture, low cost and high accuracy. As a result, they are suitable for a range of applications, including fast ferries, all-drive automobiles, aerospace precision machinery and other instances of power transmission. However, the beveloid gear pair exhibits transmission error, uneven contact stress distribution, and other factors during operation, which result in vibration impact and a reduction in load-bearing capacity and transmission performance. Consequently, there is a need to enhance the durability of the gear pair's surface and improve its load-bearing transmission performance. The load-bearing capacity represents a crucial criterion for evaluating the performance of gear transmissions. The accuracy of the model of tooth contact stress and inter-tooth load distribution coefficient is a fundamental prerequisite for ensuring the precision of the calculated results of gear bearing contact mechanics [1]. Modification can effectively enhance the gear dynamic load change gradient, mitigate shock, vibration and noise, and thus improve the quality of gear transmission. Axial modification is one of the methods of tooth modification, which can also be termed micro-geometric design. It is an effective approach to augment such performance. The modelling of involute beveloid gears provides the foundation for related research. Chen et al. [2] proposed a method for machining a straight-toothed beveloid internal gear pair with tooth blanks parallel to the gear shaping tool. Additionally, they derived theoretical models for the beveloid internal gear and the beveloid external gear, and conducted a load-bearing contact analysis. Sun et al. [3] put forth a novel approach to modelling the tooth surface of involute beveloid gears, employing the theory of minimum potential energy in conjunction with the slicing method. They also investigated the impact of design parameters on the contact characteristics of parallel axis beveloid gears. Sentürk and Fetvacyi [4] developed a mathematical method for the prevention of root cuts on the model of beveloid gears. They also developed a mathematical method for the prevention of root cuts on the model of parallel shaft variable thickness gears. Furthermore, they developed a mathematical method for the prevention of root cuts generated on beveloid gear. The meshing stiffness is also a frequently studied topic in the digital modelling of gears. Wen et al. [5] derived the contact line equations of a parallelled beveloid gear pair, proposed an analytical algorithm for calculating the meshing stiffness of a beveloid gear based on the slicing method, and analyzed the effect of changing parameters on the meshing stiffness of a beveloid gear. Song et al. [6] put forth a methodology for calculating the meshing stiffness of parallelled beveloid gears based on the potential method. They also investigated the impact of parameters such as pressure angle, pitch cone angle, gear displacement coefficient, and others. Zhou et al. [7] developed an alternative meshing stiffness model that considered the influence of parameters like the direction of inter-tooth friction. The influence of specific parameters on the meshing stiffness was examined, including pressure angle, pitch cone angle, and gear displacement coefficient. Mao et al. [8] enhanced the existing Weber energy method, which is based on the gear slicing method, constructed a time-varying mesh stiffness solution model for involute beveloid gears, and calculated and verified the time-varying meshing stiffness of beveloid gears. Liu et al. [9] constructed a three-dimensional model of a beveloid gear transmission and a gear dynamics model of a parallelled beveloid gear, based on the processing principle of a beveloid gear. They then proceeded to analyse the influence factors of the transmission error of a beveloid gear. The microgeometric design of the beveloid gear, also referred to as the modification of the beveloid gear tooth, has constituted a significant field of investigation in recent years. Wang et al. [10] put forth a methodology for calculating modifications to tooth profiles for the purpose of analysing tooth contact. Fuentes et al. [11] proposed an enhanced solution for intersecting beveloid gears with two distinct types of tooth profile bulge modification, with the objective of enhancing load-carrying capacity and reducing noise and vibration response. Sentürk et al. [12] developed a computer programme to obtain generating and generated surfaces of beveloid gears with modification. Morikawa et al. [13] examined the impact of tooth profile modifications on tooth surface damage through the utilization of a tooth fatigue test. In a further contribution to the field, Brecher et al. [14] put forward a method for the design of the tooth surface of intersected beveloid gears, taking into account the tolerance field for functional tolerance. Brimmers et al. [15] also made a significant contribution to the design of tooth flanks of variable thickness gears, investigating the possibility of free flank modification by means of a weighted objective function of beveloid gears with respect to the operating behavior. Ni et al. [16] designed a rack cutter with a parabolic modification to enhance the contact characteristics of helical beveloid gears. They then proceeded to investigate the impact of the parabolic modification on the contact path and contact ratio. Finally, they conducted an in-depth analysis to determine the sensitivity of beveloid gears to mounting errors. Liu et al. [17] put forth a numerical design methodology to augment the meshing characteristics of alternated beveloid gears through modification on the contact ellipse, contact path, transmission error, and relative curvature. Cao et al. [18] introduced the rack cutter with parabolic modification into the design of intersected beveloid gear pairs and investigated the effects of the modification coefficients on the mesh characteristics of the gear pair, including contact mode, transmission error, mesh stiffness, contact force, and tooth root stress. Zhang et al. [19] proposed a differentiated modification method based on a sine function, analyzed the dynamic characteristics and noise of gear pairs, comparing unmodified, normally, and differentially modified gears. Presently, the majority of modification designs for involute beveloid gears are based on previous design experience. Consequently, the innovative modification method for beveloid gears is of great significance, as it offers a promising avenue for enhancing the meshing characteristics of beveloid gear pair. The involute parallelled beveloid gear is distinguished by a distinctive tooth surface configuration, which exhibits the phenomenon of automatic backlash adjustment. However, modifications to this configuration are not recommended for significant alterations in the central portion of the tooth surface. Instead, high-order curves are employed to the axial modification, resulting in a notable discrepancy between the two ends and a relatively minor discrepancy in the central region. This paper derives the mathematical equations of the tooth surface of the beveloid gear under the corresponding design parameters, including the maximum amount of thickened gears. It clearly outlines the design steps of the modelling method of the high-order curve axial modification model and carries out a tooth contact analysis (TCA) for the modification of the parallelled beveloid gear. This aims to establishing a perfect modification theory of the involute parallelled beveloid gear transmission and promoting the development and application of the beveloid gear. 2 METHODS & MATERIALS The tooth profiles of the involute beveloid gear, as they exist in their unmodified state, can be considered to be generated under the tooth profile envelope of a usual rack cutter. In accordance with this principle, the tooth surface equations of the beveloid gear can be derived by undertaking a series of coordinate transformations through the rack cutter surface equations, as illustrated in Fig. 1. In the figure, the angle formed by the normal coordinate system Sn(xn, yn, zn) and the end coordinate system So(xo, yo, zo) of the rack cutter is the helix angle ß, while the angle formed by the end coordinate system So and the pitch coordinate system Sc(xc, yc, zc) is the pitch angle .. Coordinate system Sj(xj, yj, zj) is a movable coordinate system attached to the blank of the beveloid gear, while coordinate system Sb(xb, yb, zb) is a stationary coordinate system with its origin located at the center of the end face of the blank of the beveloid gear. The coordinate origin position On of Sn is jointly determined by the variable u = OoOn and the helix angle ß. The angular velocity during the rotation of the beveloid gear blank is ., the radius of the dividing cylinder of the beveloid gear is r, the translation speed of the rack cutter is v = .·r, and the real-time rotation angle of the beveloid gear during the machining process is fj. Fig. 1. Coordinate system for the generation of the beveloid gear In the normal coordinate system Sn of the rack cutter, the normal tooth profile vector of the rack cutter is given by r = [xyz 1]T (1) n nnn The rack cutter surface in the pitch coordinate system, designated Sc, can be obtained by combining the normal tooth profile vectors with the following coordinate transformations: rc = Mcnnr = [ xc yc zc 1], (2) where rc is the position vector of the rack cutter surface in the coordinate system Sn; The transformation matrix, designated as Mcn, is employed for the purpose of effecting a change from the coordinate system Sn to the coordinate system Sc, and the transformation matrix is given by: . cos.°sin . sin . cos . sin . ucos . sin .. . 0 cos . sin . usin . . M ˜ . (3) cn .°sin .°cos . coscos.. u . c coscos . . 000 1 ˆ The position vector rc of the tooth face of the rack cutter in the section coordinate system Sc is combined with the relative motion relationship between the rack cutter and the gear blank during the machining process. Consequently, the equation representing the tooth face of the beveloid gear in coordinate system Sj is derived as: xj ˜ xc cos. j ° yc sin. j . r(cos. j .. j sin)j .. . .yj ˜ xc sin . j ° yc cos. j . r((sin . j .. j cos)j . rj : ˆ , (4) z ˜ z . jc . .. j ˜ nxccy ° nycyc / nxcr where xj, yj, zj represent the position vectors of the tooth surface of the beveloid gear, and nxc, nyc and nzc denote the unit normal vectors of the rack cutter surface. In accordance with the tooth equation, a computer program has been developed to generate the point set of the tooth profile of an involute beveloid gear. Figure 2 illustrates the point cloud model of the tooth profile of the beveloid gear with the parameters outlined in Table 1. Table 1. Basic parameters of involute beveloid gears Parameters Pinion Gear Normal modules [mm] 2.5 Center distance [mm] 80 Pressure angle [°] 20 Tooth number 27 37 Tooth width [mm] 26 24 Pitch angle [°]) 6 6 Fig. 2. Point cloud model of beveloid gear tooth surface In the event of a high-order curve axial modification being applied on the involute beveloid gear, the amount of modification at each location in the axial direction can be expressed as follows: 2 . . 1 ° and after modification, respectively. The maximum modification magnitude of the beveloid gear is represented by .max. 3 RESULTS AND DISCUSSION 3.1 Tooth Contact Analysis of Involute Beveloid Gear Based on High-Order Curve Axial Modification Following the axial modification of the involute beveloid gear, the tooth face of the beveloid gear assumes a drum-shaped configuration, and the equations governing this tooth face exhibit greater complexity compared to those of the ordinary involute gear. In order to study the meshing characteristics, a meshing coordinate system is established for the beveloid gear. The TCA mathematical model of the involute beveloid gear pair is then obtained by transforming the coordinate system so that the modification of the pinion and the gear can achieve the correct meshing under this coordinate system. In accordance with the tooth surface equation of the high-order curve axial modification of the beveloid gear outlined in the preceding section, the tooth position vectors rp and rg and normal vectors np and ng of the pinion and gear can be derived using the following expressions: r uv ˜ xi° yj° zk (,) . ppp ppp r uv ˜ xi° yj° zk . g(,)gg ggg . , (6) n uv ˜ ni ° nj° (,) nk p pp px pyp pz . ˆng(,)uvg g ˜ ngxi° ngyj° ngzk . where the functions xp, yp, zp, npx, npy, and npz are defined with respect to two independent variables, up and vp. Similarly, the functions xg, yg, zg, ngx, ngy, and ngz are defined with respect to two independent variables, ug and vg. The beveloid gear pair transmission introduces two new variables, designated as .p and .g, which represent the position angles of the pinion and the follower, respectively. These variables satisfy the conditions at a ratio of i: ˜˜/ °i. (7) Following the implementation of the high-order curve axial modification scheme, the gp configuration of tooth contact is ˜ u a m0 ˆ. 2 B . u .  am1, characterized by point contact. Furthermore, the two position vectors (5) and the two normal vectors of the two conjugate tooth surfaces at the where ami represents the modification factor of gear thickness modification (i = 0, 1), B denotes the length of the tooth width of the beveloid gear, and u signifies the distance from the end face in the axial direction at that specific point. Fig. 3. Beveloid gear with higher order curve axial modification Figure 3 illustrates the modification of the involute beveloid gear under conditions of high-order curve axial modification. The dotted and solid lines indicate the tooth surface of the beveloid gear before contact point are equal [20]. In consequence, the mathematical model for contact analysis of a beveloid gear pair is given by: r uv˜° r (,,˜ ) . (,, ) uv . pppp gggg . . (8) n uv˜° n (,,˜ ) (,, ) uv . pppp gggg ˆ The partial vectors in the x, y, and z directions corresponding to the position and normal vectors in Eq. (8) are equal. Furthermore, by coupling Eq. (7), four relatively independent sets of equations can be obtained. In the contact analysis calculation of the beveloid gear pair, the position angle .p of the pinion can be regarded as a known quantity and solved separately by substituting a number of values within the meshing range. This results in the original set of equations becoming composed of four independent equations with four independent variables. Consequently, the solution to the set of contact equations can be achieved. The initial meshing point is designated as the reference point for calculation, and the resulting tooth contact trajectory points are shown as discrete red points in Fig. 4. In order to analyze the contact area of the tooth surface of the beveloid gear pair, it is necessary to calculate the curvature of the tooth surface. The first-order and second-order partial derivatives of the position vector r can be derived from the first and second fundamental quantities of the surface, as follows: .°r(,)uv r ˜ .u °u. .°r(,)uv rv ˜ . °u. .°2r(,)uv .ruu ˜ 2, (9) °u. 2 .°r(,)uv r ˜ .uv °° uv . .°r(,) 2 uv r ˜ .vv ˆ°v2 E ˜ r2, F ˜ r°r, G ˜ r2, u uvv L ˜ rn° , M= rn° , N=rn° , (10) uu uvuu rr n= t × v , (11) rr t × v where the ru and rv represent the first-order partial derivatives of the position vector r with respect to u and v, respectively. Meanwhile, ruu, ruv and rvv denote the second-order partial derivatives of the position vector r. E, F, and G are the first fundamental quantities of the surface, and L, M, and N are the second fundamental quantities of the surface. The unit normal vector n at point (u, v) is given by Eq. 11. Fig. 4. Contact path for the tooth surface with high-order curve axial modification The variation of curvature at the tooth contact point can be obtained from the principal curvature, which in turn allows for the calculation of the size and location of the contact area. The tooth profile direction principal curvature (K1) and the tooth width direction principal curvature (K2) are solved by the Gaussian curvature (K) and the mean curvature (H). This process is described by Eqs. (12) to (15): LN ° M 2 K ˜ , (12) EG ° F2 LG ° 2MF . NE H ˜ , (13) 2(EG ° F2) K1 ˜ H ° H 2 . K , (14) K2 ˜ H ° H 2 ° K. (15) The tooth contact form of the beveloid gear with high-order curve axial modification is characterized by point contact. However, in actual conditions, this contact will extend into an elliptical contact area due to elastic deformation, as illustrated in Fig. 5. Fig. 5. Tooth contact area In Figure 5, a and b represent the lengths of the semi major axis and semi minor axis of the contact ellipses, respectively. The angle a denotes the angle between the normal vector z at the contact point and the principal direction e of the tangent plane at that point. The dimensions and orientation of both the semi major axis and semi minor axis of the contact ellipses can be inferred from the curvature characteristics of the tooth surface: ˜ a = , (16) 122 2 K . K . g . 2gg cos 2°. g .. 112 2 ˜ b = , (17) 122 2 K . K . g . 2gg cos 2°. g .. 112 2 1 g1 . g2 cos 2° ˜ = arccos (18) 22 2 g . 2gg cos 2°. g 112 2 K1 ° K2 . g . g cos 2.. (19) ˜˜ 12 Given that tooth contact analyses are typically conducted under light load conditions, the elastic deformation d in Eq. (16) to Eq. (17) is typically assumed to be 0.00025 inch, or 0.00632 mm. K. 1 and K. 2 are obtained by adding K1 and K2 for the pinion and the gear, g1 and g2 are obtained by subtracting K2 from K1 for the pinion and the gear, respectively. Figure 6 illustrates the variation in contact ellipses of the modified beveloid gear pair throughout the meshing process. Fig. 6. Tooth contact ellipses Figure 7 illustrates the variation in the area of the contact ellipses over time for different modification magnitudes. It can be observed that an increase in the magnitude of the modification results in a reduction in the contact ellipse area, which consequently gives rise to an augmentation in the contact stress experienced by the tooth surface. In conjunction with the tooth contact trajectory, the theoretical transmission error of the modified beveloid gear pair in the transmission process can be derived [21], as illustrated in Fig. 8. Figure 9 illustrates the cumulative transmission error of the beveloid gear pair following the superimposition of the transmission error of a single pair of teeth under varying modification curves. Figure 10 depicts the corresponding peak-to-peak transmission error. It can be observed that an increase in the order of the modification curve is associated with a reduction in the transmission error of the beveloid gear pair. This is accompanied by a decrease in the peak-to-peak value of the transmission error, which is conducive to the smooth transmission of the gear pair. 3.2 Influence of Assembly Error on Contact Characteristics of Involute Beveloid Gear with High-Order Curve Axial Modification In involute beveloid gear pair, three principal forms of assembly errors may occur: center distance error, shaft staggering error, and shaft intersection angle error. These are illustrated in Fig. 11. In the event of a solitary assembly error, the positive and negative assembly errors exert an opposing influence on the meshing path, while the impact of center distance errors is comparatively negligible. In the event of equality between the error values, the shaft intersection angle error exerts a more pronounced influence on the meshing path than the shaft staggering error. Fig. 10. Peak-to-peak value of the transmission error of beveloid gear pair In order to investigate the impact of assembly errors on the meshing characteristics of the beveloid gear pair, a series of assembly errors have been introduced into the calculation of tooth contact analysis, with the resulting contact paths and contact ellipses illustrated in Fig. 12. Figure 13 illustrates the variation in transmission error for the beveloid gear pair in the presence of a center distance error of 0.5 mm, a shaft staggering error of 0.3°, and a shaft intersection angle error of 0.3°, respectively. Figure 14 depicts the peak-to-peak value of transmission error for the aforementioned three cases. It can be observed that any assembly error results in an increase in the transmission error of the beveloid gear pair, with the peak-to-peak value also being affected. The largest increase in transmission error is caused by the presence of shaft intersection angle error, while the largest increase in the peak-to-peak value of transmission error is caused by the presence of shaft staggering error. The combination of the results of the aforementioned analyses indicates that the beveloid gear pair with high-order curve axial modification exhibits superior tolerance performance with regard to center distance error and shaft staggering error. Conversely, the tolerance performance with respect to shaft intersection angle error is relatively poor. Based on the aforementioned analysis, the transmission characteristics of the variable thickness gear pair were simulated using simulation software. The resulting transmission error, as depicted in Fig. 15, exhibits a trend that is largely consistent with the theoretical transmission error. This finding corroborates the accuracy and reliability of the analytical method presented in this study. The beveloid gear pair with high-order curve axial modification shows strong tolerance to center distance error and shaft staggering error, but it is relatively poorly tolerated for shaft intersection error. The shaft intersection error can cause the contact point of gear meshing to shift, leading to significant transmission errors and uneven load distribution. This not only affects the stability and accuracy of the system but also reduces the service life of the gears. Consequently, while beveloid gears with high-order curve axial modified can tolerate other forms of assembly errors to some extent, their sensitivity to shaft intersection errors necessitates more precise control and compensation strategies in practical engineering applications. 4 CONCLUSIONS This paper presents a modelling and analysis of a beveloid gear pair with high-order curve axial modification. The conclusions drawn from this analysis are as follows: a) center distance error; b) shaft staggering error; and c) shaft intersection angle error 1. Following the axial modification, the tooth surface of the beveloid gear pair exhibits point contact. The contact trace is concentrated in the middle of the tooth surface, and the contact ellipse gradually increases from the root to the top of the tooth, before decreasing with the increase of the maximum amount of modification. This affects the smoothness of transmission of the gear pair. 2. It can be observed that an increase in the order of the modification curve is associated with a reduction in the transmission error of the beveloid gear pair. Consequently, the peak-to-peak value of the transmission error is also diminished. 3. The transmission error of the beveloid gear pair and its peak­to-peak values are observed to increase in the presence of assembly error. Conversely, the effect is observed to be smaller in the presence of center distance error and shaft staggering error. The largest effect is observed to occur in the presence of shaft intersection angle error. References [1] Wang, C., Cui, H., Zhang, Q., Zhang, B. Theoretical research progress of gear meshing efficiency. J Univ Jinan (Sci. & Tech.) 29 229-235 (2015) DOI:10.13349/j. cnki.jdxbn.2015.03.014. [2] Chen, Q., Song, C., Zhu, C., Du, X., Ni, G. Manufacturing and contact characteristics analysis of internal straight beveloid gear pair. Mech Mach Theory 114 60-73 (2016) DOI:10.1016/j.mechmachtheory.2017.04.002. [3] Sun, R., Song, C., Zhu, C., Wang, Y., Liu, K. Numerical study on contact force of paralleled beveloid gears using minimum potential energy theory. J Strain Anal Eng Des 56 249-264 (2020) DOI:10.1177/0309324720936894. [4] Sentürk, B., Fetvaci, M. Modelling and undercutting analysis of beveloid gears. J Fac Eng Archit Gaz 35, 901-916 (2020) DOI:10.17341/gazimmfd.544038. [5] Wen, J., Yao, H., Yan, Q., You, B. Research on time-varying meshing stiffness of marine beveloid gear system. Mathematics 11, 47-74 (2023) DOI:10.3390/ math11234774. [6] Song, C., Zhou, S., Zhu, C., Yang, X., Li, Z., Sun, R. Modeling and analysis of mesh stiffness for straight beveloid gear with parallel axes based on potential energy method. J Adv Mech Des Syst Manuf 12, 1-14 (2018) DOI:10.1299/ jamdsm.2018jamdsm0122. [7] Zhou, C., Dong, X., Wang, H., Liu, Z. Time-varying mesh stiffness model of a modified gear-rack drive with tooth friction and wear. J Braz Soc Mech Sci 44 213 (2022) DOI:10.1007/S40430-022-03517-8. [8] Mao, H., Fu, L., Yu, G., Tupolev, V., Liu, W. Numerical calculation of meshing stiffness for beveloid gear with profile modification. Compute Integra Manuf Sys 28 526-535 (2022) DOI:10.13196/j.cims.2022.02.017. [9] Liu, Y., Ren, Z., Wei, Y., Ma, D., Wei, S. Error study of involute beveloid spur gear transmission. Mach Tool Hydraulic 51, 1-6 (2023) DOI:10.3969/j.issn.1001­3881.2023.04.001. [10] Wang, C., Wang, S., Wang, G. A calculation method of tooth profile modification for tooth contact analysis technology. J Braz Soc Mech Sci 40, 1-9 (2018) DOI:10.1007/s40430-018-1268-4. [11] Fuentes, A., Gonzalez-Perez, I., Hayasaka, K. Computerized design of conical involute gears with improved bearing contact and reduced noise and vibration. VDI Bericht 2108, 635-646 (2010). [12] Sentürk, B., Fetvaci, M. Modelling and undercutting analysis of beveloid gears. J Fac Eng Archit Gaz, 35, 901-916 (2020) DOI:10.17341/gazimmfd.544038. [13] Morikawa, K., Kumagai, K., Nagahara, M. lnfuence of tooth flank modfication on tooth surface damage of conical involute gears. Trans JSME 81 1-8 (2015) DOI:10.1299/transjsme.15-00311. [14] Brecher, C., Löpenhaus, C., Brimmers, J. Function-oriented torlerancing of tooth flank modifications of beveloid gears. Procedia CIRP 43 124-129 (2016) DOI:10.1016/j.procir.2016.02.155. [15] Brimmers, J., Brecher, C., Löpenhaus, C. Potenzial von freien Flanken modifikationen für Beveloidverzahnungen. Forschung im Ingenieurwesen 81 83­94 (2017) DOI:10.1007/s10010-017-0232-2. [16] Ni, G., Zhu, C., Song, C., Du, X., Zhou, Y. Tooth contact analysis of crossed beveloid gear transmission with parabolic modification. Mech Mach Theory 113 40-52 (2017) DOI:10.1016/j.mechmachtheory.2017.03.004. [17] Liu, S., Song, C., Zhu, C., Ni, G. Effects of tooth modifications on mesh characteristics of crossed beveloid gear pair with small shaft angle. Mech Mach Theory 119 142-160 (2018) DOI:10.1016/j.mechmachtheory.2017.09.007. [18] Cao, Y., Ni, G., Fang, Z., Liu, Z. Effects of cutter modification on the meshing characteristics of intersected beveloid gear and involute cylindrical gear. J Phys Conf Ser 2691 012-022 (2024) DOI:10.1088/1742-6596/2691/1/012022. [19] Zhang, Y., Zhou, H., Duan, C., Wang, Z., Luo, H. Gear differential flank modification design method for low noise. Stroj Vestn-J Mech E 70 569-581 (2024) DOI:10.5545/sv-jme.2024.1072. [20] Wang, X., Lu, J., Yang, S. Sensitivity analysis and optimization design of hypoid gears’ contact pattern to misalignments. J Zhejiang Univ-Sc A 20 411-430 (2019) DOI:10.1631/jzus.A1900021. [21] Dong, C., Yang, X., Li, D., Zhao, G., Liu, Y. Service Performance optimization and experimental study of a new W-W type non-circular planetary gear train. Stroj Vestn-J Mech E 70 128-140 (2024) DOI:10.5545/sv-jme.2023.673. Acknowledgements This research was funded by the National Natural Science Foundation of China No. 52265008; the Youth Science Foundation of Gansu Province, No. 23JRRA751; the Education Department of Gansu Province: University Teachers Innovation Fund No. 2023A-021; the Wenzhou Municipal Basic Scientific Research Project of China No. 2022G0145. Received: 2024-12-23, revised: 2025-03-14, accepted: 2025-04-09 as Original Scientific Paper. Data Availability The data that support the findings of this study are available from the corresponding author upon reasonable request. Author Contribution Yongping Liu: Data curation, Qi Chen: Writing—original draft preparation, Changbin Dong: Writing—review & editing. Analiza zobnega kontakta evolventnega stožcastega zobnika na osnovi aksialne modifikacije s krivuljo višjega reda Povzetek V tej raziskavi so preucene kontaktne znacilnosti bokov zob evolventnega stožcastega zobniškega para z vidika aksialne modifikacije zobne površine s krivuljo višjega reda. Postavljen je koordinatni sistem obdelave modificiranega zobniškega para, izpeljane pa so tudi enacbe zobnih površin, ki temeljijo na principu zobniškega ubiranja in koordinatnih transformacijah. Izvedena je analiza kontakta modificiranega zobnika, pri cemer so preuceni vplivi spreminjajocih se parametrov na kontaktno sled in kontaktne elipse ter posledice za lastnosti ubiranja ob prisotnosti montažnih napak. Ugotovitve kažejo, da je oblika kontakta pri aksialni modifikaciji s krivuljo višjega reda pri stožcastem zobniškem paru v obliki tocke. Poleg tega maksimalna velikost modifikacije in red krivulje modifikacije vplivata na lastnosti ubiranja stožcastega zobniškega para. Stožcasti zobniški par izkazuje izboljšano toleranco glede napake medosne razdalje in napake križanja osi, obenem pa kaže manjšo toleranco do napake presecišca osi. Kljucne besede evolventni stožcasti zobnik, aksialna modifikacija s krivuljo višjega reda, analiza zobnega kontakta, napaka prenosa, napaka pri montaži Geometric Design Method of Lightweight Line Gear Mechanism Chao He1 – Yangzhi Chen1,2 – Xiaoxiao Ping1 – Zhen Chen1 – Maoxi Zheng2,3 – Qinsong Zhang2 1 Guangdong Ocean University, School of Mechanical and Energy Engineering, China 2 South China University of Technology, School of Mechanical and Automotive Engineering, China 3 University of Electronic Science and Technology of China Zhongshan Institute, School of Mechanical and Electrical Engineering, China xiaoping.xiao@foxmail.com Abstract Based on the space curve meshing theory, a lightweight line gear mechanism (LLGM) was proposed in this paper. The LLGM can achieve the objective of lightweight design by reducing the radial dimension of the gear, which has an improvement in terms of size reduction. The outstanding advantage of the LLGM is its ability to achieve a high transmission ratio. Three aspects were proposed to design the LLGM: first, the primary design method was obtained; second, an approach to simply and effectively establish the analytical model of the LLGM was presented, which showed that the lightweight characteristic of the LLGM is mainly reflected in the radial dimension; third, the discriminant condition of the LLGM was built. The simplicity and effectiveness of the LLGM were demonstrated by a design example, while gear contact simulation and kinematics experiments were carried out to verify the theoretical basis. The lightweight design method proposed in this paper belonged to structural lightweight design, which can effectively solve the lightweight design issue in gear transmissions characterized by light loads and high transmission ratios. Keywords lightweight design, space curve meshing, gear design, kinematics experiments, high transmission ratio Highlights . A novel lightweight line gear mechanism (LLGM) was proposed. . The outstanding advantage of the LLGM is the high transmission ratio. . The LLGM contributes to material amount savings. . Kinematics experiments have been carried out on a homemade test rig. 1 INTRODUCTION Gear transmission is widely recognized for its high meshing efficiency and reliability [1-4]. With the explosive growth of electromechanical products, the application of gears has become increasingly widespread [5]. Green transmission has emerged as a mainstream trend in the development of gears for the future. Achieving energy conservation and emission reduction has become a critical issue in the gear industry [6]. Among various solutions, structural lightweight design of gears has drawn the attention of many researchers and engineers [7-9]. Additionally, the miniaturization of electromechanical products has imposed higher requirements for lightweighting [10,11]. At present, lightweight gears have been important components in mini-electromechanical products, such as unmanned aerial vehicles and mini-robots [12,13]. Designing lightweight gears has been a key research issue for decades. In order to achieve the goal of lightweight, many design methods have been proposed. Some gears adopt a small module [14], others undergo topological optimization [15], and some are designed as specific gear devices [16,17]. Among the lightweight approaches, the simplest lightweight design method is to reduce the gear modulus and the number of teeth [18]. However, the load capacity of the gear decreases as the gear modulus decreases [19]. Besides, the teeth number of gears is normally not less than 17 to avoid undercut [20]. Therefore, the simplest lightweight gear design method is not universally effective. On the other hand, there are typically more than three gears when the transmission ratio of a gear device is more than 6 — a factor that obviously increases the difficulty of lightweight design due to the number of gears [21,22]. The lightweight gear pairs are also restricted by the radial dimension limitation of the two gears, that is, the diameter of the driven gear is approximately equal to the product of the diameter of the driving gear and the transmission ratio. Once the diameter of one gear and the value of the transmission ratio are determined, the radial dimension of the other gear cannot be reduced. In summary, the previous lightweight designs leave room for improvement in terms of size reduction, especially under high transmission ratios. To address the challenge of lightweight, this paper presents a lightweight design method based on a new theory. The method proposed in this paper belongs to structural design, and it can effectively address the challenge of lightweight design in mechanical transmissions with low-load and high-transmission­ratio characteristics. The research foundation of this paper primarily relies on the space curve meshing theory [23,24], called line gear (LG) [25]. Different from involute gears, LG can realize stable meshing transmission through two contact curves. Also, LG has the characteristics of a small number of teeth and no undercutting, making it suitable for lightweight design. In this paper, a novel LG, named lightweight line gear mechanism (LLGM) is proposed. The LLGM can break through the radial dimension limitation of the two gears and achieve the objective of lightweight design by directly reducing the radial dimension of the gear, thereby achieving better size reduction performance than previous lightweight gears and LG. The design of the LLGM was proposed from three aspects: first, the primary design method was obtained through research; second, a simple approach for establishing the LLGM was presented; third, the discriminant condition of the LLGM was built. After that, a design example of the LLGM was put forward and a comparison was carried out between the LLGM and an ordinary LG pair. Gear contact simulation was carried out to verify the contact characteristics. Finally, kinematics experiments ° sin  .ˆ 00 cos a a were carried out to verify the proposed design theory. The LLGM is .... ....  cos  sin 00 suitable for mini-electromechanical products, offering a new solution , (4) M ˜ a a 1a in mechanical transmissions characterized by high transmission 0 010 ratios, lightweight design and compact size. . 0 001 .  sin  ° .ˆ 00 cos b b .... .... 2 METHODS AND MATERIALS 2.1 Primary Design Parameters of the LLGM  sin 00 cos b b M ˜ , (5) 2b 0 010 . 0 001 . The LLGM is build in accordance to the space curve meshing theory. Therefore, the LLGM mainly involves the design of two conjugate curves that present cylindrical helices [25], as shown in Fig. 1. Two cylindrical helices denoted as A and B respectively, are rotated around their respective axes. Two fixed coordinate systems determine the relative position of curves A and B, and are named O1– x1y1z1 and O2– x2y2z2, respectively. The distance between them is expressed as a. The parameter equations of curves A and B were represented in the coordinate system Oa – xayaza and Ob – xbybzb, respectively. The O1– x1y1z1 and O2– x2y2z2 coincided with the Oa – xayaza and Ob – xbybzb at the initial position. The curve A rotated through an angle fa at a uniform angular velocity .a and the curve B rotated through an angle fb at a uniform angular velocity .b while the two curves conducted conjugate motion. Curves A and B can be expressed as Eqs. (1) and (2) in the coordinate system Oa – xayaza and Ob – xbybzb, respectively. .xM ..a ˜°ma cos t a ...a Ra ˜.yM ˜ma t sin, (1) ...a z ˜nt .M a ˆ ˆ.b ° xM ˜mb cos t where M21, M1a and M2b represent the transformation relationship between the coordinate systems O1– x1y1z1 and O2– x2y2z2, O1– x1y1z1 and Oa – xayaza, O2– x2y2z2 and Ob – xbybzb, respectively. The equation of curve A can also be expressed in the coordinate system O1– x1y1z1 by Eq. (6), which defines rotation of curve A around the z1 axis. R1 ˜ M ° Ra . (6) 11aa The equation of the curve B can be similarly expressed in the coordinate system O2– x2y2z2 by Eq. (7), defining rotation of curve B around the z2 axis. R2 ˜ M b ° Rb . (7) 22 b According to the conjugate condition, the curve A is tangent to the curve B while they are rotating, which can be expressed as: R2 = MR1. (8) 2 211 The two sides of Eq. (8) represented the equations of the curve A and B in the coordinate system O2– x2y2z2. The solution to Eq. (8) is the position of the contact points of the two curves. Further, Eq. (9) is obtained by substituting the known parameters into Eq. (8). ˜ in ˜ °. m b ˆ cos t .ˆ . . 00 ˆ.b cos s Rbb ˜ b b (2) ˜ b sin,t y m ° . . . ... 1 . . . ... .... .... M ˜˜ msint b si 00 . . . n cos b b zˆ.b ˜nbt  t M 0 010 n b wherein t is the independent variable, ma and mb represent the helix radius, na and nb represent the pitch parameter. .. 0 001 .° acost ˜˜ˆ .. ˆˆ The transformation matrix between the above coordinate systems m 100a in 00 cos s a a . . . ... . nta 011 . . . ... .... .... .... .... can be expressed as: ˜˜ masit 0100 0010 si 00 nn cos a a  . (9) ° . 100a 0 010 .... .... 0100 .. . 0001 00 , (3) M ˜ 21 0010 0001 . ˆ Fig. 1. Primary theory of the LLGM At the initial position, fa and fb are both equal to 0. At other positions, the relationship between fa and fb is expressed as Eq. (10) by the gear transmission principle. ˜a ° i˜b , (10) where i is the transmission ratio. By substituting the parameter i into Eq. (9), Eq. (11) can be obtained according to the space curve meshing theory. n i = b . (11) na Equation (11) reveals the main factor influencing the transmission ratio of the LLGM, is the pitch ratio of the two curves. 2.2 Lightweight Design Method of the LLGM Since the transmission ratio of the LLGM is directly related to the pitch ratio according to Eq. (11), the lightweight design of the LLGM can be mainly reflected in the radial dimension. A simple approach to establish the LLGM model was proposed for analysis, as shown in Fig. 2. In Figure 2, the LG tooth was generated by sweeping, where the cylindrical helix was used as the path and the tooth profile was as used as the profile. There were two gear pairs, one was the ordinary LG and the other was the LLGM. As Figure 2 shows, both of the two gear pairs have the same driving gear, but the driven gear of the LLGM is smaller than that of the ordinary LG, while they have the same transmission ratio. In other words, the LLGM can achieve lightweight design by reducing the radial dimension. 2.3 Lightweight Condition of the LLGM 2.3.1 The Discriminant Condition of the LLGM It can be seen from Fig. 2 that the LLGM and the ordinary LG have a similar structure. Therefore, it is necessary to introduce the discriminant condition of the LLGM. The main difference between the LLGM and the ordinary LG is the radial dimension. The change of the radial size brings a different sliding rate simultaneously. Therefore, the lightweight condition of the LLGM can be illustrated in terms of the sliding rate. According to the design theory of LG [26], the sliding rate can be calculated by Eq. (12). nb 2 .mb 2 Sr˜° , 1 (12) ina 2 .ma 2 where Sr represents the sliding rate. For the ordinary LG, the parameter mb is approximately equal to ima [26]. The relationship between the parameters na and nb can be expressed as Eq. (11). Therefore, the sliding rate of the ordinary LG can be derived from Eq. (12). Sro ˜ 0. (13) For the LLGM, the parameter mb is smaller than ima and the parameter nb is equal to ina. Therefore, nb 2 + mb 2 is smaller than in a 2 + ma 2 , and the sliding rate of the LLGM can be expressed as the inequality: Srl > 0, (14) where Srl is the sliding rate of the LLGM. As shown in Eq. (13) and inequality in Eq. (14), the sliding rates of the LLGM and the ordinary LG are different. In short, the discriminant condition of the LLGM can be derived as follows. 1. mb 0. 2.3.2 Degree of Lightweight of the LLGM The degree of lightweight of the LLGM can be described by comparing the radial sizes of the LLGM and the ordinary LG. The degree of lightweight of the LLGM can be expressed as: . ima . mb . Li ˜.100 °%, (15) m . im ˆ aa  The larger Li, the greater the degree of lightweight. According to Eq. (12), the sliding rate increases with the degree of lightweight. Generally, the degree of lightweight of the LLGM is limited by sliding speed, which can be calculated by: 22  n r . 22 nb ° mb . Sv ˜ˆ na ° ma ., (16) 30 ˆ i  .. where Sv is the sliding speed of the LLGM, and nr is the rotating speed of the driving wheel. The sliding speed can be selected according to the lubrication conditions [27]. Also, Eqs. (15) and (16) Fig. 2. Lightweight schematic of the LLGM demonstrate that the LLGM maintains a sliding ratio greater than 0, with the sliding rate increasing as the degree of lightweight increases. In contrast, the sliding rate of LG gears remains constant and theoretically approaches zero, while the sliding rate of involute gears continuously varies throughout the meshing process and reverses direction at the pitch point. In short, the basic design method of the LLGM can be derived as follows: 1. Provide the conjugate curves according to the lightweight requirements; 2. Design the two gears; 3. Check the sliding speed. 2.4 Experimental 2.4.1 Design Example of the LLGM An example of the LLGM was proposed to show the basic design method of the LLGM directly. For comparison, an ordinary LG pair was also proposed. The transmission ratios of the two gear pairs were both set to 8. The equations of the driving contact curves were set the same: .xM ..a ˜°10cos t a ...a Ra ˜.yM ˜10sin t , (17) ...a z ˜°12 7. t .M ˆ where Raa represents the equation of the driving contact curves. In this example, to achieve over 50 % lightweighting, the equations of the contact curves for the driven gears are as follows: .xM ..b ˜27 5.cos t b . ..b Rb ˜.yM ˜27 5.sin t , (18) ...b z ˜°101.6t .M ˆ .xM ..c ˜°80 cos t c . ..c Rc ˜.yM ˜80sin t , (19) ...c z ˜°101 6. t .M ˆ where Rbb represents the equation of the contact curves on the driven gear of the LLGM and Rcc represents the equation of the contact curves on the ordinary driven LG. Next, the tooth profile was set as a rounded equilateral triangle with a height of 2.5 mm and a fillet diameter of 1 mm; the tooth width was set to 40 mm; the tooth numbers of the ordinary LG pair were set to 4 (driving gear) and 32 (driven gear), while those of the LLGM were set to 2 (driving gear) and 16 (driven gear), respectively. The geometry of the LG resembles that of a worm, and among its parameters, the pitch parameter is the most critical. Specifically, the pitch values of the driving LG and the driven LG were set to 40 mm and 320 mm, respectively. The 3D models of the two gear pairs can be obtained according to the above settings. The main difference between the two gear pairs lies in the two driven gears, as shown in Fig. 3. The comparison between the ordinary LG pair and the LLGM was conducted, some comparative parameters were tabulated as shown in Table 1. Table 1. Comparisons between different gear pairs Comparative items The ordinary LG pair The LLGM Transmission ratio 8 8 Number of gears 2 2 Material amount [cm3] 681 120 Fig. 3. Model of the design example Fig. 4. The three reducers According to Eq. (15), the degree of lightweight of the LLGM design example was equal to 58.33 %. As shown in Table 1, the LLGM achieves significant lightweighting, saving 82.37 % material usage compared with the ordinary LG. 2.4.2 Testing of the LLGM Reducer For the convenience of analysis, an involute gear pair was also proposed. The modulus of the involute gear was set to 2 mm since the LG tooth size was equivalent to 2 mm modulus. For the involute gear, the tooth width was set to 40 mm; the tooth numbers were set to 10 (driving gear) and 80 (driven gear). The three reducers were designed on the basis of the above settings and then manufactured by 3D printing technology, as shown in Fig. 4. The three reducers (LLGM reducer, LG reducer and involute gear reducer) were successively installed on a homemade test rig for kinematics experiments, as shown in Fig. 5. Subsequently, kinematics experiments were carried out to test the performance of the three reducers on the homemade test rig, following the experimental methodology detailed in [28]. During the testing, the servo motor provided rotational speed to the reducer while the magnetic brake applied the load. The speed data were collected by the two encoders fixed on the input and output shafts, respectively. Fig. 5. Testing for the three reducers The input speed was set to 750 r/min, the sampling frequency of the encoder was set to 10 Hz, and the load was set to 1 N·m. 3 RESULTS AND DISCUSSION The experimental results were obtained after conducting the tests with the given experimental parameters. The transmission ratio data was obtained from the collected speed data. The transmission ratio diagram was drawn after using the six-point median filtering in data processing [25], as shown in Fig. 6. Then, according to the collected transmission ratio data, the error analysis was carried out, as shown in Table 2. Furthermore, simulation analysis was carried out for the three gear drives in ANSYS software, and the tooth surface contact stress results were shown in Fig. 7. Fig. 6. Transmission ratio diagram of the three reducers Table 2. Some transmission ratio errors of the three reducers Comparative items LLGM reducer LG reducer Involute gear reducer Average transmission ratio 8.006 7.988 8.031 Standard deviation 0.192 0.423 0.340 Range error 0.903 1.811 1.970 According to the simulation results in Fig. 7, both the LLGM and LG exhibit point contact, while involute gears feature line contact. The contact points of the LLGM are distributed along the tooth direction position parallel to the axis, and the trajectory of the meshing points is consistent with that shown in Fig. 2. The main experimental results (Fig. 6 and Table 2) indicate that the three reducers exhibit similar transmission ratio errors, demonstrating comparable kinematic performance. The results also show that the LLGM can realize stable transmission with a high transmission ratio, which explains the correctness of the geometric design method of the LLGM directly. During the testing, the input speed was given based on the rated speed of the motor. According to the Eq. (16), the sliding speed of the LLGM reducer was calculated as 0.314 nr, which can inform speed limitations or lubrication design. However, there were considerable transmission errors during the testing mainly due to manufacturing error. These errors could be mitigated through improved manufacturing precision. According to the simulation results, the contact stress of the LLGM is greater than that of the other two gears, which is unfavorable for its load bearing capacity. There are two key factors accounting for the relatively high contact stress. One is the presence of point contact, and the other is the lack of tooth profile optimization in the LLGM design. Despite using a simple tooth profile, we proved the geometric feasibility of its design. The simplicity and effectiveness of the LLGM were demonstrated by kinematics experiments and simulation analysis. Furthermore, floor cleaning robot gearbox composed of four-stage involute gear pairs was used for comparison with the LLGM reducer, as shown in Fig. 8. The main comparisons between the LLGM reducer and the floor cleaning robot gearbox were shown in Table 3. As shown in Table 3, the LLGM reducer achieves a 48% volume reduction and a 50 % reduction in the number of gears compared with the floor cleaning robot gearbox. Overall, the LLGM demonstrates significant advantages in lightweight design and high transmission ratio capability, which are especially suitable for small-scale electromechanical products. Table 3. Comparisons between the LLGM reducer and the floor cleaning robot gearbox Comparative items The LLGM reducer The floor cleaning robot gearbox Modulus [mm] 0.5 0.5 Transmission ratio 64 62.24 Number of gears 4 8 Overall size [mm] 43 × 28 ×21 85 × 36 ×16 Fig. 7. Simulation analysis of tooth surface contact conditions Fig. 8. The floor cleaning robot gearbox 4 CONCLUSIONS In this paper, a lightweight line gear mechanism (LLGM) was proposed, and the main work is summarized as follows: 1. Based on the theory of space curve meshing, the geometric design method of the LLGM was studied, while the lightweight objective was achieved by directly reducing the radial dimensions of gears. 2. An approach to establish a simple and effective analytical model of the LLGM was presented. The outstanding advantage of the mechanism is its ability to obtain a high transmission ratio in a compact size. 3. An LLGM example was given, and experimental tests were conducted. The results proved that the LLGM reducer can achieve as stable transmission as other gear transmissions under light loads. It has been verified that the geometric design method of the LLGM is correct and that the methodology can provide new solutions for the lightweight design mechanical transmissions. However, numerous challenges remain; such as torque measurement, lubrication optimization, strength characterization, dynamic response analysis, and precision manufacturing requirement. We aim to carry out more in-depth studies in the future to solve these problems. Nomenclature [12] Zhao, X., Zhou, Z., Zhu, X., Guo, A. Design of a hand-launched solar-powered O1– x1y1z1 Reference coordinate system for driving gear O2– x2y2z2 Reference coordinate system for driven gear Oa – xayaza Body coordinate system of driving gear Ob – xbybzb Body coordinate system of driven gear .a Uniform angular velocity of driving gear, [rad /s] .b Uniform angular velocity of driven gear, [rad /s] ma Helix radius of driving space curve, [mm] mb Helix radius of driven space curve, [mm] na Parameters related to driving space curve pitch, [mm] nb Parameters related to driven space curve pitch, [mm] i Transmission ratio nr Rotating speed, [r/min] t Independent variable fa Rotation angle of driving gear, [rad] fb Rotation angle of driven gear, [rad] M21 Transformation matrix between O2– x2y2z2 and O1– x1y1z1 M1a Transformation matrix between O1– x1y1z1 and Oa – xayaza M2a Transformation matrix between O2– x2y2z2 and Ob – xbybzb Raa Equation of driving space curve in Oa – xayaza Rbb Equation of driven space curve in Ob – xbybzb R11 Equation of driving space curve in O1– x1y1z1 R22 Equation of driven space curve in O2– x2y2z2 Sro Sliding rate of the ordinary LG Sr Sliding rate Srl Sliding rate of the LLGM Li Degree of lightweight Sv Sliding speed, [m/s] Rc Equation of driven space curve for ordinary LG References [1] Wang, C. Study on dynamic performance and optimal design for differential gear train in wind turbine gearbox. Renew Energ 221, 119776 (2024) DOI:10.1016/j. renene.2023.119776. [2] Okorn, I., Nagode, M., Klemenc, J. Operating performance of external non-involute spur and helical gears: A review. Stroj vestn-J Mech E 67 (2021) DOI:5545/sv­jme.2020.7094. [3] Zhang, Y., Zhou, H., Duan, C., Wang, Z., Luo, H. Gear differential flank modification design method for low noise. Stroj vestn-J Mech E 70 569-581 (2024) DOI:10.5545/sv-jme.2024.1072. [4] Dong, C., Yang, X., Li, D., Zhao, G., Liu, Y. Service performance optimization and experimental study of a new WW type non-circular planetary gear train. Stroj vestn-J Mech E 70 128-140 Stroj vestn-J Mech E 70 128-140 (2024) DOI:10.5545/sv-jme.2023.673. [5] Litvin, F.L., Fuentes, A. Gear geometry and Applied Theory. Cambridge University Press (2004) New York, DOI:10.1017/CBO9780511547126. [6] Pingkuo, L., Huan, P. What drives the green and low-carbon energy transition in China?: An empirical analysis based on a novel framework. Energy 239 122450 (2022) DOI:10.1016/j.energy.2021.122450. [7] Cosco, F., Adduci, R., Muzzi, L., Rezayat, A., Mundo, D. Multiobjective design optimization of lightweight gears. Machines 10 779 (2022) DOI:10.3390/ machines10090779. [8] Ide, T., Kitajima, H., Otomori, M., Leiva, J.P., Watson, B.C. Structural optimization methods of nonlinear static analysis with contact and its application to design lightweight gear box of automatic transmission of vehicles. Struct Multidiscip O 53 1383-1394 (2016) DOI:10.1007/s00158-015-1369-y. [9] Wang, D., Shen, H., Dong, H., Yu, S. A novel approach for conceptual structural design of gearbox. Mechanisms, Transmissions and Applications: Proceedings of the 3rd MeTrApp Conference (2015) 177-185, DOI:10.1007/978-3-319-17067­1_19. [10] Tezel, T., Schultheiss, U., Hornberger, H., Kovan, V. Operational wear behaviour of 3D-printed lightweight metal gears: EDS and oil analysis comparison. Mater Test 66 830-834 (2024) DOI:10.1515/mt-2023-0222. [11] Wu, J., Wei, P., Zhu, C., Zhang, P., Liu, H. Development and application of high strength gears. Int J Adv Manuf Tech 132 3123-3148 (2024) DOI:10.1007/ s00170-024-13479-x. unmanned aerial vehicle (UAV) system for plateau. Appl Scie-Basel 10 1300 (2020) DOI:10.3390/app10041300. [13] Zheng, C., Lee, K. WheeLeR: Wheel-leg reconfigurable mechanism with passive gears for mobile robot applications. 2019 International Conference on Robotics and Automation (2019) 9292-9298 DOI:10.1109/ICRA.2019.8793686. [14] Concli, F. Tooth root bending strength of gears: dimensional effect for small gears having a module below 5 mm. Appl Sci 11 2416 (2021) DOI:10.3390/ app11052416. [15] Ramadani, R., Pal, S., Kegl, M., Predan, J., Drstvenšek, I., Pehan, S., Belšak, A. Topology optimization and additive manufacturing in producing lightweight and low vibration gear body. Int J Adv Manuf Tech 113 3389-3399 (2021) DOI:10.1007/s00170-021-06841-w. [16] Chang, J., Chang, I. Design of high ratio gear reducer using vernier differential theory. J Mech Design, 143 023401 (2021) DOI:10.1115/1.4047347. [17] Mo, J., Gong, X., Luo, S., Fu, S., Chang, X., Liao, L. Geometric design and dynamic characteristics of a novel abnormal cycloidal gear reducer. Adv Mech Eng 15 16878132231158895 (2023) DOI:10.1177/16878132231158895. [18] Hsueh-Cheng, Y., Huang, Z.W. Using variable modulus to modify a rack cutter and generate a gear pair. P I Mech Eng C-J Mec 236 5342-5359 (2022) DOI:1177/09544062211061711. [19] Mo, S., Zhou, C., Zhang, W., Liu, Z., Hu, X., Bao, H., et al. Study on small modulus modified gear damage fault and electromechanical coupling vibration characteristics under multiple working conditions. P I Mech Eng C-J Mec 238 3746-3774 (2024) DOI:1177/09544062231207440. [20] Alipiev, O., Antonov, S., Grozeva, T. Generalized model of undercutting of involute spur gears generated by rack-cutters. Mech Mach Theory 64 39-52 (2013) DOI:1016/j.mechmachtheory.2013.01.012. [21] Xin, W., Zhang, Y., Fu, Y., Yang, W., Zheng, H. A multi-objective optimization design approach of large mining planetary gear reducer. Sci Rep-UK 13 18640 (2023) DOI:1038/s41598-023-45745-5. [22] Sun, G.B., Chiu, Y.J., Zuo, W.Y., Zhou, S., Gan, J.C., Li, Y. Transmission ratio optimization of two-speed gearbox in battery electric passenger vehicles. Adv Mech Eng 13 16878140211022869 (2021) DOI:1177/16878140211022869. [23] Chen, Z., Chen, Y.Z., Ding, J. A generalized space curve meshing equation for arbitrary intersecting gear. P I Mech Eng C-J Mec 227 1599-1607 (2013) DOI:1177/0954406212463310. [24] Chen, Y.Z., He, C., Lyu, Y.L. Basic theory and design method of variable shaft angle line gear mechanism. Stroj vestn-J Mech E 67 352-362 (2021) DOI:5545/ sv-jme.2021.7216. [25] Chen, Y.Z. Huang, H., Lv, Y. A variable-ratio line gear mechanism. Mech Mach Theory 98 151-163 (2016) DOI:1016/j.mechmachtheory.2015.12.005. [26] Chen, Y., Li, Z., Xie, X., Lyu, Y. Design methodology for coplanar axes line gear with controllable sliding rate. Stroj vestn-J Mech E 64 (2018) DOI:5545/sv­jme.2017.5110. [27] Chen, Y., Chen, H., Lyu, Y., Zhu, A. Study on geometric contact model of elastohydrodynamic lubrication of arbitrary intersecting line gear. Lubr Sci 30 387-400 (2018) DOI:1002/ls.1429. [28] He, C., Chen, Y., Lin, W., Lyu, Y. Research on the centre distance separability of parallel shaft line gear pair. P I Mech Eng C-J Mec 236 3261-3274 (2022) DOI:10.1177/09544062211035809. Acknowledgements This research was supported by program for Zhanjiang Science and Technology Plan Project (No.2024B01054 and 2024B01042); and scientific research start-up funds of Guangdong Ocean University (Grant No. YJR23010, YJR23002, YJR22016, and YJR24019); National Natural Science Foundation of China (Grant No. 52275073 and 52405056); and Open Fund of Guangdong Provincial Key Laboratory of Precision Gear Flexible Manufacturing Equipment Technology Enterprises, Zhongshan MLTOR Numerical Control Technology Co., LTD & South China University of Technology (No. 2021B1212050012-08). Received: 2024-10-04, revised: 2025-03-15, accepted: 2025-04-10 as Original Scientific Paper Data Availability The data supporting the findings of this study are included in the article. Author Contribution Chao He: Methodology, Formal Aanalysis, Validation, Writing – original draft. Yangzhi Chen: Conceptualization, Supervision. Xiaoxiao Ping: Conceptualization, Methodology, Resources, Formal Analysis. Zhen Chen: Writing – review & editing, Supervision. Qinsong Zhang: Formal analysis, Investigation, Writing – review & editing. Geometrijska metoda zasnove lahkega linijskega zobniškega mehanizma Povzetek Na osnovi teorije ozobljenja prostorskih krivulj je v tem clanku predstavljen lahek linijski zobniški mehanizem (LLZM). LLZM omogoca dosego lahke konstrukcije z neposrednim zmanjšanjem radialne mere zobnika, kar zagotavlja pomembno izboljšavo glede zmanjšanja velikosti. Posebna prednost LLZM je možnost doseganja visokih prenosnih razmerij. Za zasnovo LLZM so bili predlagani trije vidiki: prvic, dolocena je bila osnovna metoda zasnove; drugic, predstavljena je bila enostavna in ucinkovita metoda vzpostavitve analiticnega modela LLZM, ki je pokazala, da se lahka zasnova LLZM odraža predvsem v zmanjšanju radialnih mer; tretjic, postavljen je bil kriterij za presojo ustreznosti LLZM. Preprostost in ucinkovitost LLZM sta bila prikazana s konkretnim primerom zasnove ter preverjena s simulacijo zobniškega stika in kinematicnimi poskusi, ki potrjujejo teoreticne osnove. Metoda zasnove, predstavljena v tem clanku, spada med strukturne metode lahke zasnove in ucinkovito rešuje problem lahkih zobniških prenosov, ki prenašajo majhne obremenitve in zagotavljajo visoka prenosna razmerja. Kljucne besede lahka zasnova, ozobljenje prostorske krivulje, zasnova zobnikov, kinematicni poskusi, visoko prenosno razmerje Simulation Analysis and Experimental Study on Vibration Reduction Performance of Groove-Textured Friction Pair Surfaces Wusheng Tang – Yufei Nie – Zhuo Zhang – Wei Lin – YanKai Rong – Yaochen Shi – Ning Ding School of Mechanical and Vehicle Engineering, Changchun University, China 3765750755@qq.com Abstract To further investigate the influence of surface texture on the vibration characteristics of friction pair surfaces, this study fabricated groove textured on 45# steel surfaces using laser marking technology, forming a friction pair with chloroprene rubber. Numerical analysis of their friction processes under varying speeds was conducted via the finite element analysis software ABAQUS/Explicit (explicit dynamic solver).The results indicate that with increasing speed, the grooved-texture surface enhances its capability to reduce both the intensity and the continuity of contact stress concentration on the friction plate surface, while simultaneously accelerating the diffusion speed of contact stress from the leading edge to the trailing edge of the friction block, thereby better suppressing its concentration at the leading edge. Meanwhile, friction tests were conducted at varying speeds on the HCM-5 reciprocating friction and wear testing machine for the 45# steel-chloroprene rubber friction pair. The results show that at all speeds, the system damping of the freely decaying oscillation component on the surface of the groove texture after "groove-crossing fluctuations" is significantly increased compared to the non-textured surface. As the speed increases, the damping effect of the groove texture on the vibration of the friction pair surface gradually enhances, and the reduction in energy density at the main frequency of the friction pair surface becomes increasingly evident. This corresponds to the numerical analysis results, illustrating the influence of speed on the improvement of the vibration characteristics of the friction pair surface by the groove texture. Keywords groove textured, friction characteristics, numerical analysis, velocity Highlights . Groove textures reduce stress concentration at the friction block's front corner. . Groove edges disrupt continuous stress concentration on the friction plate surface. . As speed increases, textures enhance the reduction in stress concentration continuity and intensity. . Groove textures decrease vibration levels on friction surfaces. . Friction surfaces exhibit a dominant frequency of 1116 Hz, with energy density decreasing as speed rises. 1 INTRODUCTION Relative motion between friction pairs inevitably generates friction and wear. However, in practice, the failure of most mechanical components is not solely caused by friction and wear alone— friction-induced vibration also plays a critical role [1]. Such vibration not only accelerates component wear but may also lead to loosening, fracture, and even systemic instability. Therefore, it is imperative to prioritize vibration issues arising from friction and implement effective damping measures to enhance the stability and reliability of mechanical systems [2-4]. Over the past few decades, surface texturing has gained significant attention as an effective surface modification technique for vibration suppression and stress distribution optimization. By designing specific micro-scale geometric patterns on friction pair surfaces, surface texturing can substantially modify the tribological behavior at the contact interface, thereby reducing vibration, improving stress distribution, and enhancing overall system performance [5-7]. Currently, research on surface texturing primarily employs a combination of experimental studies, numerical simulations, and theoretical analysis. Experimental studies utilize equipment such as friction and wear testers and vibration testing platforms to directly measure the friction coefficient and vibration acceleration of textured surfaces, thereby validating their vibration-damping effects. For instance, experiments have investigated the influence of laser surface texturing on the tribological behavior of gray cast iron [8]. Other studies have experimentally demonstrated the potential of surface texturing in non-metallic sliding-element bearings, revealing its role in reducing friction, minimizing wear, and optimizing lubrication mechanisms [9]. Additionally, experimental research has examined the fretting wear behavior of texture Ti-6Al-4V alloy under oil lubrication [10]. Numerical simulations, including finite element analysis (FEA) and computational fluid dynamics (CFD), model the contact stress, friction force, and vibration characteristics of textured surfaces to uncover underlying mechanisms, often in combination with experimental studies. For example, simulations and experiments have been jointly used to study the friction-reducing effects of fan-shaped surface textures [11]. Similarly, laser-textured surfaces in reciprocating line contacts have been analyzed through combined numerical and experimental approaches to evaluate lubricant film thickness and friction forces [12]. Theoretical analysis involves establishing mathematical models to explore the influence of texture parameters (e.g., shape, size, distribution) on vibration reduction and stress distribution. For instance, theoretical studies have calculated frictional energy dissipation in mechanical contacts under random vibrations [13]. In terms of research objects, studies have primarily focused on metal-polymer friction pairs (e.g., steel-rubber, aluminum­polytetrafluoroethylene (PTFE)). For instance, research has investigated the effect of metal surface texture on the friction and wear performance of shale oil pumps in highly abrasive media [14], as well as the influence of lead-free bearing materials and surface textures on the performance of sliding bearings [15]. For metal-metal friction pairs (e.g., steel-steel, copper-aluminum), studies have examined the tribological behavior of gradient nanostructures on TB8 titanium alloy surfaces [16], the lubrication performance of femtosecond laser-textured stainless steel surfaces [17], and the tribological properties of single- and multi-shape laser-textured surfaces under dry friction conditions [18]. Study the concentration of stress and strain on gear tooth surfaces and observe the reduction in wear [19]. Metal-polymer pairs have garnered significant attention due to their widespread engineering applications (e.g., bearings, seals). Common surface texture types include micro-dimples, grooves, and wave-like patterns, with groove textures being a research hotspot due to their ease of fabrication and excellent vibration-damping effects. Additionally, the influence of texture parameters (e.g., depth, width, spacing) on vibration reduction and stress distribution remains a key research focus. In recent years, researchers have proposed various novel mechanisms for surface texturing to reduce vibration and improve stress distribution. For instance, the lubrication theory highlights that surface textures facilitate the formation of lubricating films, thereby reducing friction-induced vibration [20]. Additionally, synergistic mechanisms suggest that the combined effect of textures and lubricants can further enhance wear resistance. For example, the integration of laser texture with solid lubrication has been shown to significantly improve the tribological performance of titanium alloys [21]. Similarly, the synergistic interaction between self-lubricating materials and laser-textured surfaces—such as MoS2 and graphene coatings combined with laser texture—has demonstrated a remarkable enhancement in friction performance [22-23]. Tools with microstructures on both the front and rear cutting edges combined with solid lubricants can effectively improve cutting performance [24]. In summary, extensive research has demonstrated the critical role of surface texturing in reducing friction and wear across various tribological systems. However, existing studies have primarily focused on the dimensional parameters of textures, with limited investigation into how the vibration-damping performance of textured surfaces varies under different sliding velocities. To address this gap, this study employs laser processing to fabricate groove textures on 45# steel surfaces. Using a face-to-face contact configuration with chloroprene rubber (synchronous belt material) as the counter face, reciprocating friction tests are conducted on a HCM-5 tribometer at varying speeds to examine the influence of velocity on the vibration-suppression performance of textured surfaces. Furthermore, explicit dynamic simulations are performed using the finite element software ABAQUS/Explicit to numerically analyze the friction process and validate the experimental findings. This work lays a foundation for the future application of textured surfaces in synchronous belt tensioners. 2 METHODS AND MATERIALS This study employs a combined approach of numerical analysis and experimental research to enhance the accuracy and reliability of the investigation. Firstly, a simulation model is established using numerical analysis to explore the variation patterns of contact stress and friction force in the friction system under different speeds. Secondly, to validate the accuracy of the numerical analysis and investigate subtle differences in practical operation, multiple experiments were designed and conducted to record changes in vibrational acceleration at the friction interface in detail, while comparing and analyzing the results of numerical simulations and actual tests in real time. This approach not only enables timely identification and optimization of potential issues in the numerical model but also accumulates practical experience, thereby enhancing the research's practical value. 2.1 Simulation, Numerical Models The numerical model is designed to simplify the complex three-dimensional system of subsequent experiments into a relatively straightforward three-dimensional model. Its primary purpose is to predict general trends during the experimental process, such as variations in friction force and contact stress on the friction pair surface. The three-dimensional model comprises two main components, as shown in Fig. 1, a friction block and a friction plate (specific dimensional parameters are listed in Table 1). The textured friction plate features groove textured with a spacing of 4000 µm, width of 400 µm, and depth of 30 µm. The figure also illustrates the boundary conditions for the numerical analysis: at t = 0.1 s, a fixed Z-direction load of F = 400 N is applied to the friction block; five degrees of freedom of the friction block are constrained (except for Z-translation); five degrees of freedom of the friction plate are constrained (except for X-translation); a constant translational velocity v is applied along the X-direction; and the friction type is set to "hard" friction with a coefficient of 0.4. Table 1. Part Size Parameters Part L [mm] W [mm] H [mm] Friction block 20 20 5 Friction plate 200 50 10 Fig. 1. Numerical model of friction system 2.2 Experimental 2.2.1 Test Methods and Sample Preparation To systematically investigate the influence of surface textures on friction performance, this study conducts comparative experiments using multiple sets of data. Smooth surfaces are used as the control group, while textured surfaces serve as the experimental group. During the experiments, consistent texture parameters and loads are maintained, with a focus on comparing the effects of surface textures on friction performance at different speeds. The test conditions are a dry state under atmospheric conditions (temperature: 22 °C to 25 °C). To account for the randomness of the test process, each test group is repeated at least 5 times to ensure accuracy. The friction plates for both the control group and the experimental group are made of 45# steel, with a hardness of 650 HBS, elastic modulus E = 200 GPa, Poisson's ratio of 0.31, and density of Fig. 3. Friction plate 7850 kg/m³. The surfaces are polished to maintain smoothness. The control group has no texture processing, and the friction plate dimensions are 200 mm × 50 mm × 10 mm, serving as the baseline reference. The experimental group's friction plates are processed with groove textures using a laser marking machine, as shown in Fig. 2, with a depth of 30 µm, width of 400 µm, and spacing of 4000 µm. The actual processed friction plate is shown in Fig. 3. The counter friction blocks are made of chloroprene rubber, with an elastic modulus a) E = 6 MPa, Poisson's ratio of 0.4, density of 1240 kg/m³, and dimensions of 20 mm × 20 mm × 5 mm. 2.2.2 Test Equipment and Parameters Friction tests at varying speeds for the 45# steel-chloroprene rubber friction pair were conducted using the HCM-5 universal reciprocating friction and wear testing machine. A schematic diagram of the experimental system is shown in Fig. 4a. The testing machine, as shown in Fig. 4b, has a maximum test force of 5000 N, a maximum reciprocating stroke of 200 mm, and a reciprocating frequency range of 1 to 60 cycles/min. Additionally, a KISTLER 3D accelerometer (8763B100AT) is used, with a range of ±100 g, Integrated electronics piezo-electric (IEPE) output sensitivity of 50 mV/g ± 3 mV/g, and a response frequency of 0.5 Hz to 7 kHz, as shown in Fig. 4c. The accelerometer is fixed at 10 mm to 20 mm above the friction plate in the direction of motion. Signal acquisition is performed using the DEWESOFT 8-channel data acquisition system (X-DSA) with a sampling frequency of 10 kHz, as shown in Fig. 4d. During the tests, a constant load of 200 N is applied vertically to the friction block, while the friction plate undergoes reciprocating motion at different speeds: 2 mm/s, 4 mm/s, 6 mm/s, and 8 mm/s, with a reciprocating stroke of 100 mm. 3 RESULTS AND DISCUSSION 3.1 Simulation results: Contact Stress Analysis of Friction Sub-Surface The contact interface stress variations when the friction block is positioned at the center of the friction plate (t = 0.5 s) were analyzed, and the numerical results are presented in Fig. 5. The surface contact stress contour maps of the friction block at various speeds for smooth and textured friction pairs are shown in Figs. 5a and b, respectively. b) c) d) Fig. 4. Test equipment; a) schematic diagram of the experimental system, b) HCM-5 testing machine, c) KISTLER sensor, and d) X-DSA Data Acquisition System SV-JME . VOL 71 . NO 5-6 . Y 2025 . 209 It can be observed that the maximum contact stress in the smooth friction pair is primarily concentrated at the front corner and front edge of the friction block. In contrast, for the textured friction pair, the maximum contact stress is only concentrated at the front corner and shows a tendency to diffuse toward the rear corner, with minimal stress concentration at the front edge. Additionally, the maximum contact stress value decreases from 1.77×107 Pa to 7.97×105 Pa, a reduction of approximately 95.5 %. This demonstrates the role of groove textures in disrupting and dispersing interfacial contact stresses. Moreover, this effect becomes more pronounced as speed increases, as the diffusion of stress from the front to the rear corner accelerates. To more intuitively observe the effect of groove texture edges, Figs. 5c and d show the contact stress contour maps of smooth and textured friction plate surfaces at various speeds. Significant stress concentrations are observed on all smooth surfaces, but as speed increases, the continuity and intensity of these concentrations decrease noticeably, as seen at points 1, 2, and 3 in Fig 5. Although the concentration intensity increases slightly when speed rises from 2 mm/s to 4 mm/s, the stress concentration at point 4 has already disappeared. Overall, as speed increases, both the area and intensity of stress concentration on smooth surfaces gradually decrease. Similarly, as shown in Fig. 5d, the contact stress contour maps of the textured friction plate surface reveal that, compared to the smooth surface, the maximum contact stress decreases from 1.099×106 Pa to 4.885×105 Pa, a reduction of approximately 55.56 %, with an overall significant decrease in stress values. Furthermore, at various speeds, as indicated at points 5, 6, 7, and 8 in the Fig 5d, high-intensity stress concentrations are significantly reduced compared to the smooth surface, and this reduction trend is proportional to speed. This demonstrates that groove textures effectively reduce stress concentration and intensity at all speeds, and overall, the effect becomes more pronounced as speed increases. 3.2 Simulation results: Friction Analysis of Friction Sub-Surfaces To further observe the "disruptive" effect of groove textures, the friction force on the friction plate surface over time at various speeds (with a friction coefficient of 0.6 for smooth surfaces and 0.4 for textured surfaces) is obtained through numerical analysis, as shown in Figs 6a and b. Firstly, all figures exhibit fluctuations in friction force caused by collisions between the edges of the friction block and the groove textures (as shown in the magnified views at point 1 in each figure). At speeds of 2 mm/s to 8 mm/s, the fluctuation ranges are approximately ±10 N, ±20 N, ±23 N, and ±26 N, respectively. This is because, during the collision process, the system's momentum is conserved, meaning the total momentum of the system remains constant before and after the collision. The higher the speed, the greater the momentum of the object, resulting in a larger impact force during collisions and more pronounced fluctuations in friction force. However, as the friction block continues to slide, the friction force returns to its previous level due to the area loss caused by the grooves. It is precisely this intermittent friction state that achieves Fig. 5. Contact stress contour maps; a) smooth-friction block, b) textured friction block, c) smooth-friction plate, and d) textured-friction plate the goal of suppressing contact stress concentration on the friction pair surface, thereby preventing sustained vibrations in the friction system. Secondly, due to differences in speed, the frequency of edge collisions varies. Although higher speeds result in higher frequencies and more pronounced "disruptive" effects, they also imply greater impacts. Therefore, analyzing the vibrations of the friction system is essential. 3.3 Test Results: Time-Domain Analysis of Vibration Acceleration Due to the presence of groove edges, the friction force can momentarily surge, which may lead to increased damping in the system during motion, thereby reducing vibration acceleration. Conversely, a decrease in friction force may result in increased vibration acceleration. In reciprocating friction tests, the periodic variation of friction force can cause periodic fluctuations in the vibration acceleration of the friction surface, directly affecting the amplitude and frequency of these fluctuations. Therefore, analyzing the vibration acceleration of the friction surface is essential. The average vibration acceleration amplitudes in three directions for groove-textured friction surfaces at different relative speeds are shown in Fig. 7. It is evident that at all speeds, the average amplitude in the Z-direction (normal to the groove edges) is significantly larger than those in the Y-direction (parallel to the groove edges) and X-direction (perpendicular to the groove edges). This is because the Z-direction is the load application direction during friction, making the friction force changes most pronounced during impacts, leading to more surface vibrations. Furthermore, as shown in Fig. 7b, the Fig. 6. Friction time-domain diagram at different velocities: a) 2 mm/s, b) mm/s, c) 6 mm/s, and d) 8 mm/s variation trends of the average amplitudes in the three directions are highly consistent with speed changes. Therefore, this study primarily focuses on analyzing the Z-direction, where the vibration acceleration signal changes are most significant. a) b) Fig. 7. Mean value of vibration acceleration; a) three-dimensional, and b) line graphs The time-domain graph of the Z-direction vibration acceleration signals of the surface of the belt textured in two cycles at each speed is shown in Fig. 8a, which is processed by the short-time Fourier transform as well as its inverse transformation with low-pass filtering. From 0 to T1 is the friction block going, T1 to T2 is the friction block returning, and after T2 is the second reciprocal cycle. In the low-speed stage (2 mm/s), the system may be in a static friction-dominated state with significant periodic stick-slip effects, and although there is relative motion in the macroscopic view, there may be intermittent sticking and sliding, or stick-slip phenomenon, in the microscopic level. However, when the speed is increased to 4 mm/s, the system is in the critical region of the transition from static friction to dynamic friction. The dynamic change of friction (static friction alternating with dynamic friction) is the most intense, resulting in a concentrated release of vibration energy and the formation of amplitude peaks. At higher speeds (6 mm/s to 8 mm/s), dynamic friction dominates, sliding tends to be continuous and stable, and the stick-slip effect is weakened. Therefore, the amplitude fluctuation is most obvious when the speed is 4 mm/s from the point of view of the overall amplitude of acceleration. Figure 8b provides a magnified view of region 1 in Fig. 8a. The positions marked as G1 to G2 in the figure represent the amplitude surges caused by edge impacts at each speed, referred to as "groove­crossing fluctuations" These fluctuations are prominent at 6 mm/s and 8 mm/s, while at 2 mm/s, the impact is minimal due to the slower groove-crossing speed. However, at 4 mm/s, in addition to the groove-crossing fluctuations caused by edge impacts, there are also amplitude surges marked at positions I and II, which occur when passing over the non-textured surface due to vibrations. Although vibrations are present on the non-textured surface at other speeds, the amplitude surges at positions I and II demonstrate that the vibration intensity is significantly lower compared to that at 4 mm/s. 1 . Xn ˆ ˜° ln .., (1) kX . nk. ° ˜. . (2) 4. 2 .° 2 To observe the damping effect of the grooves more obviously, the free decaying oscillating part (.-.) after the "groove-crossing fluctuations" at each speed is targeted. The damping ratios were calculated according to Eqs. (1) and (2) and compared with the non-texture surface, where Xn and Xn+k are the peaks separated by k cycles, and if k = 1, it means the adjacent peaks, and d is the logarithmic decreasing amount. The damping ratios of .-. positions at each velocity are in order .2(.-.)˜0.5674, .4(.-.)˜0.1587, .6(.-.)˜0.5734, .8(.-.)˜0.5884. The damping ratios of the non-woven surface at each speed are .2(n–t) ˜0.1441, .4(n–t) ˜0.1248, .6(n–t) ˜0.1341, and .8(n–t) ˜0.1178.At a speed of 4 mm/s, the system damping after passing through the groove increased by approximately 21.3 % compared to the non-textured surface, while at other speeds, the damping increased by nearly 75 % compared to the non-textured surface, indicating that texture can improve the vibration level of the surface at all speeds. And if we want to further assess the change rule of groove vibration damping ability with velocity, it is inaccurate to use only the change of damping ratio after passing the groove. So, the mean value of the amplitude when the friction block passes through the non-textured surface was calculated, to visualize the effect of velocity on the vibration damping effect of the groove textured. Thus, the vibration acceleration data during the t1– t2 period when the friction block passes over the non-textured surface at each speed were extracted to calculate the average values, which were then compared with those of the smooth surface. The results are shown in Fig 9a. At speeds of 2 mm/s to 8 mm/s, the average values for the non-textured surface are 0.02897 m/s², 0.05132 m/s², 0.03245 m/ s², and 0.03479 m/s², respectively, while the corresponding values for the smooth surface are 0.03484 m/s², 0.06391 m/s², 0.04296 m/s², and 0.04828 m/s². Therefore, as speed increases, the vibration of the friction surface gradually intensifies, which is directly related to the "groove-crossing fluctuations" of friction force observed in the numerical analysis—higher speeds result in larger fluctuations. However, compared to the smooth surface, the average amplitudes of the textured surface are all reduced. The reduction percentages are shown in Fig 9b at speeds of 2 mm/s to 8 mm/s, the reductions are 16.81 %, 19.69 %, 24.46 %, and 27.94 %, respectively. This demonstrates that as speed increases, the vibration damping capability of groove textures gradually improves, which is highly consistent with the enhanced disruption of contact stress concentration observed in the numerical analysis. 3.4 Test Results: Frequency-Domain Analysis of Vibrational Acceleration Through short-time Fourier transform and filtering processing of the time-domain vibration acceleration signals, the frequency-domain diagram shown in Fig. 10 was obtained. Firstly, the Z-direction Fig. 8. Time-domain of vibration acceleration; a) of texture surfaces, and b) at enlargement at “1” Fig. 9. The mean value of vibration acceleration; a) Average amplitude, and b) percentage reduction vibration acceleration of friction surfaces at all speeds exhibited a dominant frequency around 1116 Hz. This dominant frequency is generally associated with the system’s natural frequency or external excitation frequency, and therefore remains stable without changing with speed variations. 2 mm/s, 28.43 dB at 4 mm/s, 24.77 dB at 6 mm/s, and 20.59 dB at 8 mm/s. Overall, the energy suppression effect of groove texture at the dominant frequency enhances with increasing speed. Additionally, all speed conditions exhibited elevated energy levels in low-frequency ranges (113.8 Hz, 125.6 Hz, 138.9 Hz, and 313.8 Hz) induced by the groove texture. These frequency components directly correlate with friction fluctuations and amplified vibration amplitudes resulting from collisions between the friction block and groove edges. As speed increases, the collision frequency accelerates, consequently shifting the energy concentration frequencies associated with groove edge impacts to higher values. 4 CONCLUSIONS To investigate the influence of speed on the improvement of vibration characteristics of friction pair surfaces by groove textures, this study combines numerical analysis and experimental research to analyze the variations in contact stress, friction force, and vibration acceleration on the friction surface at different speeds. The following conclusions are drawn: Firstly, numerical analysis verifies that, at all speeds studied, groove textures can diffuse the concentration of contact stress from the front corner to the rear corner of the friction block, thereby suppressing excessive concentration at the front corner. Simultaneously, the presence of groove edges disrupts and disperses the contact stress on the friction plate surface, reducing its concentration intensity and continuity. Compared to the smooth surface, the maximum contact stress decreases from 1.099·106 Pa to 4.885·105 Pa, a reduction ofapproximately55.56 %. Secondly, as speed increases, the diffusion of contact stress from the front corner to the rear corner of the friction block accelerates. The effectiveness of groove textures in reducing the continuity and intensity of stress concentration on the friction surface gradually improves. By analyzing the variations in friction force at different speeds, the differences in the "disruptive" effect of groove edges are validated. The final test shows that, in the time domain analysis of the surface vibration of the friction sub-surface, the system damping of the freely decaying oscillating part of the groove surface after the “groove-crossing fluctuations”at a speed of 4 mm/s increases by 21.3 % compared with the non-textured surface, and the increase in the rest of the speed is up to 75 %. The mean value of the acceleration amplitude on the grooved surface was lower than that on the non-textured surface at the same velocity, and the decrease was positively correlated with the velocity. Frequency domain analysis shows that the dominant frequency of all surfaces is around 1116 Hz, but due to the effect of “groove-crossing fluctuations”, there is a low-frequency high-energy density region at different speeds. The main frequency energy density decreases with the increase of speed, which confirms that increasing the speed is conducive to improving the vibration characteristics of the groove structure, verifying the numerical analysis of the law, and laying a foundation for the application of synchronous belt tensioning wheel. References [1] Braun, D.,Greiner, C.,Schneider, J., Gumbsch, P. Efficiency of laser surface texturing in the reduction of friction under mixed lubrication. Tribol Int 77 142-147 (2024) DOI:10.1016/j.triboint.2014.04.012. [2] Prem Ananth, M., Poovazhagan, L., Risikesh, J., Tanari, S.V., Sidheswaren, J. Tribological behaviour of titanium alloy due to surface modification technique. Int J Surf Sci Eng 18 130-148 (2024) DOI:10.1504/ijsurfse.2024.138503. [3] Wang, D.W., Mo, J.L., Liu, M.Q., Li, J.X., Ouyang, H., Zhu, M.H. et al. Improving tribological behaviours and noise performance of railway disc brake by grooved surface texturing. Wear 376-377 1586-1600 (2017) DOI:10.1016/j. wear.2017.01.022. [4] Wang, D.W., Mo, J.L., Zhu, Z.Y., Ouyang, H., Zhu, M.H., Zhou, Z.R. How do grooves on friction interface affect tribo-logical and vibration and squeal noise performance. Tribol Int 109 192-205 (2017) DOI:10.1016/j.triboint.2016.12.043. [5] Yang, S., Song, Z., Tong, X., Han, P. Study on the friction noise behaviour of coating textured surface based on STFT time-frequency analysis. Int J Surf Sci Eng 17 295-316 (2023) DOI:10.1504/ijsurfse.2023.135788. [6] Wang, A.Y., Mo, J.L., Qian, H.H., Wu, Y.K., Xiang, Z.Y., Chen, W. et al. The effect of a time-varying contact surface on interfacial tribological behaviour via a surface groove and filler. Wear 478-479 203905 (2021) DOI:10.1016/j. wear.2021.203905. [7] Zhang, D., Zhao, F., Wei, X., Gao, F, Li, P., Dong, G. Effect of texture parameters on the tribological properties of spheroidal graphite cast iron groove-textured surface under sand-containing oil lubrication conditions. Wear 428-429 470-480 (2019) DOI:10.1016/j.wear.2019.04.015. [8] Kumar, B.A., Babu, P.D., Marimuthu, P., Duraiselvam, M. Effect of laser surface texturing on tribological behaviour of grey cast iron. Int J Surf Sci Eng 13 220-235 (2019) DOI:10.1504/ijsurfse.2019.102381. [9] El Fahham, I.M.,Massoud, I.M.. El Gamal, H.A. Performance of non-metallic sliding element textured bearing. Int J Surf Sci Eng 16 1-16 (2022) DOI:10.1504/ ijsurfse.2022.122180. [10] Cao, W., Hu, T., Fan, H., Ding, Q., Hu, L. Fretting behaviour of textured Ti-6Al-4V alloy under oil lubrication. Int J Surf Sci Eng 15 165-164 (2021) DOI:10.1504/ ijsurfse.2021.118199. [11] Wei, Y., Yan, H., Li, S., Wang, X. Numerical and experimental study of a sector-shaped surface texture in friction reduction. Tribol Lett 72 60 (2024) DOI:10.1007/s11249-2024-01863-3. [12] Vladescu, S.-C., Medina, S., Olver, A.V., Pegg, I.G., Reddyhoff, T. Lubricant film thickness and friction force measurements in a laser surface textured reciprocating line contact simulating the piston ring-liner pairing. Tribol Int 98 317-329 (2016) DOI:10.1016/j.triboint.2016.02.026. [13] Aleshin, V.V., Papangelo, A. Friction-induced energy losses in mechanical contacts subject to random vibrations. Int J Solids Struct 190 148-155 (2020) DOI:10.1016/j.ijsolstr.2019.10.026. [14] Janssen, A., Pinedo, B., Igartua, A., Guido Liiskmann, Sexton, L. Study on friction and wear reducing surface micro-structures for a positive displacement pump handling highly abrasive shale oil. Tribol Int 107 1-9 (2017) DOI:10.1016/j. triboint.2016.11.017. [15] Sharma, V.K., Singh, R.C., Chaudhary, R. Improvement in the performance of journal bearings by using lead-free bearing material and surface texturing. Int J Surf Sci Eng 14 238-256 (2020) DOI:10.1504/ijsurfse.2020.110542. [16] Liu, Y., Cui, H., Pan, S., Xia, Y., Zhao, X. Study on the tribological behaviour of surface gradient nanostructure in TB8 titanium alloy. Int J Surf Sci Eng 16 238­257 (2022) DOI:10.1504/ijsurfse.2022.125457. [17] Wang, Z., Li, Y., Bai, F. Wang, C., Zhao, Q. Angle-dependent lubricated tribological properties of stainless steel by femtosecond laser surface texturing. Opt Laser Technol 81 60-66 (2016) DOI:10.1016/j.optlastec.2016.01.034. [18] Zhan, X.,Yi, P., Liu, Y., Xiao, P., Zhu, X., Ma, J. Effects of single- and multi-shape laser-textured surfaces on tribological properties under dry friction. P I Mech Eng C-J Mec 234 1382-1392 2020 DOI:10.1177/0954406219892294. [19] Xu, T., Guan, Q., Ma, C. The Impact of Micro-texture Distribution on the Frictional Performance of Straight Bevel Cylindrical Gears. Stroj vestn-J Mech E 70 582-594 (2024) DOI:10.5545/sv-jme.2024.1033. [20] Uddin, M.S., Liu, Y.W. Design and optimization of a new geometric texture shape for the enhancement of hydrodynamic lubrication performance of parallel slider surfaces. Biosurf Biotrib 2 59-69 (2016) DOI:10.1016/j.bsbt.2016.05.002. [21] Cao, W., Hu, T., Fan, H., Hu, L. Laser surface texturing and tribological behaviour under solid lubrication on titanium and titanium alloy surfaces. Int J Surf Sci Eng 15 50-66 (2021) DOI:10.1504/ijsurfse.2021.114347. [22] Hua, X., Puoza, J.C., Zhang, P., Sun, J. Friction properties and lubrication mechanism of self-lubricating composite solid lubricant on laser textured AISI 52100 surface in sliding contact. Int J Surf Sci Eng 12 228-246 (2018) DOI:10.1504/ijsurfse.2018.094774. [23] Arenas, M.A., Ahuir-Torres, J.I., García, I, Carvajal, H, de Damborenea, J. Tribological behaviour of laser textured Ti6Al4V alloy coated with MoS2 and graphene. Tribol Int 128 240-247 (2018) DOI:10.1016/j.triboint.2018.07.031. [24] Zhang, Y., Sun, H., Li, Q., Sun, K., Mou, Y. and Zhang, S. Research on the cutting performance of self-lubricating tools with microtexture of the front and back surfaces. Stroj vestn-J Mech E 71 136-145 (2025) DOI:10.5545/sv­jme.2024.1181. Acknowledgements This work was supported by the Natural Science Foundation of Jilin Province, the Provincial-Terrestrial Cooperation Fund Project (No.:YDZJ202401325ZYTS) and the 20th Innovative and Entrepreneurial Talent Funding Program of Jilin Province. Received: 2024-11-06, revised: 2025-03-28, 2025-05-11, accepted: 2025­05-28 as Original Scientific Paper. Data Availability The data that support the findings of this study are available from the corresponding author upon reasonable request. Author contribution Wusheng Tang: Project administration, Funding acquisition; Yufei Nie: Validation, Data curation, Writing - review & editing; Zhuo Zhang: Formal analysis, Validation, Writing - original draft, Wei Lin: Conceptualisation, Methodology, Writing - review & editing; Yankai Rong: Data curation, Project administration; Yaochen Shi: Supervision, Methodology, Writing - review & editing; Ning Ding: Methodology, writing - review & editing. All authors made significant contributions to this work andthoroughly reviewed and approved the final version of the manuscript. Numericna in eksperimentalna raziskava vpliva teksture z utori na zmanjšanje vibracij kontaktnih površin Povzetek Za podrobnejše raziskovanje vpliva površinske teksture na vibracijske znacilnosti površin kontaktnega para je bila v tej študiji izdelana tekstura z utori na površini jekla C45 s tehnologijo laserskega oznacevanja, ki je bila uporabljena v kontaktnem paru s kloroprensko gumo. Numericna analiza trenja med površinama pri razlicnih hitrostih je bila izvedena z eksplicitno dinamicno analizo v programskem paketu ABAQUS/Explicit. Rezultati kažejo, da z narašcajoco hitrostjo med kontaktnimi površinami, uporaba teksture z utori zmanjša intenzivnost in zveznost koncentracije kontaktnih napetosti ter hkrati pospešuje prehod kontaktnih napetosti od vstopnega roba proti izstopnemu robu kontakta, kar ucinkovito zmanjša koncentracijo napetosti na vstopnem robu. Eksperimentalni preizkusi trenja med jeklom C45 in kloroprensko gumo pri razlicnih hitrostih so bili izvedeni z izmenicnim gibanjem na preizkusnem stroju za trenje in obrabo tipa HCM­ 5. Rezultati kažejo, da se pri vseh hitrostih, dušenje prostega pojemajocega nihanja na površini z utori po »nihanjih ob preckanju prehodov« bistveno poveca v primerjavi z gladko površino. Z narašcanjem hitrosti se dušilni ucinek teksture z utori na vibracije kontaktnih površin postopno krepi, zmanjševanje gostote energije pri glavni frekvenci vibracij pa postaja vse bolj izrazito. To sovpada z rezultati numericnih analiz in kaže na vpliv hitrosti na zmanjšanje vibracij kontaktnih površin z uporabo teksture z utori. Kljucne besede tekstura z utori, trenje, numericna analiza, hitrost Contents 149 Anton Bergant, Zlatko Rek, Kamil Urbanowicz: Analytical, Numerical 1D and 3D Water Hammer Investigations in a Simple Pipeline Apparatus 157 Qijiang Ma, Zhenbo Liu, Sen Jiang: Study on the Influence of the Splitter Blade Length on Radial and Axial Force of a Centrifugal Pump 169 Fu-sheng Ding, Hong-ming Lyu, Jun Chen, Hao-ran-Cao, Lan-xiang Zhang: Multi-Objective Optimization Design of the Ejector Plate for Rear-Loader Garbage Trucks 179 Jintao Zhang, Zhecheng Jing, Haichao Zhou, Yu Zhang, Guolin Wang: Identification Method of Tire-Road Adhesion Coefficient Based on Tire Physical Model and Strain Signal for Pure Longitudinal Slip 192 Yongping Liu, Qi Chen, Changbin Dong: Tooth Contact Analysis of Involute Beveloid Gear Based on Higher-Order Curve Axial Modification 199 Chao He, Yangzhi Chen, Xiaoxiao Ping, Zhen Chen, Maoxi Zheng, Qinsong Zhang: Geometric Design Method of Lightweight Line Gear Mechanism 207 Wusheng Tang, Yufei Nie, Zhuo Zhang, Wei Lin, YanKai Rong, Yaochen Shi, Ning Ding: Simulation Analysis and Experimental Study on Vibration Reduction Performance of Groove-Textured Friction Pair Surfaces https://www.sv-jme.eu/