Journal of Mechanical Engineering no. 2 year 2022 volume 68 Strojniški vestnik – Journal of Mechanical Engineering (SV-JME) Aim and Scope The international journal publishes original and (mini)review articles covering the concepts of materials science, mechanics, kinematics, thermodynamics, energy and environment, mechatronics and robotics, fluid mechanics, tribology, cybernetics, industrial engineering and structural analysis. The journal follows new trends and progress proven practice in the mechanical engineering and also in the closely related sciences as are electrical, civil and process engineering, medicine, microbiology, ecology, agriculture, transport systems, aviation, and others, thus creating a unique forum for interdisciplinary or multidisciplinary dialogue. The international conferences selected papers are welcome for publishing as a special issue of SV-JME with invited co-editor(s). Editor in Chief Vincenc Butala University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Technical Editor Pika Škraba University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Founding Editor Bojan Kraut University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Editorial Office University of Ljubljana, Faculty of Mechanical Engineering SV-JME, Aškerceva 6, SI-1000 Ljubljana, Slovenia Phone: 386 (0)1 4771 137 Fax: 386 (0)1 2518 567 info@sv-jme.eu, http://www.sv-jme.eu Print: Demat d.o.o., printed in 250 copies Founders and Publishers University of Ljubljana, Faculty of Mechanical Engineering, Slovenia University of Maribor, Faculty of Mechanical Engineering, Slovenia Association of Mechanical Engineers of Slovenia Chamber of Commerce and Industry of Slovenia, Metal Processing Industry Association President of Publishing Council Mihael Sekavcnik University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Vice-President of Publishing Council Bojan Dolšak University of Maribor, Faculty of Mechanical Engineering, Slovenia Cover: Streamlines, cavitation (vapor structures in blue) and cavitation erosion (white spots on gray surface) obtained from the numerical simulation of cavitation erosion on a radial divergent test section. Image Courtesy: Kevorkijan L., Lešnik L. & Biluš I., University of Maribor, Faculty of Mechanical Engineering, Slovenia International Editorial Board Kamil Arslan, Karabuk University, Turkey Hafiz Muhammad Ali, King Fahd U. of Petroleum & Minerals, Saudi Arabia Josep M. Bergada, Politechnical University of Catalonia, Spain Anton Bergant, Litostroj Power, Slovenia Miha Boltežar, University of Ljubljana, Slovenia Filippo Cianetti, University of Perugia, Italy Janez Diaci, University of Ljubljana, Slovenia Anselmo Eduardo Diniz, State University of Campinas, Brazil Igor Emri, University of Ljubljana, Slovenia Imre Felde, Obuda University, Faculty of Informatics, Hungary Imre Horvath, Delft University of Technology, The Netherlands Aleš Hribernik, University of Maribor, Slovenia Soichi Ibaraki, Kyoto University, Department of Micro Eng., Japan Julius Kaplunov, Brunel University, West London, UK Iyas Khader, Fraunhofer Institute for Mechanics of Materials, Germany Jernej Klemenc, University of Ljubljana, Slovenia Milan Kljajin, J.J. Strossmayer University of Osijek, Croatia Peter Krajnik, Chalmers University of Technology, Sweden Janez Kušar, University of Ljubljana, Slovenia Gorazd Lojen, University of Maribor, Slovenia Darko Lovrec, University of Maribor, Slovenia Thomas Lben, University of Bremen, Germany George K. Nikas, KADMOS Engineering, UK Tomaž Pepelnjak, University of Ljubljana, Slovenia Vladimir Popovic, University of Belgrade, Serbia Franci Pušavec, University of Ljubljana, Slovenia Mohammad Reza Safaei, Florida International University, USA Marco Sortino, University of Udine, Italy Branko Vasic, University of Belgrade, Serbia Arkady Voloshin, Lehigh University, Bethlehem, USA General information Strojniški vestnik – Journal of Mechanical Engineering is published in 11 issues per year (July and August is a double issue). Institutional prices include print & online access: institutional subscription price and foreign subscription €100,00 (the price of a single issue is €10,00); general public subscription and student subscription €50,00 (the price of a single issue is €5,00). Prices are exclusive of tax. Delivery is included in the price. The recipient is responsible for paying any import duties or taxes. Legal title passes to the customer on dispatch by our distributor. Single issues from current and recent volumes are available at the current single-issue price. To order the journal, please complete the form on our website. For submissions, subscriptions and all other information please visit: http:// www.sv-jme.eu. You can advertise on the inner and outer side of the back cover of the journal. The authors of the published papers are invited to send photos or pictures with short explanation for cover content. ISSN 0039-2480, ISSN 2536-2948 (online) We would like to thank the reviewers who have taken part in the peer-review process. © 2022 with Authors. The journal is subsidized by Slovenian Research Agency. SV-JME is indexed / abstracted in: SCI-Expanded, Compendex, Inspec, ProQuest-CSA, SCOPUS, TEMA. The list of the remaining bases, in which SV-JME is indexed, is available on the website. Strojniški vestnik - Journal of Mechanical Engineering is available on https://www.sv-jme.eu. Strojniški vestnik - Journal of Mechanical Engineering 68(2022)2 Contents Contents Strojniški vestnik - Journal of Mechanical Engineering volume 68, (2022), number 2 L jubljana, February 2022 ISSN 0039-2480 Published monthly Papers Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš: Cavitation Erosion Modelling on a R adial Divergent Test Section Using R ANS 71 Anz e Sitar: Novel Unsteady State Method of Measuring Permeability Tested on Fabric Materials 82 Leilei Z hao, Yuewei Yu, Jianhu Cao, Weiwei Z hou: Nonlinear Coupled Dynamic Modelling of Driver­ seat-cab System and Biomechanical Behaviour Prediction 90 Marek Woz niak, Adam R ylski, Krz ysz tof Sicz ek: The Measurement of the Wear of Tie R od End Components 101 Guangjian Wang, Dongdong Z hu, Shuaidong Z ou, Yujiang Jiang, X in Tian: Simulation and Experimental R esearch on Electrical Control Anti-backlash Based on a Novel Type of Variable Tooth Thickness Involute Gear Pair 126 Strojniški vestnik - Journal of Mechanical Engineering 68(2022)2, 71-81 Received for review: 2021-08-10 © 2022 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2021-12-17 DOI:10.5545/sv-jme.2021.7364 Original Scientific Paper Accepted for publication: 2022-01-13 Cavitation Erosion Modelling on a R adial D ivergent Test Section U sing R ANS * Luka Kevorkijan – L uka Lešnik – I gnacijo Biluš University of Maribor, Faculty of Mechanical Engineering, Slovenia Cavitation is the phenomenon of fluid evaporation in hydraulic systems, which occurs due to a pressure drop below the value of the vapor pressure. For numerical modelling of this generally undesirable phenomenon, which is often associated with material damage (erosion), there are various mathematical vapor transfer models that have been validated in the past. There are different approaches to predicting cavitation erosion, which have mostly been experimental in the past. Recently various numerical models have been developed with the development of numerical simulations. They describe the phenomenon of cavitation erosion based on different theoretical considerations, such as Pressure wave hypothesis, Microjet hypothesis, or a combination of both. In the present paper, an analysis of the Schnerr-Sauer transport cavitation model was used, upgraded with an erosive potential energy model based on pressure wave hypothesis for cavitation erosion prediction. The extended numerical model has been applied to the case of a radial divergent test section in three different mathematical formulations. The results of simulation were compared and validated to experimental work performed by other authors. The study shows that the distribution of surface accumulated energy agrees with the experimental results, although certain differences exist between formulations. The applied method appears to be appropriate for further use, and to be extended to materials response modelling in the future. Keywords: cavitation, erosion, erosive potential energy, numerical simulation Highlights • The cavitation erosion process can be predicted by cavitation potential energy modelling. • Numerical definition of erosive power is possible with different formulations of cavitation potential. • Cavitation energy focusing accomplishes accurate prediction of cavitation erosion damage. 0 INTRODUCTION Cavitation is the phenomenon that consists of the formation, activity and collapse of vapor structures inside a liquid medium. Usually, the phenomenon is undesired, since it causes many negative effects. Among those negative effects the cavitation erosion is the most complex one, since it combines complicated hydrodynamics with solid material response at the mechanical and metallurgical perspectives. The research of cavitation phenomenon dates to the late 19 th and early 20th centuries. In 1917, Lord R ayleigh [1] derived the equation for calculation of the pressure in a liquid during the collapse of a spherical bubble. In 1949, Plesset [2] upgraded this work by writing the differential equation of the dynamics of a spherical bubble, today known as the R ayleigh-Plesset equation. Mostly experimental research of cavitation damage was conducted during the theoretical research of the cavitation phenomenon.. Hammit [3] studied cavitation damage on the case of fluid flow through a Venturi channel, and introduced the erosion model based on the energy spectrum of a vapor cavity. Vogel et al. [4] observed the dynamics of the collapse of a bubble along a solid wall, and derived the equation for calculation of the potential energy of the collapse. By the end of the 20th century, two hypotheses were known about the origin of cavitation damage. The first hypothesis followed the high value of the pressure in the liquid around the bubble as it was calculated by Lord R ayleigh [1]. The existence of a pressure wave was confirmed by many experimental studies, for example, by Harrison [5]. A second, “ microjet” hypothesis was proposed by Kornfeld and Suvorov [6], while Plesset and Chapman [7] predicted the microjet velocity numerically, and claimed that a microjet can damage the material. Based on the pressure wave hypothesis, Fortes Patella et al. [8] presented an approach to erosion assessment with an energy cascade that transfers the initial potential energy of large vapor structures to the smallest structures, and eventually to the surface walls. A key contribution was the definition of the potential power of cavitation structures based on the total derivative of potential energy. The model upgrade is described in a paper of Leclercq et al. [9], where all of the vapor structures, even those not in direct contact with the wall, were considered. In the paper by Carrat et al. [10] the results of numerical analysis, including the amplitude and frequency of the pressure waves, were presented and validated experimentally. Based on their findings, Schenke and Terwisga [11] proposed a model of forecasting the erosion threat where energy generated by the vapor structures is conserved and transferred to the wall surface via the pressure wave released after the collapse. They also introduced a continuous formulation of energy transfer from wall-distant structures to the wall surface. Li [12] introduced the erosion intensity function, which is based on pressure only. On the other hand, Dular et al. presented the semi-empirical cavitation erosion model based on the microjet hypothesis [13]. The same research group have contributed a lot mainly to the experimental field of cavitation erosion research, [14] and [15]. Based on their experimental findings and semi-empirical model, Peters et al. [16] introduced the multilevel modelling of two-phase flow, where the Euler-Euler and Euler-Lagrange approach were combined. Currently, the research contributions are focused on the use of existing models on practical problem geometries [17], which are, in fact, more complex than in the past, due mainly to the rapid progress of computational capabilities [18]. Brunhart et al. [19] introduced the use of erosion indicators on the case of gaps within a diesel fuel pump. Gomez Santos et al. [20], similarly, presented the comparison of erosion models in engine injectors operating with ethanol. Despite intensive research work and numerous publications in the field of cavitation erosion, the derived models have still not been perfected, and offer many opportunities for upgrading and improvement. The reason for this is mainly the lack of consensus on the predominant mechanism of cavitation erosion. Following this, Arabnejad et al. [21] joined both erosion mechanisms. In their work, the collapse induced kinetic energy is divided into the kinetic energy of the pressure wave and kinetic energy of the microjet. In the present work, a model by Schenke and Terwisga [11] has been applied to a radial divergent test section, which is a common test case for cavitation studies found in the literature. The mesh was generated with a meshing module within ANSYS Fluent 2020 R 2, cavitating flow and cavitation erosion were calculated using the computational fluid dynamics (CFD) software ANSYS Fluent 2020 R 2, where the model to predict cavitation erosion by Schenke and Terwisga [11] was implemented via a user defined function (UDF) which was written in the C programming language, compiled and linked to the solver. Different authors analyz ed cavitation on this flow case numerically, but they applied different cavitation erosion models [16]. In this work, we applied an erosive potential energy model based on the pressure wave hypothesis with three different mathematical formulations. The transfer of energy from the cavitation structures to the surface was accomplished with and without the energy focusing method. The numerical results show good trends and agreement with the experimental results, with certain limitations, which could be a subject of further studies and model optimiz ations. 1 CAVITATION EROSION 1.1 Homogeneous Mixture Method In the following, we focus on the general formulation, and treat the fluid as a continuous mixture consisting of the volume fraction of the vapor phase and the volume fraction of the liquid phase (1 – a). The flow is considered as isothermal, pure phases are treated as incompressible. The mixing rule applies to the properties of the mixture, the mixture density (.) and the mixture viscosity (µ) ˜°.˜v ..°˜l , 1 ˆ.(1 ) ˜°˜v ..°˜l , (2) 1 .ˆ. where the liquid and vapor densities are .= l 998.85 kg/m3 and .v = 0.01389 kg/m3 . The liquid and vapor viscosities are µ= 0.001 1 Pa· s and µv = 9.63 · l –6 10 Pa· s. 1.2 Conservation of Vapor Phase Mass The vapor mass is defined as mv = .vV v = .v aV , where V is the vapor volume and V is the mixture volume. v According to this, the mass flow of vapor is Dm D ..V vv m .v ˜dV ˜dV , (3) °° VDt VDt and the source of the vapor fraction is . .°.  v ...ˆu .dV S dV , (4) ˜V v ˜V v v °t  . where u is mixture velocity. If we write Eq. (4 ) in differential form and consider pure phases incompressible, we get . ˜.. °....u ˆ.S .v , (5) ˜t where S av presents the source of the vapor phase. 1.3 Conservation of the Mixture Mass Similarly to Eq. (5) , the mass conservation for the mixture is written as . ˜.. °....u ˆ.0. (6) ˜t °° From ˜°.˜..°˜l  .l 1 we v ˆ.and ˜. °v .°l can write ˜.. .. .. °....ˆlu (7) u ˆ.u .., ˜t .ˆ. vl and after rearrangement ˜... ... °...ˆ..u , (8) u .˜tv ˆl .. which presents the conservation equation of the mixture mass written by the vapor volume fraction a. 1.4 Conservation of the Mixture Momentum Within a homogeneous mixture, the approach for the momentum conservation is written for the mixture as ˜°u . .... .. ..ˆ°uu ......ˆ, (9) p ˜t where t is the viscous stress tensor and is evaluated as .. .. T ˜.°ˆ...uu.. (10) a1 is the model constant equaling 0.31. Correction of turbulent viscosity is achieved by modifying the mixture density function f (.) such that a lower turbulent viscosity is achieved in the presence of both the liquid and vapor phases ˜˜v .n °ˆ f °.˜.˜., (12) vn ˆ1 ˜˜ °l ˆv .where n is an arbitrary exponent with a recommended value of 10. This modification of turbulent viscosity was implemented in ANSYS Fluent via a UDF. 1.6 Cavitation Modelling For the cavitation modelling, the additional conservation equation of vapor mass is solved. . v . ˜.°...u .ˆ.S. (13) iv . v ˜t The Schnerr-Sauer cavitation model was used, which states that mass transfer source term Sav is  pp ˆ °°32 .v . vl ˜v .1 .˜v ˆif pp .v °Rb 3 °l S . , (14) ˜ v °°vl ˜v .1 .˜v ˆ32 . pp .v ˆif pp . °R 3 °v bl From Eq. (9) it is noticeable that body forces where the bubble radius is (gravity) are omitted in this work. The flow was 13 ../ considered as incompressible, while compressible . v behavior of the mixture is exhibited during the phase Rb ˜.. (15) 4 transition, as in the work of [21] and [22]. ..n 0 .1 °v . ˆ3  1.5 Turbulence Modelling The bubble number density was specified with a value 1· 10 11 . The vapor pressure was set to pv = 1854 To model turbulent flow, the R ANS approach was adopted, with the SST k – . model. In order to address the transient cavitation effects that arise due to the compressibility of a liquid-vapor mixture, the turbulent viscosity µt is modified according to R eboud et al. [23], and Coutier-Delgosha et al. [24] as follows f °1 ˆ.k ˜t ., (1 1) ..1 SF 2  max ,  .* a 1 .  where k is the turbulence kinetic energy, . is the specific rate of dissipation of the turbulence kinetic energy, a* is a damping coefficient, S is the strain rate magnitude, F 2 is the second blending function and Pa. 1.7 Cavitation Erosion Potential Hammit [3] and Vogel et al. [4] described the potential energy of a cavitation bubble in terms of work done as a consequence of the pressure difference between the vapor inside the bubble and the surrounding liquid. If we analyz e the spherical cavitation bubble with volume V , then the potential energy of the cavitating V bubble E pot equals E ˜p °pV , (16) pot .dv .V where p is the pressure in the surrounding liquid and dpv is the vapor pressure inside the bubble. Authors [11] highlighted the problem of defining p for the d purposes of numerical simulations, where the pressure in the surrounding liquid is evaluated locally in each computational grid cell. In that case pequals pv and d the calculated potential energy equals z ero. In this work p is taken as a time-averaged value in each d computational grid cell as proposed in [11] pd ˜* tpt dt , (17) 1 .* °. 0 t where p(t) is an instantaneous pressure that is then time-averaged over the current simulated time t* . The cavitation potential power P pot is introduced as the total derivative of the potential energy E pot Dpd °pv V ..DV P ˜V ..p °p .. (18) pot Vdv Dt Dt Since pv = const. (isothermal flow) and . pd ˜f xt , const. the potential power equals °.. Dp DV Ppot ˜d VV °V .pd .pv .. (19) Dt Dt The first term in Eq. (19) is neglected, as it has been found that it is at least an order of magnitude lower than the second term [22]. Finally, the cavitation potential power is written as DVV P ˜.p °p .. (20) pot dv Dt To calculate the cavitation potential power numerically, it is appropriate to calculate the instantaneous change of the volume specific potential energy [11] or cavitation erosion potential e .pot = Ppot , (21) Vcell where V presents the volume of the computational c ell grid cell. A positive value of e.pot is then a consequence of an increase of the vapor cavity volume, and a negative value is a consequence of the reduction in cavity volume. It is possible to rearrange Eq. (20) and rewrite it in modified forms, dependent on different variables. In the present study we used three equation forms, namely, a form dependent on the vapor volume fraction epot,a, a form dependent on the velocity divergence epot. and a form dependent on the vapor ,divu phase source term epotS v . , a 1.7.1 Cavitation Erosion Potential . epot,a If we calculate the total derivative in Eq. (16) , consider the relations V . V and V = const. and combine it with Eq. (21) we get c ellc ell ° ..˜° e ˜pot , ˜°..ˆ˜pd pv . u (22) ..t  Eq. (22) defines the cavitation erosion potential with vapor volume fraction a. It describes the change of potential due to the bubble volume decrease during collapse. 1.7.2 Cavitation Erosion Potential e˜pot divu° , Combining the mass conservation Eq. (8) and °° .l ˜. with Eq. (21) we get the cavitation °v .°lerosion potential e˜pot divu° , ° e ˜°˜.p °p ...ˆu °. (23 ) pot divu , dv .l °.v Eq. (23) presents the second form of cavitation erosion potential, which is dependent on velocity divergence. 1.7.3 Cavitation Erosion Potential . epotS v , a Combining the conservation Eqs. (5) and (8) we get ˜ˆ.. . °...S , (24 ) u ˆ .ˆ ˜t lv and a third form of equation for erosion potential ° e ., ..p .p ˆS . (25 ) pot S ˜v dv ˜v ° l In Eq. (25) the volume change is now expressed via the mass transfer source term Sav . 1.8 Collapsing Bubble Energy Balance If it is assumed that, as the vapor cavity is collapsing, reduction of it’ s potential energy is converted instantaneously into the radiated power, then evaluation of Eqs. (22), (23) and (25) gives this instantaneous radiated power e . rad ˜°e . pot , (26 ) where a negative sign is needed to obtain a positive value for instantaneous radiated power, as only a reduction in cavitation erosion potential e.pot <0 represents vapor cavity collapse, and only the collapse stage is considered in the model. To predict the cavitation erosion on a solid surface better, the potential energy conservation hypothesis was introduced by Schenke et al. [25] and Mellisaris et al. [22]. It states that during the vapor cavity collapse potential energy is first converted into the kinetic energy of the vapor cavity interface, and is accumulating until the last moment of collapse, when it is finally released. This delayed release of energy from the collapsing cavity is named energy focusing, as the release of energy is more focused towards the center of the collapse compared to the case where energy is released instantaneously as described by Eq. is described by the additional conservation equation for collapse induced kinetic energy e ˜ ˜˜ °...ˆ..°.ˆt , ue (27) i rad ˜t . where ui presents the collapse induced velocity, and e .rad ˜°t presents the radiated power of the pressure wave released in the bubble center at the final stage of collapse. The discretiz ed explicit model for transport Eq. (27) could be written as: .˜t ..| t . ˜| .t ˜t .t t .ˆ°K 1 e . t ˆ˜P ˆ1  . (28) t .ˆ..u  .1 ..pot t P in Eq. (28) presents the projection operator, u which ensures that the specific kinetic energy of the collapse accumulates at the inner side of the bubble interface, K presents the conservation parameter that ensures global energy conservation and ß is the criterion which determines whether the induced energy is released at a given moment. at the exit from the domain was used for its value. The criterion av = 0 represents the end of the collapse when there is no more vapor phase in the computational cell. The energy released after the collapse (at the end of time step t) is expressed as the specific power of the pressure wave .ˆ..| t °. (32) e .rad t ˜.t .t An infinite wave propagation speed is assumed in this model. The energy received by a surface element S is then expressed in terms of the instantaneous surface specific impact power °°° cell ˆ .  .. ° ...   rad 4 . where xcell is the position vector of the computational . cell center, xS is the position vector of the wall surface . element center, and n is the surface element normal. By integrating Eq. (33) over the simulated time (t), we obtain instantaneous surface specific impact energy with units [ J/m2] °t 0 eS ˜e . Sdt . (34) The described model to predict cavitation was implemented inside ANSYS Fluent via a UDF, which was written in the C programming language, compiled, and linked to the solver. Eqs. (22), (23) , (25) and (28) to (34) were evaluated at the end of every time step. Energy accumulated on a given face of the surface is then stored as a sum of the previously stored value and newly computed value by Eq. (3 3) . 2 NUMERICAL TEST CASE AND SETUP 2.1 Geometry 3 V ° x cell x S x x n (26) . This accumulation and delayed release of energy e ˜S ˜ 1  e ˜ S dV , (33) °° . The conventional radial divergent test section was .. °. . . chosen for the test and validation of the cavitation erosion models [26]. The radial divergent test section u max ˆˆ. . , (29) Pu ˜ 0 . u ,  . geometry is shown in Fig. 1. . °PdV u V .t 2.2 Computational Grid and Boundary Conditions K ˜ , (30) edV °. pot V A computational grid of a 45° sector was generated in ANSYS Fluent Meshing using ANSYS Mosaic °.0 . pp .1, if and . (31) ˜. meshing technology, which resulted in Poly-Hexcore v ˆ . . 0,else. grid topology. Additionally, a prismatic inflation layer p8 in Eq. (29) presents the external pressure was generated on the walls. The final mesh consisted of infinitely far from the bubble. This pressure represents 1,396,703 cells with a minimum orthogonal quality of the threshold for the propagation of the pressure wave 0.28 and maximum aspect ratio of 65. The maximum ++ outwards from the point of collapse, and the pressure wall y value was 51.9 and the minimum wall y value + was 0.174. Special care was taken to achieve y < 5 around the radius before the radial divergent test + section, and to achieve 30 < y< 300 on the walls of the extended inlet and outlet parts of the computational domain. The computational domain was extended downstream with the outlet boundary positioned 100 mm in the radial direction from the radial divergent test section axis to improve numerical convergence. The mesh was refined around the radial divergent test section radius for improved resolution of cavitation dynamics, and, therefore, cavitation erosion prediction. Fig. 1. Radial divergent test section The results of the grid independence study are presented in Table 1. A medium grid was chosen as it provided us with an acceptable computation time and numerical uncertainty, as given by the grid convergence index (GCI). All simulations were performed on a personal computer using 10 CPU cores of the Intel Core i9-10940X stock CPU with 64 GB of R AM. The time needed to obtain convergence on the medium mesh was approximately 10 days. For GCI calculations we chose the force exerted on the bottom plate, the fluid velocity at the outlet and turbulence kinetic energy at the outlet. In all three cases the estimated numerical error was less than 5 % for both the medium and the fine grids [27]. Table 1. Computational grid independence study Grid Number of cells Force on bottom wall [N] Outlet velocity [m/s] Outlet turbulence kinetic energy [J/kg] Coarse 475271 3802 3.64 0.196 Medium 1396703 3798 3.57 0.195 Fine 5175737 3799 3.51 0.187 GCImedium [%] - 0.02 4.42 0.51 GCIfine [%] - 0.001 2.59 1.99 Kevorkijan, L. – Lešnik, L. – Biluš, I. Fig. 2. Computational grid; a) isometric view of the whole grid with shown boundary conditions, b) detailed isometric view of the refined grid, and c) detailed front view of the refined grid 2.3 Numerical Setup Solution of the Navier-Stokes equations was achieved using the SIMPLE algorithm, gradients were evaluated using the least squares cell-based method, interpolation of pressure cell face values was achieved by the PR ESTO! interpolation scheme and the QUICK spatial discretiz ation scheme was adopted for all equations. Temporal discretiz ation of equations was achieved by the Bounded Second-Order Implicit time integration scheme. A time step of 1· 10 –6 s was chosen, which resulted in the maximum Courant number C = 1.59. The iterative solution of equations within each time step was limited to 100 iterations, however, a scaled residuals convergence criterion of 1· 10 –6 was achieved before this limitation. Overall, 0.012249 s of physical time were simulated. 3 RESULTS The extent of cavitation predicted by the Schnerr-Sauer cavitation transport model at different times of numerical simulation is shown in Fig. 3. The black dotted line shows the cavitation closure line positioned 32 mm in a radial direction from the radial divergent test section axis, which was determined experimentally [26]. From the cavitation extent comparison, it is visible that the predicted cavitation sheet closure coincides with the experimental position, as cavitation number s = 0.9 was achieved. According to this, it can be concluded that the used cavitation model, boundary conditions and physical properties were determined properly to predict the cavitation dynamics correctly. Qualitative assessment of cavitation erosion by the numerical model is presented in Fig. 4, with three approaches for calculating cavitation erosion potential compared to the picture of the sample from the experiment by [26]. Both energy focusing, where radiated power was obtained from Eq. (32) , and non-focusing approaches, where radiated power was obtained from Eq. (26 ), are Fig. 3. Transient sheet cavitation displayed as an iso-surface of a= 0.9 at different times; a) top view; b) three- dimensional display of cavitation structures viewed from the side Fig. 4. Comparison between the cavitation erosion models and experiment [26]; a) accumulated energy on the surface without energy focusing; b) accumulated energy on the surface with energy focusing presented. Energy focusing resulted in a sharper, more pitted image compared to a more blended picture without energy focusing. Since in the case of energy focusing all previously accumulated collapse induced kinetic energy is released only at the final stage of collapse, it is reasonable to expect accumulation of radiated energy on the surface in a smaller area, and higher magnitude compared to the non-focusing case. In both cases, with energy focusing and without energy focusing, it is apparent that modelling cavitation erosion potential based on partial derivative of vapor volume fraction (a t) and source term in vapor transport equation (S av) produced a similar erosion pattern and location of accumulated energy on the bottom plate surface, which are comparable with the experiment (shown in the background of Fig. 4) . However, calculation based on the divergence .. of velocity .˜°u. did not predict a similar erosion pattern and location of accumulation of energy on the bottom plate. The comparison of the eroded z one position was similar for both (focusing and non-focusing) cases. The deviation of the position of the damaged area between the results of the numerical simulation and the experiment visible in Fig. 4 indicates a strong sensitivity of the numerical erosion prediction model on correct prediction of the cavitation dynamics. Fig. 5. Total vapor volume time evolution in the domain Fig. 6. Pressure time evolution, 20 mm from the radial divergent test section axis in a radial direction The total vapor volume time evolution in the numerical domain is shown in Fig. 5 for the time interval 0.002 s to 0.012 s. It is apparent that the total volume of vapor in the domain changes periodically, which confirms the transient nature of the extent of cavitation in the radial divergent test section. Additionally, Fig. 6 displays the local transient behavior of the vapor cloud shedding, with minimums in pressure between the peaks, indicating the presence of cavitation. It was observed that vapor cloud shedding occurs with higher frequency than the global (total) change in vapor volume in the domain. In a half of a total vapor volume change period (the dashed line frame in Figs. 5 and 6) , at least three transitions occur between the liquid and vapor phases. Within the simulated time of 0.012249 s multiple vapor cloud shedding events occurred, which allowed for a consistent analysis of cavitation erosion. Fig 7. shows the time evolution of the total cavitation potential power in the domain within this time period, calculated with a formulation based on the partial derivative of vapor volume fraction . Fig. 7. Total cavitation potential power in the domain Within the observed time window (the dash lined frame in Figs. 5 and 6) release of accumulated kinetic energy for the energy focusing case is presented in Fig. 8. Comparison between Figs. 7 and 8a further demonstrates and explains the difference in extent and magnitude of the predicted accumulation of energy between the energy focusing case and the non-focusing case. In the non-focusing case cavitation potential power was released instantaneously, so Fig. 7. also represents the total radiated power of pressure waves in the domain, for the non-focusing case. It is again apparent that energy focusing results in fewer events with higher magnitude of pressure wave power radiation (at the end of vapor collapse). Because of Fig. 8. a) Total radiated power of pressure waves generated by cavitation implosions in the domain; b) total vapor collapse induced kinetic energy of vapor in the domain the energy conservation by the model, kinetic energy reduces at the end of vapor collapse when wave energy is released. Two pronounced events of release of energy are visible in Fig. 8. at around 0.01 15 s, first a smaller peak occurred, followed by a larger peak (Fig. 8a ). Hence, there are two pronounced corresponding drops in kinetic energy (Fig. 8b) . 4 DISCUSSION The presented study results show that the used mathematical-physical model with different formulations for prediction of vapor phase potential power, and, thus, cavitation erosion, predicts the position of damaged surfaces on the geometry of the radial divergent test section satisfactorily, which confirms the findings of other researchers on different flow cases [22] and [25] and different cavitation erosion models [16]. A more detailed comparison of the positions of damaged surfaces indicates discrepancies between the individual formulations. Given the fact that a verified cavitation transfer model was used, it can be concluded that the largest deviation of the damage positions from the experimentally observed ones occurred in the case of using a formulation with velocity divergence, both in the case of focusing and without focusing. The phenomenon was a consequence of numerical error, and is consistent with the observations of researchers who analyz ed the flow over an isolated wing and marine propeller [22] and [25]. The formulations with the source term and the derivative of the vapor phase fraction confirm the importance of a correct description of the dynamics of cavitation structures. There was practically no difference in the calculated positions of damage between the applied formulations with energy focusing and without energy focusing. Conserving the potential energy of cavitation by focusing indicates an appropriate step towards a better agreement of the shape of the damaged surface with the experimental observations, where the damage was more concentrated and the eroded surface was characteristically pitted, while, in the case without focusing, the damage was more continuous and did not reflect the experimentally observed damage of eroded surfaces [26]. It has to be noted however, that to directly link the material pitting prediction to the cavitating flow simulation a change of scales would be necessary. Furthermore, material work hardening would have to be taken into account via proposed material strengthening models [26] and [28]. 5 CONCLUSIONS The performed analysis confirms that the formulation of cavitation erosion potential with the source term is comparable to the formulation with the derivative of vapor phase volume fraction in the case under consideration. An appropriate description of cavitation dynamics is crucial. The Schnerr-Sauer cavitation model predicts the extent and dynamics of cavitation correctly. In relation to that, it would be sensible to use cavitation transfer models of other authors, such as Singhal et al. [29] and Kunz et al. [30]. However, from the thermodynamic point of view this approach to model cavitation is not as appropriate, as discussed in [31]. The applied model transfers potential energy from cavitation structures to the surface, is energy conservative, which makes it a good basis for predicting the response of various materials exposed to the phenomenon of cavitation. In modelling their response to cavitation, the model used in this work could be upgraded by taking into account the mechanical properties of the material, which would allow the calculation of the actual depth of damage on different surfaces. This would allow a quantitative evaluation of the model and it’ s different formulations used in this work. 6 ACKNOWLEDGEMENTS The authors wish to thank the Slovenian R esearch Agency (AR R S) for the financial support in the framework of the R esearch Programme P2-0196 R esearch in Power, Process and Environmental Engineering. 7 REFERENCES [1] Lord Rayleigh, O.M.F.R.S. (1917). VIII. On the pressure developed in a liquid during the collapse of a spherical cavity. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, vol. 34, no. 200, p. 94-98, DOI:10.1080/14786440808635681. [2] Plesset, M.S. (1949). The dynamics of cavitation bubbles. Journal of Applied Mechanics, vol. 16, no. 3, p. 277-282, DOI:10.1115/1.4009975. [3] Hammitt, F.G. (1963). Observations on cavitation damage in a flowing system. Journal of Basic Engineering, vol. 85, no. 3, p. 347-356, DOI:10.1115/1.3656601. [4] Vogel, A., Lauterborn, W., Timm, R. (1989). Optical and acoustic investigations of the dynamics of laser-produced cavitation bubbles near a solid boundary. Journal of Fluid Mechanics, vol. 206, p. 299-338, DOI:10.1017/S0022112089002314. [5] Harrison, M. (1952). An experimental study of single bubble cavitation noise. The Journal of the Acoustical Society of America, vol. 24, p. 454-454, DOI:10.1121/1.1917513. [6] Kornfeld, M., Suvorov, L. (1944). On the destructive action of cavitation. Journal of Applied Physics, vol. 15, no. 6, art. ID 495, DOI:10.1063/1.1707461. [7] Plesset, M.S., Chapman, R.B. (1971). Collapse of an initially spherical vapour cavity in the neighbourhood of a solid boundary. Journal of Fluid Mechanics, vol. 47, no. 2, p. 283­290, DOI:10.1017/S0022112071001058. [8] Fortes Patella, R., Reboud, J.-L., Briancon-Marjollet, L. (2004). A Phenomenological and numerical model for scaling the flow agressiveness in cavitation erosion. Numerical Prediction of Cavitation Erosion, Delft. [9] Leclercq, C., Archer, A., Fortes-Patella, R., Cerru, F. (2017). Numerical cavitation intensity on a hydrofoil for 3D homogeneous unsteady viscous flows. International Journal of Fluid Machinery and Systems, vol. 10, no. 3, p. 254-263, DOI:10.5293/IJFMS.2017.10.3.254. [10] Carrat, J.-B., Fortes-Patella, R., Franc, J.-P. (2019). Experimental and numerical investigation of the erosive potential of a leading edge cavity. International Journal of Fluid Machinery and Systems, vol. 12, no. 2, p. 136-146, DOI:10.5293/IJFMS.2019.12.2.136. [11] Schenke, S. van Terwisga, T.J.C. (2019). An energy conservative method to predict the erosive aggressiveness of collapsing cavitating structures and cavitating flows from numerical simulations. International Journal of Multiphase Flow, vol. 111, p. 200-218, DOI:10.1016/j.ijmultiphaseflow.2018.11.016. [12] Li, Z. (2012). Assessment of Cavitation Erosion with a Multiphase Reynolds-averaged Navier-Stokes Method, PhD Thesis, TU Delft, Delft. [13] Dular, M., Stoffel, B., Širok, B. (2006). Development of a cavitation erosion model. Wear, vol. 261, no. 5-6, p. 642-655, DOI:10.1016/j.wear.2006.01.020. [14] Dular, M. Petkovšek, M. (2015). On the mechanisms of cavitation erosion - Coupling high speed videos to damage patterns. Experimental Thermal and Fluid Science, vol. 68, p. 359-370, DOI:10.1016/j.expthermflusci.2015.06.001. [15] Jian, W., Petkovšek, M., Houlin, L., Širok, B., Dular, M. (2015). Combined numerical and experimental investigation of the cavitation erosion process. Journal of Fluids Engineering, vol. 137, no. 5, art. ID 051302, DOI:10.1115/1.4029533. [16] Peters, A. el Moctar, O. (2020). Numerical assessment of cavitation-induced erosion using a multi-scale Euler-Lagrange method. Journal of Fluid Mechanics, vol. 894, DOI:10.1017/ jfm.2020.273. [17] Periasamy, G., Mouleeswaran, S., Venugopal, P., Perumal, C. (2021). Investigation of hydrodynamic flow characteristics in helical coils with ovality and wrinkles. Strojniški vestnik - Journal of Mechanical Engineering, vol. 67, no. 11, p. 570­579, DOI:10.5545/sv-jme.2021.7374. [18] Dai, Y., Zhang, Y.Y., Bian, J.N., Han, K., Zhu, X., Huang, Z.H., Xie, Y. (2021). CFD Simulation on hydrodynamics of underwater vehicle with ducted propellers. International Journal of Simulation Modelling, vol. 20, no. 3, p. 595-605, DOI:10.2507/IJSIMM20-3-CO14. [19] Brunhart, M., Soteriou, C., Daveau, C., Gavaises, M., Koukouvinis, P., Winterbourn, M. (2020). Cavitation erosion risk indicators for a thin gap within a diesel fuel pump. Wear, vol. 442-443, art. ID 20324, DOI:10.1016/j. wear.2019.203024. [20] Gomez Santos, E., Shi, J., Venkatasubramanian, R., Hoffmann, G., Gavaises, M., Bauer, W. (2021). Modelling and prediction of cavitation erosion in GDi injectors operated with E100 fuel. Fuel, vol. 289, art. ID 119923, DOI:10.1016/j. fuel.2020.119923. [21] Arabnejad, M.H., Svennberg, U., Bensow, R.E. (2021). Numerical assessment of cavitation erosion risk using incompressible simulation of cavitating flows. Wear, vol. 464­465, art. ID 203529, DOI:10.1016/j.wear.2020.203529. [22] Melissaris, T., Schenke, S., Bulten, N., van Terwisga, T.J.C. (2020). On the accuracy of predicting cavitation impact loads on marine propellers. Wear, vol. 456-457, art. ID 203393, DOI:10.1016/j.wear.2020.203393. [23] Reboud, J.-L., Coutier-Delgosha, O., Benoit Stutz, B. (1998). Two-phase flow structure of cavitation: Experiment and modeling of unsteady effects. Third International Symposium on Cavitation, p. 1-7. [24] Coutier-Delgosha, O., Fortes-Patella, R., Reboud, J.L. (2003). Evaluation of the turbulence model influence on the numerical simulations of unsteady cavitation. Journal of Fluids Engineering, vol. 125, no. 1, p. 38-45, DOI:10.1115/1.1524584. [25] Schenke, S., Melissaris, T., van Terwisga, T.J.C. (2019). On the relevance of kinematics for cavitation implosion loads. Physics of Fluids, vol. 31, no. 5,, art. ID 052102, DOI:10.1063/1.5092711. [26] Franc, J.-P. (2009). Incubation time and cavitation erosion rate of work-hardening materials. Journal of Fluids Engineering, vol. 131, no. 2, art. ID 021303, DOI:10.1115/1.3063646. [27] Roache, P.J. (1994). Perspective: A method for uniform reporting of grid refinement studies. Journal of Fluids Engineering, vol. 116, no. 3, p. 405-413, DOI:10.1115/1.2910291. [28] Liu, L., Guo, H., Yu, P. (2021). A model for material strengthening under the combined effect of cavitation-bubble collapse and Al2O3 particles, and its test verification. Strojniški vestnik - Journal of Mechanical Engineering, vol. 67, no. 1-2, p. 36-44, DOI:10.5545/sv-jme.2020.6937. [29] Singhal, A.K., Athavale, M.M., Li, H., Jiang, Y. (2002) Mathematical basis and validation of the full cavitation model. Journal of Fluids Engineering, vol. 124, no. 3, p. 617-624, DOI:10.1115/1.1486223. [30] Kunz, R.F., Boger, D.A., Stinebring, D.R., Chyczewski, T.S., Lindau, J.W., Gibeling, H.J., Venkateswaranm S., Govindan, T.R. (2000). A preconditioned Navier-Stokes method for two-phase flows with application to cavitation prediction. Computers & Fluids, vol. 29, no. 8, p. 849-875, DOI:10.1016/ S0045-7930(99)00039-0. [31] Goncalvès, E. Patella, R.F. (2011). Constraints on equation of state for cavitating flows with thermodynamic effects. Applied Mathematics and Computation, vol. 217, no. 11, p. 5095­5102, DOI:10.1016/j.amc.2010.07.056. Strojniški vestnik - Journal of Mechanical Engineering 68(2022)2, 82-89 Received for review: 2021-10-26 © 2022 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2021-12-29 DOI:10.5545/sv-jme.2021.7453 Original Scientific Paper Accepted for publication: 2022-01-11 Novel U nsteady S tate Method of Measuring Permeability Tested on Fabric Materials * Anz e Sitar Fondaz ione Bruno Kessler, Italy The necessity of measuring permeability is encountered in number of applications and some of them have adequate measurement techniques available. However, in the fields of freeze drying, composite scaffolds, 3D printed materials, etc., the sample materials have unique shapes and sizes, which disables the use of standardized equipment. Hence, a novel method was developed to determine the permeability of a sample material by analyzing a high frequency unsteady state pressure measurement acquired during the permeation of a working fluid. The method is also suitable for in-line measurements due to the rapid acquisition and analysis. The developed method is novel, hence a comparison of the acquired results and the permeabilities measured with a standard device designed for measuring air permeability of fabric materials was made. Air was the permeating fluid in the novel unsteady state as well as in the referential steady state measurements. The measured permeabilities ranged from approximately 8 Da to 50 Da for the analyzed five fabric materials. The comparative analysis yielded an encouraging linear fit with a R2 = 0.98 by including only one unsteady state measurement of each fabric sample. In addition, the presented method was capable of detecting the difference in permeability of the freeze-dried 5 wt% and 12 wt% aqueous mannitol solutions, which exhibit different permeabilities due to the different porosities after the same process of lyophilization. The possibility of determining permeabilities of freeze-dried pharmaceuticals is essential for lyophilization process optimization. Keywords: novel measurement method; permeability; flow resistivity; unsteady state measurements; freeze-dried pharmaceuticals Highlights • Novel unsteady state method of determining permeability is proposed. • Inexpensive, reliable and rapid acquisition of the experimental results is presented. • The method is compared to standardized steady state measurements. • Sample materials’ size, shape and permeability is arbitrary. • The first solution for determining permeability of freeze-dried pharmaceuticals is demonstrated. 0 INTRODUCTION The technics of permeability determination are commonly divided into two categories: steady state and unsteady state methods. Both are most thoroughly studied in the field of geology due to the importance of rock permeability in the oil and gas industry. The two measurement categories are further divided into subcategories, which have unique features, recently summariz ed in Gao and Li [1]. The unsteady state methods have different variants with a common denominator. Namely, the measured signal is transient and the analysis of the experimental results requires a thorough theoretical knowledge, which was initially presented by Brace et al. [2] and afterwards modified by several authors [3] to [5]. Compared to the steady state experiments, the unsteady state measurements are more rapid, as shown by Z hong et al. [6], and more frequently applied for determining lower permeabilities, whereas the analysis of the results is generally much more complex and requires a profound theoretical background [7] as well as additional parameters: material and fluid compressibility, porosity, viscosity, etc. [8] and [9], which are often difficult to acquire. Yong et al. [10] presented large discrepancies in liquid cross-plane permeability of glass fiber fabrics measured in 26 laboratories. However, the filter, membrane and fabric materials present no extraordinary challenges regarding the gas permeability measurements especially when determining the cross-plane steady state permeability as seen in [11]. Sheet materials are therefore ideal for verification of novel measurement methods, as there are several widely used and also standardiz ed measurement methods available for determining the permeability. Most commonly used standardiz ed equipment and procedures for measuring permeability are given in ASTM standards [12] and [13]. The currently available state of the art in steady or unsteady state methods of determining permeability is unsuitable to analyz e non-extractable and arbitrarily shaped materials, which is the main advantage of the proposed technique. Hence, the use of this method is foreseen in the field of freez e-drying pharmaceuticals, as the product is not suitable for *Corr. Author’s Address: Fondazione Bruno Kessler, Via Sommarive 18, 38123 Trento, Italy, asitar@fbk.eu extraction and the permeability presents one of the most important parameters for designing the drying process. Currently used approaches of determining the porosity or permeability rely on Scanning Electron Microscopy [14] to [16]; X -ray microtomography [16] to [18] coupled with numerical simulations [19] and [20]; or solely on computational fluid dynamics [21] to [23]. Direct measurements of permeability in freez e dried products are non-existent, as the materials are brittle and therefore difficult to manipulate or extract from their containers, however the need to determine permeability is explicitly given in Fissore et al. [24] and Bobba et al. [25], as it would enable a more precise modelling and optimiz ation of lyophiliz ation cycles. The permeability is the most influential material property in the scope of lyophiliz ation process, as the water vapor is required to permeate via sublimation through the material during the freez e drying process. The sublimation flow during drying is hindered by the materials flow resistance, therefore it is crucial to measure the product’ s permeability in order to evaluate and potentially improve the protocol of freez e drying. The novel method allows for measuring the fluid permeability rate of sheet materials as well as other arbitrarily shaped materials. The method is flexible to perform measurements with different working fluids and various materials with a wide range of permeabilities. Additionally, the boundary and initial conditions of the experimental measurement can be adjusted to the expected process parameters to determine the permeability, which will affect the process of interest. Hence, the temperatures, pressures and volumetric flows during the experimental measurements can imitate the boundary conditions during freez e drying and therefore the intrinsic permeability is measured. To avoid elaborate physical modeling and acquisition of additional thermodynamic properties, a comparative analysis of the experimental results is recommended, which also diminishes the effect of fluid compressibility, as it is similar in all of the compared measurements. The rapid acquisition and analysis of the experimental results allow the use of the method in in-line quality control applications. 1 EXPERIMENT AND THEORETICAL BACKGROUND 1.1 Experimental Setup and Procedure In order to validate the novel method of permeability measurement, reference steady state measurements were performed with a Mesdan Air Tronic instrument on five fabric materials with an area of 50 cm2 at a constant pressure difference of 100 Pa. The permeability was calculated in accordance with the Darcy’ s law by using 1.84× 10 5 Pa·s as the dynamic viscosity of air at 25 ° C. The unsteady state method of determining permeability allows a prompt measurement performed in less than 50 ms and more importantly, it could also be used with arbitrarily shaped samples (for example: 3D printed materials, composite scaffolds, freez e dried materials, etc.). The unsteady state measurements of the fabric materials were performed with the experimental setup comprised from the components presented in Fig. 1. A computer and a data acquisition device are required for acquiring and analyz ing the measured signal. The air supply chamber is connected with the upstream chamber through a valve to provide the required amount of the working fluid to the upstream chamber, which is equipped with a pressure signal to determine the initial pressure p0(t= 0) as it is directly correlated with the quantity of the air in the chamber. The pressure in the upstream chamber is higher compared to the pressure in the surroundings to ensure the driving force for the working fluid is available and Fig. 1. Experimental setup for the unsteady state measurements properly oriented towards the left hand side in Fig. 1. The vacuum pump is installed at the outlet of the measured material in order to fully dry the sample material before the measurement of the permeability. After the desired amount of air is transferred to the upstream chamber, the valve connecting the supply and upstream chamber is closed and the measurements of both pressure signals p(t) and p0(t) commence. At this moment both sides of the fabric material are subjected to air at the atmospheric pressure. Afterwards, the valve connecting the upstream chamber with the fabric material is opened and the air is allowed to dissipate through the fabric as the outlet side is in contact with the atmospheric pressure in the surroundings. The pressure sensor at the inlet side of the fabric material is used to acquire the pressure signal at a high frequency during the pulse of air permeating through the sample. The measured pressure signal p(t) is in direct correlation with the sample’ s permeability. 2.2 Data Reduction The permeability of a material in a viscous flow regime is defined with the Darcy’ s law, in which the permeability of the material . is determined with the equation dV ° dt ˜..Adp . (1) dx The pressure difference across the material is dp and the volumetric flow rate is dV /dt . The area and the thickness of the material are denoted with A and dx , respectively, whereas the dynamic viscosity of the permeating fluid is . The mass transfer resistance of the material is inversely proportional to the permeability in accordance with the equation dp .dx R ˜ ˜°. (2) dV .A dt The unsteady state measurements include the same volume of the permeating fluid dV , as it depends on the volume and the initial pressure in the upstream chamber, which were both constant during the experiments. Similarly, the viscosity of the fluid and the area of the sample were also kept constant during the transient measurements, whereas the thickness of the fabric samples varied from 0.135 mm to 0.5 mm. The comparative analysis of the experimental results also allows us to exclude some of the otherwise influential parameters (for example: compressibility of air). These parameters become insignificant, as we are keeping them constant during the compared experimental runs. Therefore, in the established boundary conditions only the product of pressure difference and time dpdt is the quantity required for comparative analysis of the materials’ permeability. The tested fabric materials are subjected to a pulse of permeating air at the inlet side, whereas the fabric is in contact with atmospheric conditions at the outlet side. The unsteady pressure is measured at the inlet with an absolute pressure sensor, and the experimental results are afterwards normaliz ed as stated with the equation pt p 0 ˜°.˜° pnorm ˜°t .. (3) pmax p 0 .˜° The best results of analysis is expected at high speed measurements and large pressure differences, as the induced pressure pulse is rapidly declining when measuring highly permeable materials. The pressure signal p(t) is measured at the inlet of the fabric material, whereas the pressures p(0 ) and pmax are the initial pressure prior to the air pulse and the maximum possible value measured at the inlet. The maximum pressure value is defined as the pressure p0(0) initially set in the upstream chamber. The normaliz ation of the pressure signal improves the comparability of the experimental results acquired in various measurements as there are potential differences (i) in the initial pressure in the upstream chamber pmax and (ii) in the initial pressure p(0) at the inlet of the fabric, which is equal to the atmospheric pressure of the surroundings in the current experimental setup. Similar normaliz ation procedure was used also for the time variable in order to improve the graphical distinctness of the differences between experimental results. The normaliz ed time tnor m is defined with the equation tt pmax °.. tnorm ˜., (4) t °t .p max max in which the maximum time tmax was 70 ms, defined as the time at which the pressure pulse is diminished back to the atmospheric pressure for all of the analyz ed fabric materials. The value of t(pmax ) was derived from the experimental results separately for every experimental run at the maximum achieved pressure at the inlet of the fabric material. The comparative parameter of the permeability can be defined as the maximum pressure, the slope of the pressure rise or decay, the area underneath the pressure signal in the p(t) diagram or a combination of them all. The decision on the final comparative parameter depends on the initial and boundary conditions of the measurements and especially on the sample material’ s permeability. The performed measurements and analysis of the fabric materials exhibited the best results with comparing the following parameter PT ˜.pt .ˆt °., (5) norm norm i norm i which graphically presents the area underneath the normaliz ed pressure curve in the pnor m(tnor m) diagram. 3 RESULTS AND DISCUSSION 3.1 Fabric Materials The unsteady state measurements were performed on five different fabric materials by utiliz ing the experimental setup presented in Fig. 1 . Preliminary experiments were performed to adjust the experimental conditions. The initial pressure in the upstream chamber was varied from 0.1 bar to 0.5 bar, the number of layers of fabric materials ranged from 1 to 4, whereas the upstream chamber volume and the measurement frequency remained constant. The results from the preliminary measurements showed a relatively high permeability of the fabric samples and consequentially low signal from the inlet pressure sensor, which is in viscous flow directly dependent on the following quantities dp ˜fdx , 1 1 dV , ,, .. (6) °A ..dt Therefore, the pressure response of highly permeable fabric materials can be amended by varying their thickness or surface area perpendicular to the permeating flow, or by changing the volumetric flow, which is in our experimental setup indirectly controlled by the initial pressure in the upstream chamber. Hence, the decision on measuring 4 layers of fabric at the initial pressure of 0.5 bar in the upstream chamber was made to improve the measured pressure response. Additional rise of the initial upstream pressure, which would lower the noise to signal ratio was omitted in an attempt not to deviate even further from the pressures maintained in the referential steady state measurements, which were performed at a pressure difference of 100 Pa. The pressure signal at the inlet was measured at a 4000 Hz frequency, which is presented without any filtering or conditioning in Fig. 2. All of the measurements performed on the four layers of fabric materials were set to have the maximum pressure signal at 25 ms and the pressure decayed to the initial value of the atmospheric pressure at 70 ms. The inlet pressure sensor used was a an absolute pressure sensor, therefore some pressure differences in the initial states are observed due to the changes of the atmospheric pressure in the surroundings. This difference is negligible during the measurements on fabric materials 1, 2 and 4, as all the experiments were performed on the same day with the same atmospheric pressure. Fig. 2. Unconditioned pressure measured on five fabric materials vs. time The pressure signals as well as time were normaliz ed to further improve the comparability of the experimental results. The normaliz ation was obtained in accordance with the Eqs. (3) and (4) and the corresponding results are depicted in Fig. 3. In addition, the analyz ed data was reduced to only 141 pressure measurements for each fabric. The data was acquired at 4000 Hz from approximately 25 ms to 60 ms, as shown in Fig. 2. The graphical representation of the measurements clearly shows some pressure oscillation, especially immediately after the maximum pressure is achieved. These oscillations were observed only in measurements of thin sheet materials at a relative large initial pressure in the upstream chamber, whereas experimental results of thicker less permeable samples or at lower pressure potentials hindered the pressure fluctuations. Therefore less distinct oscillations and a more smooth pressure pulse is achievable by adding additional layers of fabric materials, which inhibits the vibrations of the material during the permeating pulse of air. Sample materials’ permeability could be compared on the basis of different experimental parameters (maximum pressure, slope of the pressure signal, etc.), however the relatively high permeabilities of the analyz ed fabric materials imposed the comparison of the area below the graphs presented in Fig. 3, which was calculated in accordance with the Eq. (5) as the P T norm. The final compared parameter included the constant pressure difference across the sample. The pressure signal in the range of 0 tnor m 0.78, hich reference permeabilities were calculated with the corresponds to the absolute time span of 25 ms t Darcy’ s law and range from approximately 8 Da to 50 60 m s. Da for the analyz ed fabrics. Fig. 3. Normalized pressure measured on five fabric materials vs. normalized time The comparison of the unsteady state measured parameter dx /PT and the permeability . measured at a steady state is presented in Fig. 4. The compared parameter is the inverse value of PT multiplied nor mnor m with the sample thickness dx as this combined value is theoretically proportional to the permeability. The presented analysis shows proportional changes in the unsteady state measurements of all five fabric materials compared to the referential steady state permeability. The comparison of the unsteady state parameter and referential permeabilities is especially encouraging as: (i) the experimental conditions and pressure sensors could be further optimiz ed; (ii) the measured results are not conditioned in any way; (iii) only one measurement was performed for each sample material; (iv) utterly different steady and unsteady state measurement protocols provide completely comparable experimental results. The complete set of the experimental results are presented in Table 1. The reference steady state measurements for all five fabric materials were performed by measuring the volumetric flow at a Table 1. Summary of the experimental results and calculations Fig. 4. Comparison of steady and unsteady state measurements The unsteady state permeability comparison parameter was calculated from the normaliz ed pressure signal and normaliz ed time step. The steady state permeability in comparison with the unsteady state permeability parameter is presented in Fig. 5 with an addition of a linear curve fit given in Eq. (7) . The unsteady state permeability derived from the linear fit as well as the relative difference compared to the reference permeability are also included in Table 1. Fig. 5. Linear interpolation of the experimental results Steady state measurements Unsteady state measurements Fabric thickness Relative difference Volumetric flow Permeability Unsteady state Calculated [mm] [%] [l/min] [Da] dx /PT permeability [Da] nor m [m] Fabric 1 Fabric 2 Fabric 3 Fabric 4 Fabric 5 0.500 0.135 0.350 0.250 0.330 59.6 153.2 37.9 139.8 231.4 18.49 12.83 8.23 21.68 47.37 0.031 0.022 0.008 0.040 0.086 17.05 12.16 4.63 22.15 48.18 -7.8 -5.2 -43.8 2.2 1.7 86 Sitar, A. The relative difference is the largest at the lowest permeability, which is anticipated, whereas the absolute differences are more comparable ranging from –3.60 Da to + 0.81 Da. The unsteady state method could be further improved by optimiz ation of measurement equipment, operating conditions and acquiring a vast number of measurement results. The latter is achieved with ease, as each measurement required less than 50 m s to conclude. The simple linear regression of the experimental results yielded the following equation with the coefficient of determination calculated at R = 0.98 2 dx ˜°. 558 5 . . (7) PTnorm The trend line of the referential steady state and unsteady state measurements was linear in the range of the results. The experimental data can be fitted with an arbitrary form of the fitting curve, if the experimental results would deviate from the linear trend. The current linear form of the regression in Eq. (7) has the coefficient 5 58.5 with a unit of Da/m multiplied with dx /PT in meters, hence the calculated permeability nor m is in Da. The presented regression model included all of the unsteady state measurements, as the objective was to test the novel measurement method, which was successful due to the promising trends achieved with a limited amount of measurements. The calculated trend line equation can be used without additional referential steady state measurements as long as the measured sample’ s permeability ranges from 8 Da to 50 Da. Any extrapolation of the linear equation, Eq. (7) , is not appropriate, as we are relying on the comparative approach. Otherwise additional referential measurements are required to derive a renewed trend line valid in the broadened range of permeabilities. 3.2 Freeze-dried Pharmaceuticals After the measuring method was evaluated on fabric materials, the transition to measuring freez e-dried materials in vials was made. The main challenges of measuring the lyophiliz ed material is the brittleness of the material, which prohibits its’ extraction from the vial. Hence, all of the measurements have to be made in-situ. WE have tested three different vials: (1) empty vial, (2) vial with a freez e-dried 5 % aqueous solution of mannitol, and (3) a vial with a freez e dried 12 % aqueous solution of mannitol. The inlet and outlet to the vial were both made through the rubber cap of the vial as seen in Fig. 6 . The initial amount of the solution in the vial was 4 ml, which occupied approximately 1/ 3 of the vial volume, as it is commonly seen in lyophiliz ed pharmaceuticals. Fig. 6. Inlet / outlet for the freeze dried material in a vial To evaluate the possibility of comparing the permeabilities of lyophiliz ed products in vials we used two different aqueous concentrations of mannitol, which result in an utterly different porosities after freez e drying. Our previous research [14] indicates that higher concentrations of mannitol (or lactose) prolong the drying times during the same lyophiliz ation protocol. Therefore, the higher 12 wt% concentration of mannitol should have a lower permeability compared to the lower 5 wt% concentration. Fig. 7 depicts the normaliz ed pressure signal measured at the inlet of the vial filled with the freez e dried material. Fig. 7. Normalized pressure signal of the vials with freeze dried materials The upstream chamber initial pressure was 1 bar during all three measurements and therefore only the material inside the vial affected the pressure signal. The results show the normaliz ed pressure was the highest during the permeation of air through the freez e dried 12 wt% mannitol, which indicates that higher concentration of mannitol has a lower permeability as anticipated. The pressure rise in the empty vial occurs due to the pressure drop of the air flowing through the inlet and outlet tubing, and is therefore not in correlation with any permeability. Consequently we have made a relative comparison to the measurement of the empty vial by using the relative value of the P T nor m,r el parameter. The derivation is in accordance with Eq. (5) and is more clearly presented in Eq. (8) bellow. PTnorm ˜empty ° PT , ˜empty °..1, normrel PTnorm ˜empty ° PTnorm ˜5 wt % ° PT ˜5w% ..134 wt °., (8) normrel , PTnorm ˜empty ° PTn˜12 wt % ° norm PT ˜12 wt % °..211 .. normrel , PTnorm ˜empty ° At this moment it is not yet possible to translate the PT parameter to actual permeabilities, as we nor m have to perform the underlying measurements of a reference material in a vial. The reference material must have a known permeability and must also be placed in a vial with the same inlet/outlet arrangement to allow a straightforward derivation of the actual permeability of the lyophiliz ed products. 5 CONCLUSIONS The presented novel unsteady state method of measuring permeability was successfully evaluated with referential standardiz ed steady state measurements performed on five different fabric materials. The presented experimental setup is capable of acquiring the entire pressure pulse signal, and the comparability parameter (slope of increasing/ decreasing pressure, pressure peak, integral parameter, or a combination of the above) is best defined after the preliminary measurements, as it depends on the properties of the measured samples. The characteristics of the fabric materials imposed the use of the integral parameter dx /PT for the comparative analysis. The validation of the measurement method yielded a linear fit curve with the R 2 of 0.98. In the range of experimental conditions, the derived linear regression equation can be used to assess the materials’ permeability measured solely with the unsteady state method. This result is extremely encouraging as only one unsteady state pressure measurement of each fabric material is included in the analysis. Moreover, the comparative analysis is based on only 141 individual pressure measurement acquired at a frequency of 4000 Hz , which correspond to the measurement time of 35 ms. A rapid assessment of the permeability enables possible in-line applications for quality control or similar. The experimental results of the freez e-dried materials in vials have demonstrated that the presented method is capable of detecting the nor m differences in permeabilities of lyophiliz ed 5 wt% and 12 wt% mannitol aqueous solutions. The results show a significant change in the comparison parameter P T nor m,r el , which was 1.34 and 2.1 1 for the 5 wt% mannitol and 12 wt% mannitol. The parameter P T nor m,r el increases with increasing concentration of mannitol. The higher concentrations of mannitol have smaller pores and lower permeabilities after freez e drying. Therefore the comparative parameter P T nor m,r el increases with the decreasing permeability of the measured material, which was expected. Several improvements of the presented novel method of measuring permeability could be implemented without any hindrance: (i) optimiz ation of experimental equipment and operating conditions; (ii) larger number of analyz ed experimental results; (iii) conditioning of the measured signal. Therefore, the future work will include: (i) measurements with reference material in a vial (ii) inclusion of multiple unsteady measurements to improve the accuracy; (iii) analysis of effects in variations of experimental setup and boundary conditions, and (iv) the use of the developed method in several additional fields and potential materials: composite materials, filters, membranes, geological samples, freez e-dried materials, 3D printed materials, artificial bones scaffold, etc. All of these applications would benefit from the analysis of the fluid permeability rate of the material under consideration. The novelty of the presented method lies in: (i) the possibility of measuring samples of arbitrary shape, non-extractable and brittle materials, (ii) the analysis of the entire pressure signal response and not just the decay during the permeation of the working fluid, (iii) the comparative analysis of the experimental results, which requires less thermodynamic properties compared to the physical modeling approach. 6 ACKNOWLEDGEMENTS The author acknowledges the Department of Textiles, Graphic Arts and Design at the Faculty of Natural Sciences and Engineering, University of Ljubljana and their head prof. Barbara Simoncic for the referential steady state measurements of the fabric materials’ air permeability. 7 REFERENCES [1] Gao, H., Li, H.A. (2016). Pore structure characterization, permeability evaluation and enhanced gas recovery techniques of tight gas sandstones. Journal of Natural Gas Science and Engineering, vol. 28, p. 536-547, DOI:10.1016/j. jngse.2015.12.018. [2] Brace, W.F., Walsh, J.B., Frangos, W.T. (1968). Permeability of granite under high pressure. Journal of Geophysical Research, vol. 73, no. 6, p. 2225-2236, DOI:10.1029/JB073i006p02225. [3] Jones, S.C. (1997). A technique for faster pulse-decay permeability measurements in tight rocks. SPE Formation Evaluation, vol. 12, no. 1, p. 19-25, DOI:10.2118/28450-PA. [4] Dicker, A.I., Smits, R.M., (1988). A practical approach for determining permeability from laboratory pressure-pulse decay measurements. International Meeting on Petroleum Engineering, p. 8, DOI:10.2118/17578-MS. [5] Cao, C., Li, T., Shi, J., Zhang, L., Fu, S., Wang, B., Wang, H. (2016). A new approach for measuring the permeability of shale featuring adsorption and ultra-low permeability. Journal of Natural Gas Science and Engineering, vol. 30, p. 548-556, DOI:10.1016/j.jngse.2016.02.015. [6] Zhong, W., Tao, G., Li, X., Kawashima, K., Kagawa, T. (2011). Determination of flow rate characteristics of porous media using charge method. Flow Measurement and Instrumentation, vol. 22, no. 3, p. 201-207, DOI:10.1016/j. flowmeasinst.2011.02.002. [7] Song, W., Yao, J., Ma, J., Li, Y., Han, W. (2018). A pore structure based real gas transport model to determine gas permeability in nanoporous shale. International Journal of Heat and Mass Transfer, vol. 126, p. 151-160, DOI:10.1016/j. ijheatmasstransfer.2018.05.012. [8] Heidari Sureshjani, M., Ahmadi, M., Fahimpour, J. (2019). A generalized approach for fast estimation of reservoir permeability distribution from analysis of pressure transient data. Journal of Petroleum Science and Engineering, vol. 182, art. ID 106240, DOI:10.1016/j.petrol.2019.106240. [9] Ghanizadeh, A., Gasparik, M., Amann-Hildenbrand, A., Gensterblum, Y., Krooss, B.M. (2014). Experimental study of fluid transport processes in the matrix system of the European organic-rich shales: I. Scandinavian Alum Shale. Marine and Petroleum Geology, vol. 51, p. 79-99, DOI:10.1016/j. marpetgeo.2013.10.013. [10] Yong, A.X.H., Aktas A., May, D., …, Willenbacher, B. (2021). Out-of-plane permeability measurement for reinforcement textiles: A benchmark exercise. Composites Part A: Applied Science and Manufacturing, vol. 148, art. ID 106480, DOI:10.1016/j. compositesa.2021.106480. [11] Swery, E.E., Allen, T., Comas-Cardona, S., Govignon, Q., Hickey, C., Timms, J., Tournier, L., Walbran, A., Kelly, P., Bickerton, S. (2016). Efficient experimental characterisation of the permeability of fibrous textiles. Journal of Composite Materials, vol. 50, no. 28, p. 4023-4038, DOI:10.1177/0021998316630801. [12] ASTM D737-18 (2015). Standard Test Method for Determining Gas Permeability Characteristics of Plastic Film and Sheeting. ASTM International, West Conshohocken. [13] ASTM D737-18 (2018). Standard Test Method for Air Permeability of Textile Fabrics. ASTM International, West Conshohocken. [14] Sitar, A., Škrlec, K., Voglar, J., Avanzo, M., Kocevar, K., Cegnar, M., Irman, Š., Ravnik, J., Hriberšek, M., Golobic, I. (2017). Effects of controlled nucleation on freeze-drying lactose and mannitol aqueous solutions. Drying Technology, vol. 36, no. 10, p. 1263-1272, DOI:10.1080/07373937.2017.1399903. [15] Arsiccio, A., Sparavigna, A.C., Pisano, R., Barresi, A.A. (2019). Measuring and predicting pore size distribution of freeze-dried solutions. Drying Technology, vol. 37, no. 4, p. 435-447, DOI:10.1080/07373937.2018.1430042. [16] Metzger, T. (2019). A personal view on pore network models in drying technology. Drying Technology, vol. 37, no. 5, p. 497­512, DOI:10.1080/07373937.2018.1512502. [17] Zhou, N., Matsumoto, T., Hosokawa, T., Suekane, T. (2010). Pore-scale visualization of gas trapping in porous media by X-ray CT scanning. Flow Measurement and Instrumentation, vol. 21, no. 3, p. 262-267, DOI:10.1016/j. flowmeasinst.2010.05.002. [18] Siebert, T., Zuber, M., Hamann, E., Baumbach, T., Karbstein, H.P., Gaukel, V. (2020). Micro-CT visualization of structure development during freeze-drying processes. Drying Technology, vol. 38, no. 3, p. 376-384, DOI:10.1080/073739 37.2019.1572619. [19] Caglar, B., Orgeas, L., Du Roscoat, S.R., Sozer, E.M., Michaud, V. (2017). Permeability of textile fabrics with spherical inclusions. Composites Part A: Applied Science And Manufacturing, vol. 99, p. 1-14, DOI:10.1016/j.compositesa.2017.03.031. [20] Zakirov, T., Galeev, A. (2019). Absolute permeability calculations in micro-computed tomography models of sandstones by Navier-Stokes and lattice Boltzmann equations. International Journal of Heat and Mass Transfer, vol. 129, p. 415-426, DOI:10.1016/j.ijheatmasstransfer.2018.09.119. [21] Puszkarz, A.K., Krucinska, I. (2018). Modeling of air permeability of knitted fabric using the computational fluid dynamics. Autex Research Journal, vol. 18, no. 4, p. 364-376, DOI:10.1515/aut-2018-0007. [22] Zhu, G., Kremenakova, D., Wang, Y., Militky, J., Mishra, R. (2015). Study on air permeability and thermal resistance of textiles under heat convection. Textile Research Journal, vol. 85, no. 16, p. 1681-1690, DOI:10.1177/0040517515573407. [23] Ngo, N.D., Tamma, K.K. (2001). Microscale permeability predictions of porous fibrous media. International Journal of Heat and Mass Transfer, vol. 44, no. 16, p. 3135-3145, DOI:10.1016/S0017-9310(00)00335-5. [24] Fissore, D., Pisano, R., Barresi, A.A. (2018). Process analytical technology for monitoring pharmaceuticals freeze-drying - A comprehensive review. Drying Technology, vol. 36, no. 15, p. 1839-1865, DOI:10.1080/07373937.2018.1440590. [25] Bobba, S., Harguindeguy, M., Colucci, D., Fissore, D. (2020). Diffuse interface model of the freeze-drying process of individually frozen products. Drying Technology, vol. 38, no. 5-6, p. 758-774, DOI:10.1080/07373937.2019.1710711. Strojniški vestnik - Journal of Mechanical Engineering 68(2022)2, 90-100 Received for review: 2021-10-07 © 2022 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2021-11-22 DOI:10.5545/sv-jme.2021.7429 Original Scientific Paper Accepted for publication: 2021-12-16 Nonlinear Coupled D yn amic Modelling of D river-seat-cab Sys tem and Biomechanical Behaviour Prediction Leilei Z hao1 – Yuewei Yu1,* – J ianhu Cao1 – Weiwei Z hou2 1 Shandong University of Technology, School of Transportation and Vehicle Engineering, China 2 Z ibo Vocational Institute, School of Automotive Engineering, China The biomechanical responses of the driver-truck system in a dynamic environment have become a significant concern in the design and control of trucks. When evaluating the riding comfort of the seat-cab system, it is necessary to predict the biomechanical responses of the driver’s different parts and directions. However, there is no reliable model and method for effectively predicting the response characteristics of the driver-seat-cab system. Aiming at such problem, firstly, based on the 7 DOF (degree-of-freedom) seated human biodynamic model established previously, a 10 DOF non-linear dynamic model of the driver-seat-cab system was created, and its vibration differential equations were established. Secondly, the vibration signals for simulation and verification were collected through the road test using a truck. Thirdly, based on the Newmark-ß integration method, the specific solution process of the model was given. The non-linear damping coefficients of the front and rear dampers for the cab suspensions were measured with a bench test. Moreover, the simulations were conducted based on the measured model parameters, taking the collected frame vibration signals as the inputs. The results show that the simulation results agree with the test results, proving that the dynamic model can effectively predict the driver’s biomechanical responses. Finally, some useful conclusions were obtained through the simulation analysis. The established model and conclusions can provide technical support for comfort evaluation, optimization design, and control of the seat-cab suspension system. Keywords: driver-seat-cab system, driver’s health, dynamic modelling, biomechanical responses Highlights • A 10 DOF non-linear dynamic model of the driver-seat-cab system was created. • The proposed dynamic model can effectively predict the driver’s biomechanical responses. • The driver’s vertical vibration caused by the road excitation mainly concentrates in 0 Hz to 10 Hz. 0 INTRODUCTION When trucks run on real roads, the vibration caused by road roughness directly affects ride comfort and human health [1]. When the vibration reaches a certain level, it will make drivers feel uncomfortable, even cause traffic accidents [2]. Thus, it is imperative to study the vibration characteristics of the driver-vehicle under the dynamic environment. Moreover, it has significant application values in truck design, dynamic simulation, the auxiliary analysis of the man-machine interface, and the prediction of the driver’ s biomechanical responses [3]. With the rapid development of the truck industry, improving ride comfort has become a significant concern of modern truck designers [4] and [5]. The seat-cab system is an essential part of modern trucks [6] and [7]. Their suspension systems directly affect the NVH (noise, vibration, and harshness) performance [8] and [9]. The optimal design of the seat-cab system through dynamic models is an important and effective means to improve ride comfort. It is the premise and basis to accurately simulate the biomechanical responses of the driver coupled with the seat-cab system for the design and optimiz ation of the cab and/or seat suspension [10]. If the biomechanical responses obtained are inaccurate, it is easy to lead to wrong optimiz ation results. The research methods of the vibration of driver-vehicle systems can be roughly divided into two categories: the experimental method and the theoretical method [11]. In recent years, with the development of testing technology, the experimental method has been widely applied. However, because of the restriction of the high cost of the test conditions, the researchers have been seeking to theoretically simulate and predict the driver’ s biomechanical responses [12]. The truck driver’ s vibration can be divided into local vibration and whole-body vibration. The local vibration refers to the vibration of the driver’ s particular parts (e.g., the head and limbs), including the vibration of the driver’ s hand or feet. The whole-body vibration refers to the vibration transmitted to the body through the driver’ s supporting surface. The driver’ s vibration mainly belongs to the whole-body vibration, which is often transmitted to drivers through their buttocks and back [13]. Generally, the whole-body vibration is the primary vibration form that causes bodily harm to drivers. In order to evaluate *Corr. Author’s Address: Shandong University of Technology, School of Transportation and Vehicle Engineering, China, yuyuewei2010@163.com the effects of the whole-body vibration on people, the International Standard Organiz ation (ISO) proposed the international standard ISO 2631 [14], which has been supplemented and revised many times. When studying the vibration of the driver-seat-cab system for trucks, the driver is usually regarded as a rigid mass [15]. Marcu [16] studied the cab suspension control using a simplified cab model with DOFs; however, he did not consider the driver-seat-cab dynamics interaction. Akç ay and Tkay [17] derived a 3 DOF (degree-of-freedom) cabin model to investigate the ride comfort of the cabin system for a truck; they regarded the driver-seat-cab as a rigid mass. Wang et al. [18] made a fatigue failure assessment based on loading reproduction for commercial vehicle cab. They also ignored the vibration characteristics of the driver-seat system. Liu et al. [19] optimiz ed the ride comfort of cab mounts but did not consider the drivers’ biomechanical characteristics. Z hao et al. [20] simulated the non-linear vibration responses of the cab system subject to suspension damper complete failure and took the driver-seat-cab as a rigid mass. These studies provide useful references for analysing the vibration of the driver-seat-cab system. However, few studies on integrated modelling and biomechanical responses simulation of the driver coupled with a seat-cab system for trucks undergoing random roads excitation exist. The present paper’ s main objective is to provide a reliable model for effectively predicting the responses characteristics of the driver-seat-cab system. The previous work [21] proposed a 7 DOF seated human biodynamic model. The test results show that the biodynamic model can effectively reproduce the seated human biodynamic responses. In this paper, as the extension of the previous work, a 10 DOF non­linear dynamic model of the driver-seat-cab system was proposed to simulate, analyse, and predict the driver’ s biomechanical responses for trucks undergoing random roads. In Section 1, combined with the 7 DOF seated human biodynamic model, a 10 DOF non-linear dynamic model of the driver-seat-cab system was created, and its vibration differential equations were established. In Section 2, the vibration signals for simulation and verification were collected. In Section 3, the specific solution process of the model was given. In Section 4, the non-linear damping coefficients of the front and rear dampers for the cab suspensions were measured. Moreover, based on the measured model parameters and the collected frame vibration signals, the model was validated. In Section 5, the biomechanical responses of the driver coupled with the seat-cab were simulated. In Section 6, some useful conclusions were obtained through the simulation analysis. 1 INTEGRATED MODELLING OF DRIVER COUPLED WITH SEAT-CAB 1.1 The 10 DOF Driver-seat-cab Model To simulate, analyse, and predict the driver’ s biomechanical responses for trucks undergoing random roads excitation, a non-linear dynamic model of the driver-seat-cab system is proposed, as shown in Fig. 1; the driver sub-model is shown in Fig. 2. This sub-model is based on the previous work [21]. The proposed model has 10 DOFs, with two DOFs for the cab ( and f), one DOF for the seat ( ), and seven cs DOFs for the driver ( 1, 2, 3, 4, 5, x b, and .). In Fig. 1, variables a and b are the horiz ontal distances from the front and rear suspensions to the cab mass centre, respectively; h is the vertical distance between the mass centres of the driver and the cab; l is the horiz ontal distance between the seat x mounting point to the cab mass centre; K and Cs are s the vertical stiffness and the damping coefficient of the seat suspension, respectively. The cab suspension system includes the front spring-damper assembly and the rear spring-damper assembly. K f and Cf are the vertical stiffness and the damping coefficient of the front spring-damper assembly, respectively. K and Cr r are the vertical stiffness and the damping coefficient of the rear spring-damper assembly, respectively. Cf and Cr nonlinearly depend on the relative motion velocities of the front damper and the rear damper, respectively. q f and q are the frame displacements at the mounting r positions of the front spring-damper assembly and the rear spring-damper assembly, respectively. The differential equation of the cab pitch vibration can be expressed as Fig. 2. The 7 DOF driver biodynamic model .... ... )] I ˜.. °[( Kz a ˜q ) Cz a ( ˜. qa fc ffc f .[( Kzb q ) Cz .b qb .. )] ...˜( . ˜. rc rrc r Kz ..z )( .zl . )] .[( Cz . scs s scs s x Kx ( ..x )( . xh .[ Cx .. )] bcs b bcs b . .k t( ..c t( ˜.. ˜.) .), (4 ) where, I is the cab moment of inertia. The differential equation of the seat vertical vibration can be expressed as mz .. ˜°Kz z ) °( . °. ) ( °Cz z ss sscs sscs °( ° Kz z ) °Cz z ( . °. ). (5 ) 1s1 1s1 The differential equation of the lower torso vertical vibration can be expressed as mz .. ˜°Kz z ) °( . °. ) °Kz z ( °Cz z ( °) 11 11s 11s 212 .. ( . °. °Cz z ( °) °Kz ° z ) °C ( zz ). (6 ) 212 41 14 4114 In Fig. 2, the driver model is composed of two sub-models: the vertical model and the fore-and-aft model [20]. The vertical model consists of five segments: head and neck, upper torso, arms, viscera, and lower torso. The spring (K 41 ) and damper (C41 ) .. connecting the upper and lower torsos represent the body spine. In the fore-and-aft model, the main body mass (m1 + m2+ m3 + m4 ) is treated as a single combined mass. The head and neck and the main body are connected by a rotational degree. These rigid masses are coupled by linear elastic and damping elements. 1.2 Non-linear Motion Equations 1.2.1 Differential Equations The differential equation of the cab vertical vibration can be expressed as mz .. ˜°Kz a °q ) °Cz a ( °. °q ) ( °.. .. ccfc ffc f °Kzb q Czb °) °( . ..°q . ( ... ) rc rrc r °Kz °z ) .( . °z . )], (1) [( Cz scs s scs s where, Cf and Cr can be further expressed as: f+ a . . f- C q . z . ° ° . .ˆ. ˜ 0 , (2) f C c f C q . z . ° ° . 0 f c The differential equation of the viscera vertical vibration can be expressed as mz ˜°Kz z ) °Cz z . ° .. ( °( . ) 22 221 221 Kz z °) °C( zz . (7 ) °( . °). 424 424 The differential equation of the arms vertical vibration can be expressed as mz .. ˜°Kz z ) °( . °. ( °Cz z ). (8) 33 334 334 The differential equation of the upper torso vertical vibration can be expressed as mz ˜°Kz z ) °( . °. ) °Kz z .. ( °Cz z ( °) 44 343 343 442 °Cz z . °) °Kz ° z ) °( ° ( . ( Cz z .. ) 44 2 4141 4141 Kz z ( . °( °) °Cz z °. ). (9 ) 545 545 The differential equation of the head and neck vertical vibration can be expressed as mz .. ˜°Kz z ) °( . °. ( °Cz z ). (10) 55 554 554 The differential equation of the whole-body longitudinal vibration can be expressed as: .. . 2 mx .. ˜ml .cos .°ml .sin .. bb 5 n 5 n .°Kx °x ) °( °x . ( Cx . ), (1 1) bbcs bbcs where, mb= m1 +m2+ m3+ m4 . The differential equation of the head and neck pitch vibration can be expressed as ° a b 2 .. ml ˜.mlx .. cos ˜.mlg sin ˜.k ( ˜.°) 5 n 5 n b5 n t  . b . r-c r C r ˆ.. ˜ C C q . z . . . 0 . . (3) c ˜° . . ). (12) r+ . c r ( q . ° . . z . t 0 In addition, the cab vertical displacement cs at The matrix Cq can be expressed as the seat installation position can be expressed as: T ˜c °lx .. (13) C z z cs The longitudinal (x -direction) displacement x of the seat backrest can be expressed as: cs x ˜°h .. cs (14) 1.2.2 Matrix Equation The displacement, velocity, and acceleration response vectors of the proposed model can be expressed as T u ˜[,z °,,,,,,,,zzzzzzx .], (15) c s12345 b . T . z . °,,,,,,,,.. . (16) u ˜[,. zzzzzzx .. .. .], c s12345 b .. T z .. ,,,,,,,,.. .... .. .. .. .. u .. ˜[,.. °zzzzzzx .]. (17) c s12345 b The displacement and velocity excitations of the proposed model can be expressed as q = q f q r], [, T (18) . = q .. ]. q [, q T (19) fr The matrix equation of the proposed model can be expressed as .. . MuCuKu Cq Kq+ . + ˜q °, (20) uu q where, M is the mass matrix, K is the system ustiffness matrix, Cu is the system damping matrix, K q is the excitation stiffness matrix, Cq is the excitation damping matrix. The matrix Cu can be expressed as C aC 00000000 ° ˆ. .. From Eqs. (21) and (22), it can be seen that both Cq and Cu contain Cf and Cr. From Eqs. (2) and (3) , it can be seen that Cf and Cr are non-linear. Thus, there are non-linear terms in Cu and Cq. Therefore, matrix Eq. (20) is a set of non-linear dynamic equations. 2 VIBRATION SIGNALS ACQUISITION .. To obtain the vibration signals for simulation and validation, according to the national standard GB/T 4970-2009 (China) [22], the road test was carried out, and the vibration signals were collected. The test was conducted on an automobile comprehensive test field. An asphalt road and a gravel road were selected as the test roads. For each road, the pavement is straight, and the longitudinal slope is not more than 1 % . A medium-siz ed truck was selected as the test vehicle. The driver is a male, and his weight is 72.0 kg. The test equipment is LMS Test.Lab test system, as shown in Fig. 3. Before the test, one accelerometer was installed on the seat surface and one accelerometer was installed on the seat backrest. In addition, two accelerometers were installed on the frame at the mounting positions of the front spring-damper assembly and the rear spring-damper assembly. The four accelerometers are PCB 356A 26 accelerometers. During the test, the driver fastened his seat belt and put his hands on the steering wheel . (22) ff ˜ q C bC 00000000 rr . ˆ ccc 000000 0 1112 13 and remained still. Moreover, the running speed was 60 km/h. The sample frequency was set as 512 Hz . The total duration of the collected signal is 120 s for each measurement. The measured vibration signals on .... ............... ° ccc 00000 cc 212223 29 t c 31 c 32 c 33 c 34 0 0 0 0 0 0 0 0 c 43 c 44 c 45 0 c 47 0 0 0 0 0 0 c c 0 c 0 0 0 5455 57 0 0 0 0 0 c 66 c 67 0 0 0 0 0 0 c 74 c 75 c 776 c 77 c 78 0 0 0 0 0 0 0 0 c 87 c 88 0 0 0 c 92 0 0 0 0 0 0 c 99 0 ° 0 c 0000000 c tt . .......... ˜ the asphalt road are shown in Fig. 4. The measured vibration signals on the gravel road are similar to ,(21) C u those on the asphalt road. where, Cf+ Cr+ Cs, c 44 = C1+ C2+ C , 41 c = 11 c = C, c = a2Cf+ b 2Cr+ l 2Cs+ h2Cb+ ct, 66 322 x c 33 = Cs+ C1, c 55 = C4+ C2, c 77 = C3+ C4+ C41 + C5, c 88 = C5, c 99 = Cb, c 21 = c 12 = aCf+ bC r+ lxCs, c 13 = c 31 =– Cs, c 23 = c 32 = lxCs, c 34 = c 43 =– C1, c 45 = c 54 =– C2, c 47 = c 74 =– C41 , c 57 = c 75 =– C4, Fig. 3. The LMS Test.Lab system for the test c 67 = c 76 =– C3, c 78 = c 87 =– C5, c 29 = c 92 = hCb. of the vibration signals acquisition 3 2 1 0 -1 -2 -3 0 20 40 60 80 100 120 Time [ can easily fulfil the needs of solving complex non­linear dynamic differential equation set in a large number of projects. Therefore, this method is used to solve the numerical solution of Eq. (20). The iterative equations used in the Newmark-ß method are as follows: . u t ˜t . °. u t ˜[( 1 ..) .. u t ˜... u ] t .t ˜t . , (23 ) u °u ˜u . t .˜[( 1 ..u ) .. ˜.u .. .]() 2 t . (24 ) ˜.˜. t t t t t t t Acceleration [ m/s2] Acceleration [ m/s 2] Acceleration [ m/s2] Acceleration [ m/s 2] ]s 1 .0 where, ut is the displacement vector at time t, u.t and u..t are the velocity vector and the acceleration vector 0.5 at time t, respectively; utt, u.t ˜°t , and u..t ˜°t are the 0.0 vectors of the displacement, the velocity, and the acceleration at time t+ t, respectively. a and ß are -0.5 constant. Based on the Newmark-ß method, combined -1 .0 with the Eq. (20) and its non-linear characteristics, ]s a) b) 2 Time [ the specific solution process was given to solve Eq. 12 8 4 0 -4 -8 -1 2 Time [ (20) for the time length [ 0, T ]. T is divided into N parts (n t = T ), meaning the integral step is t. In this study, .t was set as the reciprocal of the sampling frequency 512 Hz . The solution process based on Newmark-ß method is shown in Fig. 5. ) t.n=Load the excitation signals (T c) 12 8 4 0 -4 -8 -1 2 ]s Fig. 4. The measured acceleration signals on the asphalt road: a) the seat surface vertical acceleration, b) the seat backrest longitudinal acceleration, c) the vertical frame acceleration at the mounting position of the front spring-damper assembly, and d) the vertical frame acceleration at the mounting position of the rear spring-damper assembly Time [ d) No 3 NUMERICAL SOLUTION METHOD OF THE NON-LINEAR EQUATION SET Yes The solution of Eq. (20) is the solving problem of the non-linear differential equation set. The Newmark-ß Fig. 5. The specific solution process for the equation set of the method has good stability and high accuracy [23]. It proposed model 4 MODEL VALIDATION In order to verify the effectiveness of the proposed model and equations for the driver-seat cab system, the numerical simulations and comparison of the driver’ s vibration responses were carried out according to the specific solution process. 4.1. Model Parameters To obtain the damping coefficients of the dampers for the front suspension and the rear suspension of the cab, using the Z JS-30 hydraulic test bench produced by Changchun Test Machine R esearch Institute, the damping characteristic tests of the dampers were carried out by bench test, as shown in Fig. 6a . The master computer is shown in Fig. 6b. a) b) Fig. 6. Damper test; a) the test bench, and b) the master computer Using the bench test, a sine wave excitation with an amplitude of 30 mm and a frequency of 1.0 Hz was exerted on each damper. At the same time, the damper damping force was measured by an NS-WL1 tension-compression sensor. The damper displacement was measured with a SDVG20 displacement sensor. Based on the measurement data, the master computer automatically calculated the damping coefficients of the tested dampers, as shown in Table 1. Other parameters of the seat-cab system provided by the manufacturer, as shown in Table 2. Table 1. The values of the damping coefficients Moving stroke Front damper Rear damper Compression damping coefficient [Ns/m] 2317 1574 Rebound damping coefficient [Ns/m] 6256 4013 Table 2. The parameters values of the seat-cab system Parameter K Value Parameter Value a [m] 0.87 I [kg· m2] 751 b [m] 1.00 Cs [Ns/m] 592 h [m] 0.25 K s [N/m] 15000 l x [m] 0.03 f [N/m] 52000 mc [kg] 810 K r [N/m] 32000 According to the curve fitting method in reference [21], by minimiz ing the errors of the seat-to-head transmissibility function between the computed values using the model and the tested values, the parameters of the human biomechanics model were re-estimated. The previous model parameter values in reference [21] were selected as the initial values in this research. When the root mean square errors converged to 0.01, the parameter values were determined, as shown in Table 3. Table 3. The parameters values of the driver model Parameter Value Parameter Value C1 [Ns/m] 2375.4 K 1 [N/m] 120103.1 C2 [Ns/m] 676.8 K 2 [N/m] 5301.9 C3 [Ns/m] 144.8 K 3 [N/m] 13172.5 C4 [Ns/m] 1791.7 K 4 [N/m] 9150.8 C41 [Ns/m] 4020.2 K 41 [N/m] 128193.0 C5 [Ns/m] 979.4 K 5 [N/m] 292014.6 Cb [Ns/m] 621.2 K b [N/m] 9925.3 c t [Ns/m] 18.3 k t [N/m] 779.0 m1 [kg] 10.2 m4 [kg] 20.2 m2 [kg] 9.0 m5 [kg] 5.8 m3 [kg] 8.1 l n [m] 0.19 4.2 The Effectiveness Validation of the Proposed Model In this section, the simulations were conducted and the simulation results were verified with the tested results. For the simulations, the measured frame vertical acceleration signals at the mounting positions of the front spring-damper assembly and the rear spring-damper assembly were taken as the inputs. The simulation time T was set as 120 s. The lower torso vertical frequency-weighted R MS (root mean square) acceleration values were calculated from the simulated .. 1 and the measured .. 1 , respectively. The driver body longitudinal frequency-weighted R MS acceleration values were calculated from the simulated ..xb and the measured .. xb , respectively. A comparison of the frequency-weighted R MS acceleration values calculated from the measured data and from the simulated data is shown in Table 4. From Table 4, it can be seen that the relative errors of the driver body longitudinal frequency-weighted R MS acceleration values calculated from the simulated .. xb and the measured .. xb are less than 20 % ; the relative errors of the lower torso vertical frequency-weighted R MS acceleration values calculated from the simulated .. 1 and the measured .. 1 are less than 10 % . In the longitudinal direction (x -direction), the relative errors are larger than those in the vertical direction. Moreover, the driver body longitudinal frequency-weighted R MS acceleration values calculated from the simulated ..xb are less than those calculated from the measured .. xb . The main reasons are as follows: firstly, in the modelling process, the cab longitudinal acceleration ..x c was not considered. According to Eq. (14) , the longitudinal acceleration ..x cs of the seat backrest can be expressed as .. x ˜°h ... . In fact, ..x should be written as cs cs .. x cs ˜°h ... ... x c . Secondly, there are some errors in the 7 DOF driver biodynamic model. From the previous work [21], the model shows a higher average goodness-of-fit in the vertical direction than in the fore-and-aft direction. Moreover, in the fore-and-aft direction, the average goodness-of-fit is slightly larger than 80 % . It is noted that, although the driver model accuracy is not ideal, the fitting results show an improved fit compared with other existing human models. Table 4. A comparison of the frequency-weighted RMS acceleration values Measured Simulated Absolute Relative Road Response [m/s2] [m/s2] error [m/s2] error [%] 0.34 0.37 0.03 8.8 Asphalt .. 1 road .. xb 0.16 0.13 0.03 18.7 Gravel 1 road .. xb 0.26 0.21 0.05 19.2 By comparison, the results show that the simulation values agree with the test values. Although the proposed model should be further improved, the model is workable, and the model accuracy is acceptable. The proposed model can effectively predict the driver’ s biomechanical responses and basic dynamic characteristics for trucks running on random roads. 5 BIOMECHANICAL SIMULATION OF DRIVER COUPLED WITH SEAT-CAB 5.1 Simulation Analysis of the Time Domain Responses In the prediction and analysis of ride comfort, it is helpful to understand the changes of the human vibration in time series. The time-domain analysis of the driver’ s biomechanical responses is conducive to the design and the control of vehicle suspension systems. Moreover, it can help designers to determine the problems in design and to make corresponding improvements. The simulations were conducted using the proposed model to investigate the driver’ s biomechanical responses for trucks running on random roads, taking the collected frame vibration signals as the inputs. The simulated force and acceleration responses for trucks running on the asphalt road are shown in Fig. 7. The simulated driver’ s biomechanical responses for trucks running on the gravel road are similar to those on the asphalt road. The R MS force values and the frequency-weighted R MS acceleration values of the driver’ s different parts and different directions are shown in Table 5 and Table 6, respectively. Table 5. The RMS force values Position Asphalt road Gravel road The head for the vertical motion [N] 3.79 4.55 The upper torso for the vertical motion [N] 13.71 17.82 The arms for the vertical motion [N] 5.87 7.04 The viscera for vertical motion [N] 6.25 7.50 The lower torso for the vertical motion [N] 6.56 8.20 The whole-body for longitudinal motion [N] 6.15 7.38 The head for the pitch motion [Nm] 0.10 0.13 Table 6. The frequency-weighted RMS acceleration values ..0.56 0.51 0.05 8.9 Position Asphalt Road Gravel road The head for the vertical motion [m/s2] 0.37 0.53 The upper torso for the vertical motion [m/s2] 0.37 0.52 The arms for the vertical motion [m/s2] 0.43 0.57 The viscera for vertical motion [m/s2] 0.37 0.53 The lower torso for the vertical motion [m/s2] 0.36 0.51 The whole-body for longitudinal motion[m/s2] 0.13 0.24 The head for the pitch motion [rad/s2] 0.24 0.32 From Fig. 7, it can be seen that the time histories of the forces of the lower torso, the viscera, the arms, the upper torso, and the head are very similar. Table 5 shows that the R MS force values of the upper torso are larger than those of other organs. For trucks running on the asphalt road, the R MS force value of the upper torso is nearly twice as large as that of the lower torso. The R MS force values illustrate a similar rule for trucks running on the gravel road. Table 6 shows that the frequency-weighted R MS acceleration values of the viscera, the arms, the upper torso, and the head are all larger than that of the lower torso. For trucks running on the asphalt road, the frequency-weighted R MS acceleration value of the arms is 19.4 % larger than that of the lower torso. For trucks running on the gravel road, the frequency-weighted 40 30 20 20 0 0 Force [N] Torque [Nm] Force [N] Force [N] Force [N] -1 0 -20 -20 -4 0 -3 0 0 20 40 60 80 100 120 Time [ f) ]s 0.5 30 Time [ a) 20 10 0 -1 0 -20 0.0 -0.5 -3 0 ]s 0 20 40 60 80 100 120 Time [ g) ]s 40 40 30 Time [ b) 20 0 -20 20 10 Force [N] 0 10 - -20 -4 0 -3 0 ]s Time [ h) Time [ c) 60 3 40 2 1 ] ] m/s 2 Force [N] 20 0 -20 Acceleration [ 0 -1 -2 -4 0 -3 -6 0 ]s Time [ Time [ i) d) 30 1. 0 20 m/s 2 0.5 Force [N] 10 0 Acceleration [ 0.0 -0.5 -1 0 -20 -1. 0 ]s]s Time [ e) j) Time [ Fig. 7. The simulated vibration responses for trucks running on the asphalt road: a) the lower torso vertical force, b) the viscera vertical force, c) the arms vertical force, d) the upper torso vertical force, e) the head vertical force, f) the whole-body longitudinal force, g) the head pitch torque, h) the seat vertical force, i) the upper torso vertical acceleration, and j) the whole-body longitudinal acceleration R MS acceleration values illustrate the similar rule. The results show that the driver body can amplify the vertical vibration of the lower torso for trucks running on the asphalt road and the gravel road. 5.2 Simulation Analysis of the Frequency Domain Responses To analyse the relationship between the vibration energy and the vibration frequency of each part of the driver body, based on the simulated force responses in the time domain, the PSD (power spectral density) curves were plotted, as shown in Fig. 8. From Fig. 8, it can be seen that when the truck drives on the asphalt road, the driver’ s body vertical vibration energy caused by the road excitation mainly concentrates in the low-frequency section (0.5 Hz to 10.0 Hz ). The driver’ s body vibration energy in 200 150 the high-frequency segment (10.0 Hz to 100.0 Hz ) is almost z ero. The frequency section of the most extreme concentration for the vibration energy is about from 1.0 Hz to 5.0 Hz . The largest resonance peak appears at about 1.5 Hz . This is related to the natural frequencies of the cab system. From the vibration in Eq. (20), it can be obtained that the vertical vibration natural frequency is 1.35 Hz and the pitch vibration natural frequency is 1.70 Hz for the cab system. Figs. 8a to e show that the PSD curves of the vertical forces of the lower torso, the viscera, the arms, the upper torso, and the head are very similar. Fig. 8f illustrates that the longitudinal vibration energy of the driver’ s body concentrates in the low-frequency section from 1.0 Hz to 8.0 Hz . These simulation results are consistent with the results of the human vibration characteristics test [24] and [25]. 8 00 7 00 6 00 N 2/Hz /Hz 5 00 N2 1 00 4 00 3 00 PSD [ PSD [ 50 200 1 00 Frequency [ Hz ] Frequency [ Hz ] a) d) 80 ] /Hz 200 60 150 N2 N2 40 1 00 PSD [ PSD [ 20 50 0 -1012 -1 01 2 1010101010101010 Frequency [ Hz ] Frequency [ Hz b) ] e) ] z 1 00 140 80 1 20 /H 1 00 ] z/H 60 N2 N2 80 PSD [ 40 PSD [ 60 40 20 20 0 c) Frequency [ Hz ] f) Frequency [ Hz ] Fig. 8. The PSD curves for trucks running on the asphalt road: a) lower torso vertical force, b) viscera vertical force, c) arms vertical force, d) upper torso vertical force, e) head vertical force; f) whole-body longitudinal force Zhao, L. – Yu, Y – Cao, J. – Zhou, W. 6 CONCLUSIONs This paper aims to provide a reliable model and method for effectively predicting the responses characteristics of the driver-seat-cab system. Based on the 7 DOF seated human biodynamic model established previously, a 10 DOF non-linear dynamic model of the driver-seat-cab system was created, and its vibration differential equations were established. The vibration signals for simulation and verification were collected through the road test. Based on the Newmark-ß integration method, the specific solution process of the model was given. The non-linear damping coefficients of the front and rear dampers for the cab suspensions were measured with a bench test. The proposed model was validated, and the simulations were numerically conducted. Some useful conclusions were drawn: (1) By comparison, the results show that the simulation values agree with the test values. The model is workable, and the model accuracy is acceptable. The proposed model can effectively predict the driver’ s biomechanical responses and basic dynamic characteristics for trucks undergoing random roads excitation. (2) The time histories of the vertical forces and accelerations of the lower torso, the viscera, the arms, the upper torso, and the head are similar for trucks undergoing random roads excitation. (3) When the truck runs on the asphalt road, the driver’ s body vertical vibration energy caused by the road excitation mainly concentrates in the low-frequency section (0 Hz to 10.0 Hz ). The driver’ s body vibration energy in the high-frequency segment (10.0 Hz to 100.0 Hz ) is almost z ero. To reasonably evaluate the design effect of the cab suspension system, it is necessary to assess the cab comfort objectively. According to the ISO 2631-1 standard [14], the cab comfort is usually evaluated on the basis of the vibration accelerations at the driver’ s backrest and hip. However, the existing cab model can only simulate the vibration of the cab floor. Therefore, based on the existing cab model, the cab comfort cannot be accurately and comprehensively simulated. The proposed model in this paper can overcome the shortcomings of the existing model, and can more effectively simulate the vibration accelerations at the driver’ s backrest and hip. It can be seen that the proposed model is more helpful to guide the matching and designing of the cab suspension system than the existing model. In addition, the proposed model can also be used to analyse the driver’ s whole-body vibration when the truck is impacted, which can provide a theoretical basis for taking effective measures to avoid the injury of the driver’ s viscera. It should be noted that the model is not suitable for the vibration fatigue analysis of the human body. The shortcoming of the proposed model is that it can not accurately reproduce the real longitudinal acceleration. Thus, to overcome the shortcoming, the longitudinal DOF of the cab should be considered, and the guide mechanism of the cab suspension system should be built into the dynamic model in the following study. 7 ACKNOWLEDGEMENTS This work is supported by Nature Science Foundation of Shandong (Z R 2020ME127) and the National Natural Science Foundation of China (1 1802338) . 8 REFERENCES [1] Basaran, S., Basaran, M. (2020). Vibration control of truck cabins with the adaptive vectorial backstepping design of electromagnetic active suspension system. IEEE Access, vol. 8, p. 173056-173067, DOI:10.1109/ACCESS.2020.3025357. [2] Nataraj, M., Thillikkani, S. (2020). Failure analysis of leaf spring suspension system for heavy load truck vehicle. International Journal of Heavy Vehicle Systems, vol. 27, no. 1-2, p. 1-17, DOI:10.1504/IJHVS.2020.104413. [3] Wang, G., Li, K., Yang, J., Liang, C. (2017). Simulation and vibration isolation optimization on the cab mount of a commercial vehicle. Automotive Engineering, vol. 39, no. 9, p. 1081-1086. [4] Zhao, L.L., Yu, Y.W., Zhou, C.C., Li, X.H. (2018). Throttle slice opening size and damping characteristics of cab damper for special vehicles. Acta Armamentarii, vol. 39, no. 4, p. 645­654. [5] Nguyen, V., Jiao, R., Zhang, J. (2020). Control performance of damping and air spring of heavy truck air suspension system with optimal fuzzy control. SAE International Journal of Vehicle Dynamics, Stability, and NVH, vol. 4, no. 2, p. 179-194, DOI:10.4271/10-04-02-0013. [6] Diba, F., Esmailzadeh, E. (2018). Development of hybrid electric heavy-duty truck with self-propelled trailer. International Journal of Heavy Vehicle Systems, vol. 25, no. 2, p. 203-222, DOI:10.1504/IJHVS.2018.10011869. [7] Lee, D.Y., Elgowainy, A., Kotz, A., Vijayagopal, R., Marcinkoski, J. (2018). Life-cycle implications of hydrogen fuel cell electric vehicle technology for medium- and heavy-duty trucks. Journal of Power Sources, vol. 393, p. 217-229, DOI:10.1016/j. jpowsour.2018.05.012. [8] Yang, F.X., Zhao, L.L., Yu, Y.W., Zhou, C.C. (2017). Analytical description of ride comfort and optimal damping of cushion-suspension for wheel-drive electric vehicles. International Journal of Automotive Technology, vol. 18, no. 6, p. 1121­1129, DOI:10.1007/s12239-017-0109-2. [9] Nazemian, H., Masih-Tehrani, M. (2020). Hybrid fuzzy-PID control development for a truck air suspension system. SAE International Journal of Commercial Vehicles, vol. 13, no. 1, p. 55-69, DOI:10.4271/02-13-01-0004. [10] Zhao, L.L., Yu, Y.W., Zhou, C.C., Li; X.H., Huang, D.H. (2020). Analytical simulation of dynamic characteristics of seat system with non-linear suspension considering friction effects. International Journal of Modeling, Simulation, and Scientific Computing, vol. 11, no. 6, p. art. ID 2050058, DOI:10.1142/ S1793962320500580. [11] Rossa, F.D., Mastinu, G. (2017) Analysis of the lateral dynamics of a vehicle and driver model running straight ahead. Nonlinear Dynamics, vol. 92, p. 97-106, DOI:10.1007/ s11071-017-3478-1. [12] Maradei, M.F., Quintana, L., Castellanos-Olarte, J.M. (2016). Assessment of biomechanical demands and discomfort in drivers to stablish design criteria for truck seats. International Journal on Interactive Design and Manufacturing, vol. 10, p. 431-437, DOI:10.1007/s12008-016-0325-4. [13] Kim, J.H., Marin, L.S., Dennerlein, J.T. (2018). Evaluation of commercially available seat suspensions to reduce whole-body vibration exposures in mining heavy equipment vehicle operators. Applied Ergonomics, vol. 71, p. 78-86, DOI:10.1016/j.apergo.2018.04.003. [14] ISO 2631-1:1997. Mechanical vibration and shock-evaluation of human exposure to whole-driver vibration-part 1: general requirements. International Organization for Standardization; Geneva. [15] Stein, G.J. Múcka, P., Gunston, T.P. (2009). A study of locomotive driver’s seat vertical suspension system with adjustable damper. Vehicle System Dynamics, vol. 47, no. 3, p. 363-386, DOI:10.1080/00423110802148920. [16] Marcu, F.M. (2009). Semiactive Cab Suspension Control for Semitruck Applications. Virginia Polytechnic Institute and State University, Blacksburg. [17] Akcay, H., Trkay, S. (2014) Stochastic optimal control of truck cabin with active suspension. International Journal of Heavy Vehicle Systems, vol. 21, no. 3, p. 183-207, DOI:10.1504/ IJHVS.2014.066067. [18] Wang, T., Wang, L., Li, T., Zou, X. (2017). Fatigue failure assessment based on loading reproduction for commercial vehicle cab. Journal of Huazhong University of Science and Techology, vol. 45, no. 5, p. 61-66. [19] Liu, W., Wang, T., Shen, J., Yi, C. (2017). Ride comfort optimization of cab mounts based on response surface method. Automotive Engineering, vol. 39, no. 3, p. 323-334. [20] Zhao, L.L., Guo, J., Yu, Y.W., Li, X.H, Zhou; CC. (2020). Simulation of non-linear vibration responses of cab system subject to suspension damper complete failure for trucks. International Journal of Modeling, Simulation, and Scientific Computing, vol. 11, no. 2, p. 2050017-1~2050017-13, DOI:10.1142/S1793962320500178. [21] Gan, Z., Hillis, A.J., Darling, J. (2015). Biodynamic modelling of seated human subjects exposed to uncouples vertical and fore-and-aft whole-body vibration. Journal of Vibration Engineering and Technologies, vol. 3, no. 3, p. 301-314. [22] GB / T 4970 (2009). Method of running test-Automotive ride comfort. The Standardization Administration of China, Beijing. [23] Pourzeynali, S., Zhu, X., Zadeh, A.G., Rashidi, M., Samali, B. (2021). Comprehensive study of moving load identification on bridge structures using the explicit form of newmark-ß method: Numerical and experimental studies. Remote Sensing, vol. 13, no. 12, art. ID 2291, DOI:10.3390/rs13122291. [24] Duarte, M.L.M., de Melo, G.C. (2018). Influence of pavement type and speed on whole body vibration (WBV) levels measured on passenger vehicles. Journal of the Brazilian Society of Mechanical Sciences and Engineering, vol. 40, no. 3, DOI:10.1007/s40430-018-1057-0. [25] Huang, Y., Ferguson, N.S. (2018). Identification of biomechanical nonlinearity in whole-body vibration using a reverse path multi-input-single-output method. Journal of Sound and Vibration, vol. 419, p. 337-351, DOI:10.1016/j. jsv.2018.01.002. Strojniški vestnik - Journal of Mechanical Engineering 68(2022)2, 101-113 Received for review: 2021-09-06 © 2022 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2022-01-10 DOI:10.5545/sv-jme.2021.7389 Review Scientific Paper Accepted for publication: 2022-02-02 The Measurement of the We ar of Tie R od End Components Marek Woz niak1 – Adam R ylski2 – K rz ysz tof Sicz ek1,* 1 Lodz University of Technology, Department of Vehicles and Fundamentals of Machine Design, Poland 2 Lodz University of Technology, Institute of Materials Science and Engineering, Poland This study aimed to determine the wear intensity of the bearing mating with the ball stud in the tie rod end. The course of the vehicle driving during the period between repairs of the steering assembly was assumed. The modelled scheme of cornering was assumed to be close to the real one. The model of steering mechanism allowing determination of the load of the stud of the tie rod end as a function of the steering angle was elaborated and presented. The model of the ball stud was realized and described. The values of wear intensity for components were calculated and presented. Keywords: tie rod end, wear degree, finite element method, steering mechanism system Highlights • Wear of the steel ball end of a stud was much lower than that of its polyoxymethylene (POM) bearing in the case of the tie rod end. – 194– • For the POM bearing of the stud, the wear intensity factor k was equal to 1.38· 10 m· N2 and corresponded to the linear wear – 9 intensity Ih equal to 5.9· 10. • The linear wear intensity Ih for the POM bearing can vary in a wide range depending on the contact pressure, sliding speed, temperature, humidity, and the operating and lubricating conditions for the tie rod end assembly. • The averaged values of contact pressure between the steel ball end and its POM bearing almost linearly increased with the force F loading the steering tie rod for values of such a force below 500 N. 0 INTRODUCTION A vehicle suspension is defined as a group of elastic elements and rods connecting the axes or the individual wheels of a car with a frame or directly with a vehicle body. The purpose of suspension is to diminish shocks coming from road unevenness, assuring maximum comfort and protection of transported objects against shocks and harmful vibrations. Securing the vehicle from too strong shocks also has a significant influence on the durability of various mechanisms. One of the main components of the vehicle suspension is a tie rod end. The task of the tie rod end is to enable the wheels to turn while maintaining the vehicle in the straight line. Mondragon-Parra et al. [1] asserted that the driveline and chassis components account for about 2 % of the total mechanical losses in passenger vehicles. The repair methods for various components of the suspension system, including tie rods ends, were discussed in [2]. Chavan and Patnaik [3] stated that the tie rod may fail due to varying forces and bumping of vehicle during steering. The forces from the steering are also crucial during the static condition of the vehicle. Excessive vibration and buckling of tie rod may lead to structural failure. The authors assessed buckling strength and compared the buckling performance of tie rods for different materials. Using finite element models of the tie rod, they obtained stiffness and stress distributions in each component of the tie rod. Falah et al. [4] described the vehicle tie rod as tangentially loaded because of the vertical, lateral, and horiz ontal forces acting on the suspension system when the vehicle operates in long-term conditions. The tie rod material does not exceed the elastic limit, and these varying loading effects directly influence the performance of the car linkages like the tie rods. Many conventional vehicles operate with a steering system known as rack and pinion, which incorporates tie rods. Guiggiani [5] stated that the tie rods are slender structural elements that can be used as tie-carrying tensile loads. Tie rods are part of the steering mechanism, which operates under both tensile and compressive stress. The ratio of a typical tie rod length in the cross-section is very large, and it could buckle under the influence of compressive load. Manik et al. [6] concluded that the working strength of a tie rod is the product of the allowable working stress and the minimum cross-sectional area. Each end of a steering rack is joined to the tie rod. The pinion rotates over the slotted rack. Its function is pushing and pulling the front wheels as the steering turns either to the left or the right side of the car. The presented results showed that the distribution of deformation and stress does not exceed the yield strength value and that there are neither damages nor failure of the tie rod. However, the correctness and accuracy of calculations results are dependent on selected modelling parameters. Kim et al. [7] wrote about forging aluminium alloy tie rod ends for the vehicle steering systems. In the casting experiment, the effects of additives, Ti, B, Z r, Sr, and Mg, improve the mechanical properties. Finite element analysis was performed to define an optimal configuration of the cast preform. It was observed that the highest hardness was obtained for added 0.2 % Mg. They used the finite element model to determine an optimal configuration of the cast preform. Essienubong et al. [8] determined that vehicle tie rod should resist the vertical, lateral, and horiz ontal forces acting on the suspension system when the vehicle operates in long-term conditions. The tie rod material does not exceed the elastic limit. The tie rod can be made of nickel alloys, titanium alloy, aluminium alloy, and low alloy steel. The tie rod requires a high value of modulus of elasticity for stiffness, high fracture toughness, and resistance against wear and fatigue. Aravindaraj et al. [9] reported that a tractor connecting tie rod failed due to the overload applications. Natrayan et al. [10] studied old tie rod design and materials applied in agriculture and then analysed a modified design with various stresses and deformation characteristics under critical loading conditions using ANSYS software. Also, Oz soy and Pehlivan [11] conducted the structural analysis of a tie rod end part for a van-type vehicle using the finite element model of the body, the joint, and the bearing. The analyses for the joint assembly were carried out for various tie-end orientations. They studied the stress variations and deformation characteristics of each component for various operational loading conditions. Patil et al. [12] numerically determined the natural frequency and the static stress of the tie rod end and the rod arm for a specific vehicle. They modelled the whole mechanism as a single part and found that it could operate safely. Senniappan et al. [13] studied a steering tie rod arm of a heavy commercial vehicle. Through computer-aided engineering (CAE) analysis, the authors predicted the life of the initial design proposal of the steering tie rod arm and compared it with the actual component test life. Using the failure modes, effects, and criticality analysis (FMECA), atic et al. [14] carried out a criticality analysis of failure modes of a light commercial vehicle’ s steering tie rod joint elements. Interestingly, using the material selection and the optimiz ation technique, Park et al. [15] developed a lightweight outer tie rod for an electric vehicle. They utiliz ed aluminium alloy Al6082M and improved the tie rod structural shape using metamodel-based optimiz ation. The tie rod end consists of a ball joint. During vehicle cornering, the wear of mating surfaces of the tie rod end components, such as the ball end of the cotter pin and its bearing, occurs. The wear in tribological pairs can be modelled using various models [16] to [19]. Z wierz ycki [19] proposed a general model of wear in tribological pair, in which the volumetric wear rate w is proportional to the contact area S , sliding velocity v and linear wear intensity wear intensity Ih. The latter is affected by two parameters: factor k and pressure exponent x . The latter varies usually in a range of < 1, 2> [19] and [20]. Such two parameters are functions of material properties, geometrical characteristics of surfaces, etc. An experimental determination of these parameters is described in [21]. Some values of parameter k for various tribological contacts were given in [20] and [22]. The similar linear wear intensity wear intensity Ih with exponent x equal to 1 was utiliz ed by Burdz ik et al. [23] during studies on the tribological behaviour of contact between samples made of cast iron conducted on the pin-on-disc tribometer under dry conditions and under lubrication with engine oil. The linear dependence of the linear wear intensity Ih on the average pressure to friction force was also reported in [24]. This article aims to obtain the volumetric wear of components for the tie rod end and the parameters of assumed wear model of them under specific history (conditions) of the cornering process. During analysis, different models of the tie rod end and steering mechanism were used. The authors managed to estimate experimentally the value of the wear intensity factor k for the POM bearing of the stud, which has not been presented in the literature thus far. This factor corresponded to the linear wear intensity Ih, rarely used in the literature to quantify the wear process in the joints. The effect of lubrication occurrence in contact z one on the load sharing between the one carried by direct contact of rough asperities and the one countered by hydrodynamic effect from the grease microclines between them. The effect of squeez ing the grease film on the pressure in the grease layer between the stud ball end and its bearing was also analysed. 1 TIE ROD END The ball joint of the tie rod end is a movable connection enabling a rotationally oscillating movement of one of the connected components relative to the other. The instantaneous axis of the wheel rotation is strongly affected by the cotter pin axis. Additionally, it enables angular deflection and transmits shearing and longitudinal forces. Since the ball joints perform a rotationally oscillating movement, they are lubricated with different greases [25]. The main components of the tie rod end are presented in Fig. 1, inspired by [25], and [26]. Fig. 1. Components of the tie rod end; 1 body, 2 lower cup, 3 lower part of bearing, 4 lower ring, 5 dust boot and dust boot skirt, 6 upper ring, 7 ball stud, 8 castle nut with cotter pin, inspired by [25], and [26] 2 THE FLOW CHART OF CALCULATIONS The realiz ation of this article’ s main aim needs to assume a wear model in the contact z one between mating components of the tie rod end and to identify its parameters. This needs the knowledge of parameters influencing such wear, including sliding speed mating components of the tie rod end, their geometry, and loading at the distance for which such wear occurs. The distance and sliding velocity are functions of the number and values of steering angles of vehicle wheels, and they depend on courses of vehicle velocity and millage resulted from the vehicle driving history. Because driving history has not been registered due to costs of the appropriate equipment, it can only be estimated, as described in Chapter 3. The geometry of the contact z one between mating components, namely cotter pin and main POM bearing (Fig. 1) , is the function of outer radius D 1 /2 of the ball of the cotter pin and inner radius D 2/2 and shape of the main POM bearing. The geometry of the contact z one also depends on an applied force F , so the necessary L physical and mathematical model has been elaborated and described in a further chapter. The load of the contact z one is a function of torques due to vertical, lateral, and tractive forces loading the wheels of the vehicle and kinematic relationships between components of the vehicle steering mechanism. Torques M due vertical, M due VL lateral and M due tractive forces vary mainly with Tsteering angle dand average vehicle velocity vaver f during turning. Kinematic relationships vary with steering angle d. The load has been estimated using f the model of vehicle cornering and the modelled trajectory of the cotter pin against the main POM bearing (Fig. 1) . The information obtained from calculations is complex and of a semi-iterative nature, particularly in the case of the determination of the lateral force F . L A simplified flow chart of calculations is presented in Table 1. 3 THE COURSE OF VEHICLE DRIVING Heged et al. [27] reported that the planning of proper road vehicle trajectories became a critical and complex issue due to the enhanced interest in the automation of road transportation. During such a planning, it is necessary to consider the nonholonomic dynamics of wheeled vehicles. Only nonlinear optimiz ation-based trajectory planners are able to reach valuable results, but at the expense of enhanced computational needs sometimes up to not acceptable levels. Such problems can be fully or only partly resolved using a hybrid approach based on artificial neural networks (ANNs). Similar problems can occur when recreating the trajectories of the movement of various cars that, for example, were involved in road incidents [28] and [29], during the preparation of trainers for drivers of special or high-speed vehicles or the elaboration of driving simulators [30] to [32]. Lefè vre et al. [33] reviewed various methods for motion prediction and risk assessment for intelligent vehicles. It was emphasiz ed that the choice of a risk assessment method depended on the selected motion model. According to [34], the trajectory of a vehicle strongly depends on the interactions between the vehicle, the driving behaviour, and the road environment. The knowledge of the relationship between the trajectory and road geometry facilitates the identification of high-risk driving behaviours. Table 1. The simplified flowing chart of calculations Input Input parameters Model Output parameters Output START 1 1 Time history of jth vehicle turnings, when j=1,…, N Driver decision model (course of driving) Chapter 4 vaver( j), amax( j), R ( j) Turning number j 2 3 2 vaver( j), amax( j), R ( j) Cornering scheme Chapter 5 tu( j), vb ( j) 4 Measured D 1 , D 1 , D 2, D 2 Wear calculation Chapter 8 Volumetric wear W 1 , W 2 5 3 4 Turning number j tu( j), vb ( j) Wear distance calculation Chapter 8 Wear distance M l z ˜vjtj b .() °u () j ˜1 6 Measured D 1 , D 2 Assumed Material parameters E 1 , v1 , E 2, v2 FEM model of the ball joint Chapter 8 Contact pressure p 7 9 Load of steering rod F Assumed 2 F L 0 = mvaver( j)/R ( j), a0, a1 , a2, a3 , a4 , a6 , a7 , a1 7 Magic Tire Formula Chapter 9 af 0, ar0 8 Assumed H c , R c , hc · d c , m, f t, ., Cx , A , s, ac , ., .c , g Measured 2 10 F Turning speed vaver( j) Steering angle df M Model of steering mechanisms V torque due vertical force, M L torque due lateral force, M T torque due tractive force Chapter 6 Load of steering rod F 9 13 Lateral force F L 8 Measured 2 af 0 , ar0 Vertical force F R ( j) Steering angle calculation Chapter 9 Steering angle df 10 10 Steering angle df Assumed ar, af , tr, tf , I Model of cornering vehicle Chapter 9 ., ß 11 2 Turning speed vaver( j) 11 ß Slip angle iteration Chapter 9 af , ar 12 Assumed 12 a0, a1 , a2, a3 , a4 , a6 , a7 , a1 7 af , ar Magic Tire Formula Chapter 9 Lateral force F L 13 5 Volumetric wear W 1 , W 2 M 6 7 Wear distance l z ˜vjtj b .() °u () j ˜1 Contact pressure p Volumetric wear model Chapter 8 Volumetric wear model parameter k 1 , k 2 14 9 Load of steering rod F 14 FINISH Using the field driving experiments performed on a certain number of selected highways, X u et al. [34] studied the trajectories passenger cars took on two-lane mountain roads and determined the track patterns and their relevant risks. The new area for studies on trajectory control relates to various mobile robots [35] and unmanned vehicles [36] to [40]. According to [35], the steering control of mobile robot trajectory can relate to a wide range of turning issues: from highway driving, with negligible turning for prolonged distances, up to steering for battlefield robots, requiring agile turning. The author reviewed different steering schemes that are used for robotic vehicles and battlefield robot vehicles. Skid steering and four-wheel steering were analysed. Li et al. [41] explained that overtaking is a complex and challenging manoeuvre for road vehicles and the autonomous control of the vehicle during such a manoeuvre can enhance vehicle safety. Four-wheel independent steering (4W IS) and four-wheel independent driving (4W ID) for the autonomous electric vehicles provide redundant control, enabling better control. Z heng and Yang [42] proposed a trajectory tracking strategy based on the hierarchical control method for the 4 WIS/4W ID electric vehicle with the independently controlled rotation angle and driving torque of each wheel. In the present study, researchers attempted to recreate the trajectory of vehicle motion based on data obtained from an interview with the driver driving his car in a relatively repeatable manner on similar routes over a known time. In the analysis, the driver with his Ford Mondeo, who has driven each working day from his house to the job and returned from it, was considered. Doing so, he drives 6 0 km per day and 15,600 km per year. Additionally, about 5 to 6 weekends per year he made longer journeys of 250 km to 300 km. This gives 1250 km to 1800 km per year. After six years the total millage was equal to 93,600 + (7,500 to 10,800) km. During this 93,600 km, the road distances were practically repeatable, and they were driven with an average speed 50± 10 km/h. The route required the driver to make 20 major turns of about 90 deg. Their assumed radius was 20 m. This yielded 55 00 such turns after six years. Additionally, he had to make about 30 smaller turns about 40 deg due to overtaking or bypassing other road users. Their assumed radius was 50 m. This yielded 8250 such turns for 6 years. He had also to make small turns (100 to 500 depending on the seasons of the year – averaged 30 0) of about 1 degree because of the existing ruts on the roads. Their assumed radius was smaller than 120 m. This yielded 82, 500 such turns after six years. During 7,500 km to 10,800 km of motion, the changes of speed were more dynamic. It was assumed that the average speed was about 80± 20 kph. Due to the lack of data, it was assumed that the driver made proportionally fewer big turns, namely 440 to 630 turns with the radius of 20 m to 500 m, averaged to 275 m and an angle of 90 degrees, about 660 to 950 smaller turns with the radius equal 50 m and angle 40 degrees and about 6,600 to 9,500 small turns with the radius R less than 120 m and angle of fewer than 1 degree. It was assumed that the car was moving at the maximum speed allowed without slip during cornering (below 14 m ·s– 1). 4 THE SCHEME OF CORNERING Jiang et al. [43] proposed a human-like trapez oidal steering angle model to quantitate the relationship between objective road factors and steering wheel angles generated by human drivers. They elaborated the idealiz ed curve described by the straight curve length, the radius of curve R , and curve subtended angle a, constituting three consecutive sections: straight, circular curve, and straight. The driving speed on the ideal curve was assumed to be constant. As the steering positions of human drivers disagreed with the varying geometrical position of the curve, a series of feature distances replicating the actual human steering positions on the curve was introduced in the proposed steering angle model. In this study, the simpler scheme of each jth cornering was assumed as follows. During the first one-fifth of the turn length l (j), the driver turns the max steering wheels to the needed angle value, namely amax(j) equal to 1, 40 or 90 degrees, respectively. The next three-fifths of the turn length he keeps the steering wheels with the same constant angle value. During the last one-fifth of the turn length, the driver turns the steering wheels back to the initial position. Such an algorithm was illustrated in Fig. 2. The time of the jth cornering can be estimated from Eq. (1) . l () j .() .() jRj max max tj ˜°, u () (1) v () jv () j aver aver where j = 1, ... M is number of cornering, M total number of cornering made by the driver, R ( j) the radius of the jth cornering [ m] , and vaver( j) average speed of vehicle during the jth cornering [ m·s–1 ]. The time of turning the steering wheel forward and back is equal to 0.4 t. An average speed of the bearing relative to the ball stud during the jth cornering was estimated from Eq. (2). 2 °.() jd ° max ball vj ˜, b () (2) 04. °tju () where d is the average diameter of spherical surface bal l of the ball stud [ m] . Fig. 2. The scheme of cornering 5 THE MODEL OF STEERING MECHANISM The kinematic scheme of the steering mechanism is presented in Fig. 3. The steering angle d is the independent variable. The following parameters are constants: H = 160 mm, a= 42 °, R = 120 mm, h= c ccc 200 mm, and d = 250 m m. c The force F loading the steering rod is estimated from Eq. (3) . M °M °M V LT F ˜. (3) q .cos( .C °.°ˆ) The angle e is calculated from Eq. (4) . .1 ˜.sin( hc .sin( °c )/ q ). (4) The constant distance q is calculated from Eq. (5) . q ˜ Rh .2 Rh 2 °2 cos().. cccc c (5) The angle ß is calculated from Eq. (6) . -1 .Hq .°.ˆ)  c -sin( ˜c .sin . (6) d .c  Fig. 3. The kinematic scheme of the steering mechanism Three factors affecting the steering force are vertical load, lateral force, and tractive force. Each of these conditions produces a torque acting on the steering axis, and the torque is balanced through the steering rack. The total torque about the steering axis is a summation of the individual torques. The torque resulting from the vertical load was calculated from Eq. (7) . M V ˜°( F zl .F ) s ..sin( ..)sin() ˆzr .( F zl .F ) s ..sin( v .)cos(), ˆzr c (7) where (values readout from catalogue): F l = 47 90 N is vertical load for left wheel, F = 31 10 N vertical r load for right wheel, s = -0.006 m lateral offset at the ground, . = 7° lateral inclination angle (estimated), d steer angle [ °] , and v = 3.32± 1° caster angle. The torque resulting from the lateral force was calculated from Eqs. (8) and (9) . c where C is cornering force [ N] , F yl lateral force for fleft wheel [ N] , F yr lateral force for right wheel [ N] , and r = 0.32 m tire radius. The torque due to tractive force, calculated from Eq. (10) . M ˜( F °F ) .s , (10) T xl xr where F is tractive force for left wheel [ N] , and Fxr xl tractive force for right wheel [ N] . The cornering force C was estimated from Eq. f (1 1) . C = mv 2 ()Rj j /(). (1 1) f aver The tractive force F was estimated from Eq. (12) . x 2 F ˜mgf °05 .CAv , (12 ) xt . x aver where m = 1580 kg is mass of the vehicle, f = 0.01 rolling resistance coefficient, . = 1.29 kg·m–3 t air density, C = 0.31 drag coefficient, and x A 0.84 H W = 2.055 m 2 reference frontal area. car car 6 OPERATIONAL CONDITIONS FOR THE ANALYSED STUD OF THE STEERING ROD The ball joint can be made of AISI-SAE 5140 steel, usually heat treated by quenching and tempering with the resulted microstructure of tempered martensite with acicular ferrite [44]. Such heat treatment resulted in higher material toughness [45] and [46], as cited in Ossa et al. [44], lower fracture toughness, and mechanical strength decreasing element durability, and local reduction on the material hardness lowering the fatigue endurance limit. According to Murakami [47], as cited in Ossa et al. [44], the uniaxial fatigue strength s depends on the f Vickers hardness H as Eq. (13) . v ˜f °16. Hv .. Hv . 01 (13) The calculated Vickers hardness of the failed element was of 353 H . Using Eq. (13) , the uniaxial v fatigue strength of the material is 565± 35.3 MPa level. According to Alsaran et al. [48], as cited in Ossa et al. [44], the heat treatment of AISI-SAE 5140 steel used for suspension system ball joints resulted in a fatigue endurance limit of 416 MPa. As Alsaran’ s endurance limit does not exceed the value obtained from Eq. (13) , this value is then used as the bulk fatigue endurance limit of the analysed ball joint element. The endurance limit calculated from Eq. 13 for the acicular ferrite ..v (8) M ˜°( F .F ) r tan( ), L yl yr approximates 326± 20.4 MPa. Loading for the stud of steering rod can approach up to 500 N [49], 1540 N Fyl = Fyr = Cf /,4 (9) [50] or even 2000 N [11]. The operational conditions of the analysed stud of the steering rod were highly affected by the lubrication conditions in the contact z one between the stud ball end and its insert. The lithium greases in pure form and with various additives are often used for the lubrication of ball joints in car suspensions. Miettinen and Andersson [51] found that the lubrication significantly affected the acoustic characteristics rolling bearings. Sharma et al. [52] elaborated a theoretical model for studies of the effects of operating parameters on the energy of acoustic emission (AE) generated in rolling element bearing. The model was based on the asperities-level interaction between the mating surfaces of bearing components using the Hertz ian contact approach, statistical concepts, contact load distribution during the load z one, and lubrication effects. Patil et al. [53] examined the effects of Hertz ian deformations between the mating bodies, elasto-hydrodynamic lubrication (EHL), and elastic asperity interactions of rough surfaces on combined acoustic and dynamic characteristics of a ball bearing. Wang et al. [54] developed the dynamic models allowing considering the effects of dynamic contact forces, EHL rolling friction torque, centrifugal forces, elastic hysteresis, differential sliding friction torques on the ball motion state, and slipping in the rolling bearings. Liu et al. [55] proposed an interesting dynamic model of angular contact ball bearing (ACBB), which enabled considering the effects of elastic hysteresis, differential sliding friction torques, and EHL rolling on the ball motion state. Masjedi and Khonsari [56] related the shear friction force between the raceways and ball from: the contact force and the relative sliding speed between the raceways and ball, and the ellipse contact area dimensions, features of asperities, including their friction coefficient of asperities, their load ratio, and their ultimate shear stress coefficient, features of the oil film, including the pressure sustained by it, its average viscosity, its ultimate shear stress, and the central oil film thickness. Tong and Hong [57] elaborated a friction torque model of ACBB, including the shear friction force and the ball friction force depending on the elastic hysteresis, EHL rolling, and differential sliding friction torques. Marques et al. [58] determined that in mechanical joints, a static friction model with elimination of the discontinuity at z ero velocity was usually a good choice. For such joints operating under highly variable normal loads, a dynamic friction model was required. Fortunately, apart from participation in rallies, the dynamics of loads on tie rod ends in vehicles driving mainly on the most frequently used roads is not so large. 7 WEAR MODEL FOR THE STUD OF THE STEERING ROD In the literature on abrasive wear at the contact of sliding bodies, many papers deal with various types of radial sliding bearings [59] to [70]. Fewer articles concern abrasive wear of various ball joints, including those used in medicine [71] to [78]. Watrin et al. [75] reported that the ball joint failure resulted from the interface degradation between the steel ball pin and its POM socket. Such wear was caused by two mechanisms, the abrasion of the ball socket surface with the ball pin and plastic deformation due to the axial load. In the literature, some authors combine the results concerning the resistance to motion with the results concerning the abrasive wear in a given frictional contact. For example, numerically studying the plane elastic contact problem to determine the contact strength and tribological durability of sliding bearings, Chernets et al. [79] found that both contact angles and maximum contact pressures almost linearly depended on the load, and the durability weakened nonlinearly with load enhancement. The authors proposed the characteristic function of wear resistance of materials for discrete values of specific friction forces established during tribological tests. Other researchers [16] to [19], however, believe that although both types of results may depend on the same physicochemical properties of the mating bodies, and possibly also the presence of lubricant and impurities in the contact, they are independent of each other, particularly on the macro scale. Similarly, in this paper, we have decided not to make the model of wear in contact between steel stud and its POM bearing dependent on the friction forces occurring in such a contact z one. The abrasive wear of loaded surfaces was modelled using the Archard law, described by Eq. (14) , by Ejtehadi et al. [71], as cited by [75], for the contact interface in the POM/Carbone-based-steel couple. V kFs˜°°, (14) w Anw where V is the worn volume [ mm3 ] , F the normal wnapplied load [ N] , sw the sliding distance [ mm] , and k A –1 the wear coefficient [ mm2N] . The modified Archard law in the form of Eq. (15) was used by the other authors, including [76] to [78], as cited by [75], for the evaluation of the wear depth. h (,,) k L R ˜°.p (,,), (15) R ˜°..(,,) R ˜° w stt A w st t w st t where h(Rs, .t, ft) is the evolution of radial dimension w of the ball socket surface during one cycle [ mm] , L (Rs, .t, ft) the sliding distance of the ball socket w point during the cycle [ mm] , p(Rs, .t, ft) the normal w pressure applied on the ball socket element [ MPa] , .t, ft the angle coordinates of the sphere [ rad] , and R s the ball pin radius [ mm] . Watrin et al. [75] proposed a relation Eq. (16) between the wear coefficient k and the sliding Adistance Lw. 7 058 kA ˜°.°Lw . (16) 610 . In the present study, it was assumed that the linear wear intensity Ihi for both the steel stud and its POM bearing can be calculated using Eq. (17) . I kp 2; i ˜,, (17) ˜°12 hi ia where k 1 and k 2 are the wear intensity factors for the stud and its bearing, respectively, and pa is the averaged pressure in the contact z one between the ball stud and its bearing [ Pa] . The wear intensity factor k characteriz ing the i linear wear intensity Ihi was estimated from Eq. (18) . M W wt °˜IS °.() °() . ˜°vjtj i zcorn hi b u j ˜1 Wi ki ˜M , (18) Rp °.vjtj () °() ° ba bu j ˜ 1 aver 3 –1 Sv , where w ˜I °° is volumetric wear rate [ m· s] zhi b vaver = l /t average speed of the bearing relative to b z corn the ball stud during all analysed time of driving [ m3 · s–3 ] , t summariz ed time of driving during cornering c orn M aver [ s], l ˜v °t ˜.v () jt °() j is wear distance zbcorn b u j ˜1 represented by the assumed formula [ m] , W i volumetric wear of the ball stud (i = 1) or its bearing (i = 2) [m3], S = R / pa the contact area between the ball b b stud and its bearing [ m2] , and R reaction in contact between the ball stud and its bearing [ N] . The majority of the encountered wear laws, discussed by Mukras [80], required determining contact pressure and the sliding amounts between the mating surfaces. The contact pressure can be determined using various methods, including the popular FEM [81] to [96], the less used elastic foundation model (EFM) [83] and [97] to [103], and the boundary element method (BEM) [104] to [106], and other less popular analytical procedures, like the semi-analytical method (SAM) [107]. The common problem of unwanted stress singularities is usually prevented by applying pressure instead of single point (node) loading [85]. In this study, the FEM was utiliz ed to determine contact pressure in the contact z one between the stud ball end and its bearing. Fig. 4. FEM model of contact between rotating steel ball stud 1 and its POM bearing 2 fixed in the steel body 3 The normal contact pressure in a specific contact z one can be strongly affected by a load distribution in a complex object involving several or even several doz en components. In cases of rolling bearings or complex shapes of 3D surfaces mating under frictional conditions, the effective stiffness, sometimes in the matrix form, is utiliz ed to represent the relationship between the external load and the resultant or averaged relative displacement of the object components. In the case of a rolling bearing, such stiffness is influenced by the number of rolling components, the distribution of the load on the individual bearing elements, the presence of grease, and preload. In the case of a conformal joint, the shape of the mating surfaces can be complicated by the presence of various grooves facilitating lubricant distribution, the presence and type of lubricant, the presence of a preload in the form of a tight-fitting between stud ball end, and its bearing. In the case of worn surfaces due to abrasive wear, such tight-fitting become the loose one. Some researchers conducted studies on the load distribution and stiffness characteristics of the rolling bearings [108] to [110]. Their findings might also be useful for rod ends with rolling bearing, described, int er al ia, in [111] and [112]. For these two issues, various methods can be used, including the statics analytical method (SAM), quasi-statics numerical method, and quasi-statics analytical method (QAM). The SAM hypothetically assumed the local load distribution. The effects of the preload and the variation of the contact angle are not usually considered in the traditional SAM [113]. That method was utiliz ed by Palmgren [114], Houpert [115], and Hernot et al. [116], also in problems needed considering an effect of preload on the bearing stiffness matrix [117] to [120]. The quasi-statics numerical method hypothetically assumed the rigid body displacement of the inner ring. Such an approach introduced by Jones [121] was used by Li et al. [122], Fang et al. [123], Tang [124], and Feng et al. [125]. Neol et al. [126] gave an analytical expression of quasi-static stiffness of the rolling bearings. To limit the time requirement of the quasi-statics numerical method, the quasi-statics analytical method (QAM) utiliz ing various analytical algorithms or models were applied by Sheng et al. [127], Z hang et al. [128] and [129], Fang et al. [130], Cheng et al. [131], Jedrz ejewski and Kwasny [132], Z hang et al. [133], R en et al. [134]. The QAM can also be based on models elaborated using the Finite Element Method. Such an approach was utiliz ed by Guo and Parker [135], and by Daidié et al. [136]. Wang et al. [137] proposed an improved quasi-dynamic model for the ACBBs to obtain a more accurate load distribution of them. Liu et al. [113] utiliz ed an analytical calculation method (SAMCA) for the load distribution and stiffness of a preloaded ACBB with combined loads, including the radial and axial ones. The contact pressure was obtained from the steady-state analysis using the FEM model of the contact between the rotating steel stud ball end and its POM bearing fixed in the steel body (Fig. 4) . The material model for the stud and the body utiliz ed the linear isotropic model of steel. In such a model, the density was equal to 7800 kg/m3 , the Young modulus to 210,000 MPa and the Poisson number to 0.3 [11]. The material model for the bearing utiliz ed the linear isotropic model of POM with density of 1410 kg/ m3 , the Young modulus was equal to 3000 MPa, the Poisson ratio was equal to 0.42, r espectively [11]. The linear behaviour of steel ball and POM socket was also assumed by Watrin et al. [75] for the finite element analysis of contact pressure distribution during their studies on the ball pin and plastic socket contact in a ball joint. The grid of the tetrahedral finite elements was generated automatically by the program. Each finite element comprises nodes with three degrees of freedom, such as displacements u, uy, u along x the respective main axes X , Y , Z of the global set of coordinates. The boundary conditions were as follow: C1 is the top, and the neighbouring cylindrical surfaces of the ball stud end could move only perpendicular to the axis of the stud, C2 the cylindrical surface of the steel body was fixed, and C3 the conical surface of the stood was loaded by the constant horiz ontal force F . Such a force can reach values up to 500 N [49], 1540 N [50] or even 2000 N [11]. The grid of FEs and boundary conditions are presented in Fig. 4. While driving according to the described earlier course, the necessity to make sudden turns could occur sporadically. Therefore, it was assumed for further analysis that the load on the contact z one between the stud ball end and its bearing was of quasi-static nature. However, during analysis, the non-linear character of the contact stiffness between the stud ball end and its POM bearing was also considered, by using a non-linear computational algorithm [138] to solve the problem of contact status and contact pressure for each contact element in the contact z ones between mating surfaces of the components. The necessity to consider the nonlinear character of contact stiffness is usual for rolling bearing. According to Liu et al. [113], the rating speed had a significant effect on the contact angle and stiffness of roller bearings, especially for the high-speed range. The enhancement of macroscopic loads not only increased the contact loads but also enhanced stiffness of rolling bearing. The selection of the calculation method also influenced the obtained contact pressure values for the given values of macroscopic loads and geometrical parameters of the contacting elements [113]. An unexpectedly received calculation method (e.g., a linear algorithm assuming a linear relationship between deformation and stress) resulted in obtaining the most frequently overestimated values of contact pressures, especially at high loads. The siz e of each contact element was strictly correlated to the siz e of the directly underlying volumetric finite element generated for material lying beneath the surface on which such a contact element was formed, as described in [138]. The effect of finite element siz e on the convergence of solution was studied, assuming the maximal value of mean macroscopic contact pressure pm between the stud ball end and its bearing under loading by the force F (Fig. 4) equal to 500 N, as the criterion. Such a macroscopic contact pressure related to the area comprising both contact areas of deformed loaded asperities and the nearest adjacent areas completely separated by the grease layer. After obtaining the solution, the 3D triangular finite contact elements associated with the inner spherical surface of the POM insert were recorded together with the coordinates of the corner nodes Xe, Y and Z of the given eth finite contact element, ee the contact status CS (CS = 1, when in contact and ee CS = 0, when not in contact) of this element and the evalue of the contact pressure pe. The coordinates Xe, Ye, Z of the corner nodes allowed estimating the 3D e triangle area close to area A of the given eth finite e contact element. For each finite contact element, the area of said triangle was multiplied by the contact pressure for a given contact element and the contact status. Then these products were summed to give the value of S 1 . Moreover, the area of said triangle was multiplied by the contact status of a given finite contact element, and then such products were summed to obtain the value of S 2. Dividing the S 1 by the S 2 allowed obtaining the estimated value of mean contact pressure. A more coarse but easier estimate can be obtained using Eq. (19) . The 3D triangle contact elements with the bond option were introduced between the bearing and the body, while between the stud and this bearing, similar contact elements with frictional option were used. The surface-to-surface type with asymmetric option were utiliz ed. Therefore, the flexible target elements were applied to the stud ball end 1 and the inner surface of steel body 3. The deformable contact elements were applied to the inner and outer surface of POM bearing 2 (Fig. 4) . It was used the default values of parameters characteriz ing behaviour of contact elements [138]. The friction coefficient between the steel body and the POM bearings was assumed to be equal to 0.22 [139]. Under small slipping, the friction coefficient does not change with applied force [140]. Therefore, the friction coefficient was independent of the eventual preloading of the assembly [11]. An averaged pressure was calculated from Eq. (19) . K 1 pa ˜°pe , (19) Ke ˜1 where pe is contact pressure on the eth contact element, and K number of contact elements loaded by contact pressure (based on status number: 1 – in contact, 0 – not in contact, assigned to each contact element after solution convergence). The wear W 1 of the ball end of the steel stud was calculated as a difference between the volumes of a new ball and a worn ball bordered by its maximal tangential outer surface. It was assumed that the shape of the new ball-end volume is a sphere of diameter D 1 . The worn ball end volume is a rotated ellipse with two main diameters equal to D 1 and the third main diameter equals the difference D 1– D 1 . D 1 is the change of the measured outer diameter due wear process. Therefore, the wear of the ball end of the steel stud was estimated from Eq. (20). 2 W ˜05°.°D °.D . (20) . 1 11 The wear W 2 of the POM bearing was calculated as a difference between volumes contained within the worn and new bearing/seat and bordered by its minimal tangential inner surface. It was assumed that the shape of the new bearing/seat volume is a sphere of diameter D 2. The worn bearing/seat volume is a rotated ellipse with two main diameters equal D 2 and the third main diameter equal to the sum D 2D 2. The D 2 is the change of the measured inner diameter due the wear process occurring between the ball end of the steel stud and its POM bearing. Therefore, the wear of the ball end of the steel stud was estimated from Eq. (21) . W ˜05. °.°D 2 °.D . (21 ) 2 22 7.1 The Effect of Lubrication on the Contact Pressure in Contact Zone The effect of lubrication in contact z one under mixed friction conditions on average contact pressure is rather small compared with the case of unlubricated contact, as shown in Z hang et al. [141]. According to Otero et al. [142], determination friction forces in mixed lubrication for highly loaded non-conformal contacts requires, inter alia, a load share function determining the fraction of load supported by the lubricant film. For point contacts, Z hu and Hu [143], as cited in [142], proposed such a function in the form of Eq. (22). 121 . °.061 f ? ˜.., (22) .. .126 1037 ° where . and the specific lubricant film thickness determined from Eq. (23) [142]. h0 , (23 ) ˜. °12 .°22 where h0 the minimum lubricant film thickness [142], and s1 , s2 the root mean square (R MS) roughness of the mating surfaces 1 a nd 2, respectively [142]. For line contacts, Castro and Seabra [144], as cited in [142] proposed such a function in the form of Eq. (24) . 023 . f ˜°084.˜(24) .. According to [142], both such functions did not differ significantly. The actual spherical surfaces, both of stud ball end and of its bearing, as separate mechanical components, are characteriz ed by roughness being the component of surface texture, excluding waviness and deviation of form [145]. The rough surface can be modelled using various methods. The model directly created using the FEM was proposed by Thompson [146]. In the model, the rough surface was formed first by using the nodal or elemental set, then the height of the nodes was changed randomly, which resulted in the formation of the surface with sharp peak asperities bringing extreme stress. Another approach is the application of the random Gaussian height distribution (GHD) equation for generation of the matrix of roughness points (2D), followed by developing the 3D rough surface model. Using appropriate CAD software (CATIA, SolidWorks, among others), the cloud of rough surface points is transformed into a 3D surface form and then converted into the 3 D solid form. This approach was utiliz ed in [147] to [153]. The rough asperities can be modelled as the set of the elasto-plastically deformed hemispheres [154] to [157], hemi-ellipsoids [157] and [158] or paraboloids [159] and [160] with various distributions. The contact between the stud ball end and its bearing belongs to the group of conformal joints. Therefore, to estimate the effect of lubrication occurrence in the contact z one between stud ball end and its POM bearing on the contact pressure values therein, the following assumptions were made. The content of moisture, debris, and contaminants in the grease was constant, and all these inclusions were collected in the grease outside the contact z one, e.g., in grooves to facilitate spreading grease. This assumption was made to avoid unnecessarily complicating the model and to neglect such an effect on the wear process. The temperature in the contact area between the stud ball end and its bearing was constant and did not affect the rheological properties of the grease and the wear of the mating surfaces and the grease itself. The mapping of the mixed lubrication conditions in the contact between the stud ball end and its POM bearing was made by considering the movement of a flat, rigid, and smooth surface on a rough surface, presented as a set of model asperities (Fig. 5) . Such a movement was under constant sliding speed U . The space between such surfaces was divided into two areas with different lubrication mechanisms: the first one containing modelled asperities under boundary lubrication conditions and the second one containing lubricating microclines under hydrodynamic lubrication conditions. The load on the area of a single asperity and its environment containing microclines of the lubricant is transferred by the force from the plastic deformations of the asperity in contact with the rigid plane and the lifting force of the lubricating microclines. The layer of lubricant around the centre of the asperity top was so thin that boundary lubrication processes occurring in such an area were not subject av to the laws of hydrodynamics. The thickness h of f the layer was independent on the normal load N (Fig. 5a ). The distribution of unit pressures in the presence of this thin layer was very close to the distribution of Hertz pressure, arising at the contact of non-lubricated surfaces. In the region of the microclines, the layer of lubricant was sufficiently thick, and the principles of hydrodynamic lubrication were applicable to it. The asperities were modelled as the elastic-plastic hemispheres evenly distributed over the contact z one (Fig. 5b) . The radius R asp of such a hemisphere was equal to the roughness parameter Ra equal to 5 µm. The distance R 0 between the centres of adjacent hemispheres was equal to the half of roughness parameter R Sm equal to 0.2 mm. Fig. 5. Model of contact between a rough surface and a rigid smooth plane; a) model of the grease layer separating rough surface asperities from a flat, rigid, and smooth plane, b) model of the distribution of asperities in contact with a movable rigid plane A circular isobar was created around each asperity when the maximal thickness hsav of lubricant occurred av (Fig. 5b) . It was assumed that hs = 0.5 Ra. It was avav also assumed that hs = 0.01 hf . The normal load N 1 on the contact of a single asperity with the rigid plane (Fig. 5a ) was estimated from Eq. (25) [161], as cited in [162]. N 12.HPOM °Rasp °.asp , (25) ˜°° where R asp = Ra is the radius of the hemisphere modelling the asperity, H = 150 MPa hardness of PO MPOM [163], and aasp deformation of the centre of the contact area, called ‘ roughness susceptibility’ (Fig. 5a ), defined by Eq. (26) [162]. av av ˜z ( h .h asp °.sf ), (26) where is the height of the undeformed hemisphere modelling an asperity (Fig. 5a ). The roughness susceptibility aasp was determined from Eq. (27) [164], as cited in [162]. .HPOM 22 ˜asp °Rasp ...1 .vPOM ˆ, (27) E .POM  where E = 3000 MPa is the tensile modulus of POM POM [163], and v = 0.35 the Poisson number of PO M POM, as assumed. The number N asp of asperities being in contact with the rigid plane was estimated from Eq. (28) . . °° 2 The dependency between the lubricant thickness h and pressure p for the hydrodynamic lubrication conditions was described by average R eynold’ s Eq. (32) [165] and [166], as cited in [162]. 1 °.3 °p ˆ1 °.3 °p ˆ ˜rasp ˜hav ˜..˜hav ˜. ..2 . r °r °rr °. °. asp asp .asp .asp °hav 12 ˜eq ˜U ˜, (32 ) °rasp where U is the sliding speed rigid plane relative to the rough surface (Fig. 5) . It was assumed that it can reach values up to 0.03 m/s, and .eq 0.126 Pas the equivalent dynamic viscosity for the lithium grease determined according to the procedures described in [167]. The lithium grease rheological properties were modelled using the Herschel-Bulkley model given by Eq. (33) [168]. n ˜˜...m °. (33 ) (), 01 where t0 = 61 1.26 Pa, m3.201 Pasn , n = 0.61 6 1 are parameters of the Herschel-Bulkley model for the UM185L i2 grease [168]. The dynamic viscosity of its base oil at room temperature was assumed to be about 0.1 Pas. aver For the region aasp rasp R 0 where hs was constant, the solution of Eq. (29) took the form of Eq. D 22 N ˜ (28) (34) [161], as cited in [162]. . °° where R = 0.5 · R Sm = 0.1 mm is the distance between 0 the centres of adjacent hemispheres modelling asperities. The average thickness hav of lubricant (Fig. 5a ) was described by Eq. (29) [161], as cited in [162]. 2 .r . av asp h °., for b a  f ˆˆ1 av asp rasp asp hav ˜2 °Rasp °hf , (29) .. av hs ,fasp r for a aspR 0 where rasp is the current value of the radius of the projection of the hemisphere onto the rigid plane (Fig. 5 ), and b asp the boundary radius of the area of the plastically deformed contact z one of the hemisphere modelling asperity (Fig. 5a ); it was determined from Eq. (30) [161], as cited in [162]. b ˜ 2 °R °., (30) asp asp asp where aasp is the radius of the base of the hemisphere (Fig. 5a ); considering Eq. (29) , it can be estimated from Eq. (31) . av av a ˜2 °R °.h .h ., asp asp s f (31) asp R 0 2 4 , 12 .eq .U .R rasp . 0 pp ˜m °.Cf 1 .....cos, (34 ) a aa asp ˆasp asp  where . is the current value of angle around the axis of the hemisphere modelling an asperity (Fig. 5b) , pm mean macroscopic pressure over the entire lubricated area, Cf 1 constant factor determined from Eq. (35 ) [161], as cited in [162]. aver aver °hf .°hf ...05. .03 kasp . .. .. 05 .... .aasp ˆ.aasp ˆ Cf 1 ˜,(35) 32 2 aver aver °hf .°R .°hf .3 °R .3  00 1 1 k  .1 ....... ......asp . .aasp ˆ.aasp ˆ. aasp ˆ.3 ˆ  where k asp is constant factor determined from Eq. (36) [161], as cited in [162]. a 2 kasp ˜asp , (36 ) 2 °R °haver asp f For the region aasp rasp R 0 where hsaver varied, the solution of Eq. (29) took the form (37) [161], as cited in [162]. 12 ..U .r .ar . eq asp asp asp pp ˜m °...2 °Cf 2 ..coos, (37) a 5 .ha asp ˆaver asp  where C2 is a constant factor determined from Eq. f (38) [161], as cited in [162]. °R .2 Cf 2 ˜..0 ...1 Cf 1 21 . (38) a aver .asp ˆ°hf .2  5 ....1  kasp .aasp ˆ A normal load N 2 of a single asperity arising because of hydrodynamic phenomena accompanying the displacement of a single unevenness in relation to a rigid plane was determined from Eq. (39) [161], as cited in [162]. R 02 .N 2 ˜..°asp prdrd.. (39) asp b 0 asp The total normal load of contact z one was estimated from Eq. (40) . NN °.N . ˜.N . asp 12 (40) av The average contact pressure pasp resulted from load carried out by asperities in contact z one was estimated from Eq. (41) . °SN °N 12 °SN 4 asp ° av 2 21 p ˜°˜. (41) asp 2 2 222 2 °.DN °°b .Db .° 2 asp asp 2 asp 7.2 The Squeeze Film in Contact Zone between the Stud Ball End and Its Bearing Clarke and Bell [169] elaborated a model for the prediction of the friction in a ball joint operating under very high loads and rotations at very low speeds, including z ero speed, thus, in the mixed and boundary lubrication regimes. The model considered the effect of a squeez e film. Using this model for studies on the hydraulic motors with various surface finishes on the ball joints, they found that friction in the ball joints was strongly affected by various surface qualities. In the present study, it was assumed that after the axial load on the tie rod end, the conditions necessary for the mixed friction to occur in the contact z one of the stud ball end with its bearing are created. First, during their rotation relative to each other, the lubricant layer is pressed out of the gap between them. Then a squeez e film is formed in the fracture, accompanied by the production of a normal pressure distribution. Similarly, to the case of rotation of the stud ball end relative to its bearing under lubrication, it was also assumed that the effect of content of moisture, debris, and impurities in the grease, and the temperature in contact z one on the rheological properties of grease and the wear process was neglected. Although in fact, the contact z one is closer to the sphere segment, the research on squeez e fill was carried out using the hemisphere model, assuming that this difference does not significantly affect the obtained pressure values. To estimate values of this pressure, the model like that developed by Walicki and Walicka [170] was adopted, while the then-occurring rotational movement of the stud ball end in relation to its bearing takes place under conditions of very low rotational speed, so its effect resulting from the accompanying it the hydrodynamic interactions in the grease layer on the pressure distribution may be omitted. The model of hemispherical contact between the steel stud ball end and its POM bearing operating under a squeez e film of grease was shown in Fig. 6. For the further analysis, it was assumed that squeez ing took place at an averaged speed. Fig. 6. Model of hemispherical contact between the stud ball end and its POM bearing operating under a squeeze film of grease; C = 0.5·(D – D ), ˜sq = esq / C, ˜. °d ˜/ dt . 2 1 sq sq The pressure was estimated from Eq. (42) [170]. ˆ 4 . °(.05 °D ) °3 .1 eq 1 sq p ˜2 °°.2 C  sq ..1 ..cos sq ..   * 1 . .°l °.1, (42) 4 1 .ccos . 16 ..sq ..  * where l = 0.1 is dimensionless parameter [170], C = 0.5· (D 2 – D 1 ) the radial clearance between stud ball end and its bearing, esq the eccentricity, esq = esq /C the eccentricity ratio, and ˜. °d ˜/dt the time sq sq derivative of the eccentricity ratio esq . The average pressure was determined from Eq. (43) . ./2 p °.°d . .. eq eq 0 001 p ˜. (43) av ./2 d .eq .. 0 001 The loading carried out was estimated from Eq. (44) [170]. 42 eq °(.05 °D 1) . sq 6 .ˆsq sq N ˜2 °3 °.. C 1 .2 sq ˆ.sq  2  *2 11 sq ˆ .ln 1 .sq .6 l ° 3 .1 .. (44) .° 3 1 .sq 2  ˆ   Eq. (44) enabled estimating the average value of the time derivative of the eccentricity ratio ˜. sq for the given value of force N loading the contact z one between the stud ball end and its bearing during squeez ing of grease layer between them. This allowed determining values of pressure p from Eq. (42) and of average pressure pav from Eq. (43) . 8 LATERAL FORCES LOADING WHEELS DURING CORNERING Lateral forces loading wheels are important from the point of view of driving with a varying motion trajectory of the vehicle. Using a two-degrees of freedom (DOF) planar two-track model with a nonlinear magic tire formula [171] showed that the better adhesion capacities of tires worsened both lateral vehicle stability and increased rollover propensity, which both were highly affected by both the suspension parameters and the road excitation inputs. R ega et al. [172] reported that the vehicle steering bifurcation analysis usually utiliz ed the nonlinear autonomous vehicle model based on the two degrees of freedom (2DOF) linear vehicle model. The driving effect on steering bifurcation characteristics is omitted in such a model. However, in real driving conditions, the tires can provide various lateral forces. The steering bifurcation mechanism neglecting the driving effect cannot fully reveal the vehicle steering and driving bifurcation characteristics. In the present study, the model of vehicle cornering was assumed to be similar to the 4-wheel model used by Nguyen [173]. The model is presented in Fig. 7. The geometry of the model is characteriz ed by vehicle parameters tf and tr (the front and rear track widths) equal to 1.2 m. The distance between wheel axes is equal to L = 2850 mm, and the coordinates of the centre of the gravity are a = 1 120 mm and ar = f 1730 mm, respectively. Tire force was based solely on a slip angle. Tire forces were taken as perpendicular to the velocity at each of the individual tires. The motion of the vehicle is described by the set of Eqs. (45) to (54) [173]. After their time integration, the variables . and ß were obtained. This process was iterative, as the angle d was determined also in an iterative manner. f .d . mv ˜°F cos ˜˜.F cos ˜˜ ˆ.fr fr fl fl  .dt . .F cos ˜˜.F cos ˜˜, (45)  rrrr rrlrl d 22 ˆ1 ˆ1 I ˜.F .a .. t s05 ˆ.tan(0.5 025 in(. °.ta )) . z fr ff fr ff dt 22 ˆ1 ˆ1 F .a .. t sin(. ˆ.tan(0.5 .025 05 °.ta )) . flf f fl ff 22 ˆ1 ˆ1 ˆ. t 05 °..F . a .025 sin( . ˆrr .tan (.05 ta )) rrff rr ˆF . a 2 .. t 2 sin(. °.rl taˆ1 05 taˆ1)), . ..05 n(. (46) 05 ˆ . rlfr rr .v .sin( ˜) .af °˜fr .tan ˆ1 , (47) v .cos( ˜) ˆ05. tf ° . ˆv .sin( ˜) .af °˜fl .tan 1 .., (48) v .cos( ˜) .05. tf ° . .1 ˆv .sin( ˜) .ar ° ˜fr .tan ., (49) ˜) .. t ° v .cos( 05 .r  .1 .v .sin( ˜) .ar ° ˜fl .tan , (50) ˜) ˆ. t ° v .cos( 05 .r  v ˆsin( °) .af .˜fr .tan .1 ..f , (51) °) .. tf . v ˆcos( 05  .v ˆsin( °) .af .˜fl .tan 1 .f , (52) °) .. tf . v ˆcos( 05  ˆ1 .v .sin( °) ˆar . ˜fr .tan , (53) .v .cos( °) ˆ05. tr . where for the wheel pressure equal to 35 psi, curve-fitted parameters are the following: coefficients a0 = 1.3, a1 = –0.01 15 N–1 , a2= 0.8447, a3 = –123.6505 N, a4 = 14.273 N, a6 = –0.08073 N–1 , a7 = 0.0272 and a17 = 0.0594 [171], as cited in [173]. The steering angle d was estimated from Eq. (56) f [173]. 1 1 ˜f .tan .ar cos °Rj  .ˆ.ˆ . 1 1 tan ..af cos .ˆ.ˆ°Rj .r .f , (56) 9 RESULTS For the different values of turns radius, during vehicle cornering, some courses of the force F loading the stud of the steering rod versus the steering angle d were obtained. Fig. 8 shows the course for the average value of vehicle speed v during cornering equal to 10.55 m·s–1 and Fig. 9 does this for the average value of 12.22 m·s–1 . With the increase of the cornering radius R , the values of the force F decrease almost linearly, and the changes of the force F against the steering angle decrease. The increase of the averaged velocity v with about 7 % results in the increase of the force F with about 30 % . ˜ fl . tan ˆ 1 . v . ° sin( ) ˆ .  , (54) v a r . ° cos( ) . . t . 05 r where v = vaver( j) is the average speed of vehicle [ m·s–1], and Imass inertia moment of the vehicle relative to the vertical axis [ kg·m2] . During calculations, the average values of the slip angles were used: af = 0.5( afr + afl ) and ar = 0.5 (arr + arl ). According to Nguyen [173], the Magic Tire Formula was used in general terms of Eq. (55) [171]. The index m in equation (29) can assume the following values: f for the front wheel or r for the rear wheel, respectively. Fig. 8. The force F loading the stud vs. the steering angle d for various radiuses R of cornering with the average value F ( ˜) °( aF 2 .aF ) . m 1 zm 2 zm of vehicle speed v = 10.55 m·s–1. .1 .1 ˆa sin(2 tan( Fa )) .13 zm 4 sin a 0 tan .2 ˜m .The results of the study on the effect of a finite (( aaF .aF ) .01 zm 2 zm element siz e on the convergence of solution controlled 6 zm 7 17 m by the maximum value of the mean macroscopic ( aF .a )( 1 .a sgn( ˜)) . .1 .1 contact pressure pm between the stud and its bearing a 3 sin(2 tan( Fza 4 )) 2 zm ˜m .loaded by the force F (Fig. 4) equal to 500 N were aaF .aF 01 zm 2 zm () presented in Table 2. .1 .1 .1 a 3 sin(2 tan( Fzma 4 ))  tan ˜ , (55) 2 m  a( .aF aaF ) 01 zm 2 zm  Table 2. The effect of finite element size on the maximum contact pressure between the stud and its bearing Iteration [-] Average finite element size [m] Maximum value of mean macroscopic contact pressure between the stud and its bearing [MPa] 1 0.002 2.5848 2 0.001 2.4393 3 0.0007 2.6374 4 0.0005 3.1686 5 0.0004 3.3535 It was assumed that the average finite element siz e equal 0.0005 m can allow the proper convergence of solution without unnecessarily increasing computation time and computing power involved. Fig. 9. The force F loading the stud vs. the steering angle d for different radiuses R of cornering with the average value of vehicle speed v = 12.22 m·s–1 The resulting values of von Mises stresses in the model for the stud of steering rod loaded by force F = 300 N were presented in Fig. 10a for the complete assembly and in Fig. 10b for one without the stud. The calculated maximum value of such stresses equal to 53.1 MPa occurred in the stud. Values of such stresses were below 5.9 MPa for the body and bearing. For the comparison, under loading, the tie rod end by the force F equal to 2000 N, the von Mises stresses reached values of 21 1 MPa for the stud, and 9.1 MPa for the POM bearing [11]. This non-linear tendency may indicate a complex dependence of the von Mises stresses in the stud and bearing materials on the siz e of the ball joint, materials of its components, and the latter’ s geometric conditions of cooperation determined by the inclination angle [11]. The values of contact pressure p between rotating steel ball stud and its POM bearing for the force F = 300 N are shown in Fig. 1 1a . The maximum values of them were equal to 1.87 M Pa. The values of contact pressure p between rotating steel ball stud and its POM bearing are shown for the force F = 200 N in Fig. 1 1b, for the F = 100 N in Fig. 1 1c and for F = 50 N in Fig. 1 1d, respectively. The maximum values of contact pressure relating to those force values were 1.24, 0.63 and 0.35 MPa, respectively. The averaged values of contact pressure p against the force F are shown in Fig. 12. The averaged values of contact pressure p were equal to 0.21 MPa for the force F value equal to 300 N . For the maximum sliding speed v of the stud ball b end relative to its bearing equal to 0.03 m/s, lithium grease UM185L i2 was modelled using the Herschel-Bulkley model, with its mentioned earlier parameters t0, m1 , and n , the estimated initial grease thickness equal to 7 µm, and the assumed value equal to 0.1 Pa s for the dynamic viscosity of its base oil, the equivalent dynamic viscosity .eq obtained from procedures described in [167] was equal to 0.129 P a s. The considerations on the effect of the lubricant presence on the contact pressure distributions in the contact z ones between touching asperities of the stud ball end and its bearing, and pressure distribution in lubricant microclines surrounding such asperities were carried out for the following case chosen. The motion of the stud ball end relative to its bearing under constant rotational speed 2 rad/s and loading by constant force F = 300 N, was modelled as the motion of a rigid plane relative to the rough surface with asperities modelled using evenly distributed and elasto-plastically deformed hemispheres. Such a motion was modelled with constant sliding speed U = 0.03 m/s and under loading by force F = 300 N. During modelling of the motion, the load was carried out partly by modelled asperities being in contact and partly by hydrodynamic effects from grease microclines between these asperities. The obtained value of the average contact pressure in areas participating in carrying out load via plastically deformed asperities was equal to 150 MPa, and the average value of pressure in grease microclines was equal to 0.21 M Pa. The analysis of the squeez ing process of lithium grease in a gap between the stud ball end and its bearing under the force N = F = 300 N loading the stud ball end was carried out for various initial eccentricity esq in a range from 0 up to 7 µm. The obtained average value of the time derivative of the eccentricity ratio ˜. sq reached values in a range from 0 to 0.028, respectively. The course of the average pressure values in the squeez ed grease layer in the gap between the stud ball end and its bearing under loading by force F = 300 N versus initial eccentricity ratio esq was presented in Fig. 13. Such values can reach 0.25 MPa, which was close to the average value of pressure in grease microclines when the stud ball end rotated relative to its bearing at the maximum sliding speed equal to 0.03 m/s and was loaded by force F = 300 N. It was visible that the increase of the initial eccentricity ratio esq resulted in a quick enhancement of average pressure values in the squeez ed grease layer in the gap between the stud ball end and its bearing under loading the former by a given force F value. The measured outer diameters D of the new 1 and the worn ball end of the stud made of steel differ less than the resolution of the measurement with the passameter equal to 0.005 mm. So, the wear W 1 of the ball end was omitted, as a small value. The maximal difference between measured inner diameters D 2 of new and worn POM bearings positioned in the bushing was equal to 0.4 mm. The maximum wear W 2 of the POM bearing was of about 1.88 10 – 7 m3 . The obtained value of the wear intensity factor k 2 was equal to 1.38 10 –19 m4 ·N–2 and thus the estimated value of the linear wear intensity Ih2 for the POM bearing equal to 5.9 10 –9 . It should be noted that such values were obtained when the effect of lubrication occurrence in contact z one was included via an averaging of contact pressure values over the whole macroscopic contact area. Using only the contact pressure occurred between plastically deformed asperities in the contact z one, the obtained value of the wear intensity factor k 2 was equal to 2.45 m and thus the estimated value of 10–204 ·N–2the linear wear intensity Ih2 for the POM bearing equal to 5.5 10 –4 , which was five orders higher compared to the former one. For comparison, the results from tests on the pin-on-flat tribotester [174] for the dry contact between POM pin and a plate made of cold rolled steel AISI 42C rMo4 ground to a roughness parameter Ra = 0.20 µ m were used. During the pin reciprocal motion with stroke equal to 15 mm relative to the fixed steel plate, the average sliding speed of 0.3 m·s–1 was almost eighteen-fold higher than the one for the stud-bearing contact. Under a normal load of 100 N, the average contact pressure values were equal to 2 MPa, respectively. Tests were conducted at room temperature; however, the temperature of the contact z one was not controlled. These values were almost ten-fold higher than that for the stud-bearing contact. The linear wear intensity for the POM, corresponding to the normal load of 100 N, reached values 2.63 10 –10 MPa. Therefore, it was almost six orders lower than one for the POM bearing of the tie rod end, calculated using only the contact pressure occurring between plastically deformed asperities in the contact z one. When an effect of grease lubrication was included via an averaging of contact pressure values over the whole macroscopic contact area, the linear wear contact intensity for the POM pin was only 30-fold lower than the one for the POM bearing of the tie rod end. It was because the sliding velocity, type, and shape of relative motion of mating surfaces highly affect the linear wear intensities in different contact z ones. Also, for the comparison, the results of the tests conducted using pin-on-disc tribotester [175] for the dry contact between POM-C pin and rotating disc made of 100C r steel or AISI 3004 stainless steel were used. Such discs were polished to the roughness parameters Ra 0.1 m. During tests, the constant sliding speed was equal to 0.05 m/s and was three­fold higher than the average one for the stud-bearing contact. The contact z one was loaded by the normal force 150 N corresponding to the contact pressure 1.5 MPa, which was seven-fold higher than that for the stud-bearing contact. The average temperature in contact z one was equal to 27 ° C for the AISI 3004 stainless steel and to 30 ° C for the 100 Cr steel. The linear wear intensity for the POM, corresponding to the normal load of 150 N, reached values 0.00585 for the case of 100C r steel and 0.002925 for the case of AISI 3004 stainless steel. They were by the seven orders higher than the one for the POM bearing of the tie rod end when the effect of lubrication was included via an averaging of contact pressure values over the whole macroscopic contact area. When considering that only the contact pressure occurred between plastically deformed asperities in the contact z one, the linear wear intensity for the POM pin was ten-fold and five-fold higher, respectively, compared to the one for the POM bearing of the tie rod end. It is visible that linear wear intensity in POM-steel contact can reach values from a wide range. The mentioned intensity depends on the contact pressure [174] and [176], sliding speed [175], temperature [175], humidity [177], operating conditions, particularly relative motion manner [174] to [177]. The studied POM bearing of the tie rod ends was lubricated by the lithium grease, although its amount quickly decreased during a millage of the vehicle operating in a severe environment. The temperature and humidity also varied in a wide range depending on the technical condition of the seal in the tie rod end assembly. 10 CONCLUSIONS The wear rates of mating components of the ball joints were obtained experimentally. Wear of ball end was much lower than of the POM bearing in the case of the analysed stud. Using finite element model of the stud the steering rod, mathematical equations describing kinematic dependencies for components the steering mechanism and assumed style of cornering, it was estimated value 10 –19 m4 ·N–2 of the wear intensity factor k equal 1.38· for POM bearing of the stud. This corresponds to the linear wear intensity Ih equal to 5.9· 10 –9, which can vary in a wide range depending on the contact pressure, sliding speed, temperature, humidity, operating, and lubricating conditions for the tie rod end assembly. Obtained values of the wear intensity factor k and the linear wear intensity Ih could deviate from the actual values, as they were strongly affected by the method used for determining the average values of contact pressure in the contact z one between stud ball end and its bearing. A similar observation concerning the effect of the calculation method on the values of the mean contact pressure used in the calculation of was emphasiz ed in [18]. Obtained averaged values of contact pressure between ball end and POM bearing almost linearly increased with force loading steering rod F . The force F value equal to 300 N caused the averaged values of contact pressure p equal to 0.21 MPa. Under loading, the stud ball end by the same force F values of contact pressure were close to these occurring during squeez ing of the grease layer in the gap between the stud ball end and its bearing. Additionally, such values of contact pressure were close to these of pressure in areas of hydrodynamic effect of grease microclines when the similarly loaded stud ball end rotated relative to its bearing under maximal sliding speed equal to 0.03 m/s. Simultaneously, the contact pressure in contact z ones between plastically deformed asperities reached value of 150 M Pa. 8 REFERENCES [1] Mondragon-Parra, E., Ambrose, G. (2016). Influence of Grease in Mechanical Efficiency of Constant Velocity Joints. SAE Technical Paper, 2016-01-1132, DOI:10.4271/2016-01-1132. [2] Duffy, J.E. (2017). Modern Automotive Technology, 9th ed. The Goodheart-Willcox Company. Tinley Park, IL, USA. [3] Chavan, P.M., Patnaik, M.M.M. (2014). Performance evaluation of passenger car tie rod using numerical and theoretical approach with different materials. International Journal of Research in Engineering and Technology, vol. 3, p. 92-100, DOI:10.15623/ijret.2014.0308015. [4] Falah, A.H., Alfares, M.A., Elkholy, A.H. (2007). Failure investigation of a tie rod end of an automobile steering System. Engineering Failure Analysis, vol. 14, p. 895-902, DOI:10.1016/j.engfailanal.2006.11.045. [5] Guiggiani, M. (2014). The Science of Vehicle Dynamics. 2nd ed. Springer, Netherlands, DOI:10.1007/978-94-017-8533-4. [6] Manik, A.P., Chavan, D.S., Kavade, M.V., Umesh, S.G. (2013). FEA of tie rod of steering system of car. International Journal of Application or Innovation in Engineering & Management, vol. 2, no. 5, p. 222-227. [7] Kim, H., Seo, M., Bae, W. (2002). A study of the manufacturing of tie rod ends with casting/forging process. Journal of Material Processing Technology, vol. 125, p. 471-476, DOI:10.1016/S0924-0136(02)00323-0. [8] Essienubong, I.A., Ikechukwu, O., Ebunilo, P.O., Ememobong, E.I. (2016). Static analysis on a vehicle tie rod to determine the resulting buckling displacement. International Journal of Industrial and Manufacturing Systems Engineering, vol. 1, no. 1, p. 16-24, DOI: 10.11648/j.ijimse.20160101.13. [9] Aravindaraj, E., Natrayan, L., Santhosh, M.S., Kumar, M.S. (2018). Design and analysis of connecting tie rod assembly for automotive application. Journal of Engineering Sciences, vol. 5, p. D15-D20, DOI:10.21272/jes.2018.5(2).d3. [10] Natrayan, L., Aravindaraj, E., Santhosh, M.S., Kumar, S. (2019). Analysis and optimization of connecting tie rod assembly in Agriculture application. Acta Mechanica Malaysia, vol. 2, no. 1, p. 6-10, DOI:10.26480/amm.01.2019.06.10. [11] Ozsoy, M., Pehlivan M. K. (2015). Computer aided structral analysis of a tie rod end. Acta Physica Polonica A, vol. 128, no. 2B, p. B488-B490, DOI:10.12693/APhysPolA.128.B-488. [12] Patil, M.A., Chavan, D.S., Kavade, M.V., Ghorpade, U.S., (2013). FEA of tie rod of steering system of car. International Journal of Application or Innovation in Engineering & Management, vol. 2, no. 5, p. 222-227. [13] Senniappan, M., More, R., Bhide, S., Gowda, S. (2016). Optimization of commercial vehicle’s steering tie rod arm design based on strain life approach. SAE Technical Paper, 2016-01-0381, DOI:10.4271/2016-01-0381. [14] Catic, D., Jeremic, B., Djordjevic, Z., Miloradovic, N. (2011). Criticality analysis of the elements of the light commercial vehicle steering tie-rod joint. Strojniški vestnik - Journal of Mechanical Engineering, vol. 57, no. 6, p. 495-502, DOI:10.5545/sv-jme.2010.077. [15] Park, Y.-C., Baek, S.-K., Seo, B.-K., Kim, J.-K., Lee, K.-H. (2014). Lightweight design of an outer tie rod for an electrical vehicle. Journal of Applied Mathematics, vol. 2014, 942645, DOI:10.1155/2014/942645. [16] Meng, H.C., Ludema, K.C. (1995). Wear models and predictive equations: their form and content. Wear, vol. 181-183, p. 443-457, DOI:10.1016/0043-1648(95)90158-2. [17] Szczerek, M. (1999). Wear description from the appearance and detection point of view. Scientific Problems of Machines Operation and Maintenance, vol. 34, no. 4, p. 659-675. [18] Zwierzycki, W. (1990). The Issues of a Priori Evaluation of the Durability of Sliding Friction Nodes. Poznan University of Technology, Poznan. [19] Zwierzycki, W., (1998). Forecasting the Reliability of Wearing Parts of Machines. ITE, Radom. [20] Wozniak, M., Siczek, K., Ozuna, G., Jozwiak, P. (2018). Investigations on wear and friction in the SI engine valvetrain. Combustion Engines, vol. 175, no. 4, p. 53-64, DOI:10.19206/ CE-2018-408. [21] Zwierzycki, W., Gradkowski, M. (2000). Physical Basis for the Selection of Materials for Frictionally Mating Machine Elements. ITE, Radom. [22] Rochatka, T., Stachowiak, A., Zwierzycki, W. (2006). Wear of chosen material combinations in food fats (in the combination: three rolls - cone). Agricultural Engineering, vol. 10, no. 7(82), p. 377-382. [23] Burdzik, R., Folega, P., Lazarz, B., Stanik, Z., Warczek, J. (2012). Analysis of the Impact of Surface Layer Parameters on Wear Intensity of Friction Pairs. Archives of Metallurgy and Materials, vol. 57, p. 987-993, DOI:10.2478/v10172-012­0110-8. [24] Vyshnevskyi, O.A., Davydov, A.S. (2015). Modeling the linear wear intensity with consideration of hardness and average pressure to friction surface. Information Processing Systems, vol. 4, no. 129, pp. 64-67. [25] Wozniak, M., Siczek, K., Kubiak, P., Jozwiak, P., Siczek, K. (2018). Researches on tie rod ends lubricated by grease with TiO2 and ZrO2 nanoparticles. Journal of Physics: Conference Series, vol. 1033, no. 1, 012006, DOI:10.1088/1742­6596/1033/1/012006. [26] Muscoplat, R. (2018). Tie rod - What is a tie rod? Rick’s Free Auto Repair Advice, from https://ricksfreeautorepairadvice. com/tie-rod-tie-rod/, accessed on 221-12-31. [27] Hegeds, F., Bécsi, T., Aradi, S., Gáspár, P. (2019). Motion planning for highly automated road vehicles with a hybrid approach using nonlinear optimization and artificial neural networks. Strojniški vestnik - Journal of Mechanical Engineering, vol. 65, no. 3, p. 148-160, DOI:10.5545/sv­jme.2018.5802. [28] Subirats, P., Goyat, Y., Jacob, B., Violette, E. (2016). A new road safety indicator based on vehicle trajectory analysis. Transportation Research Procedia, vol. 14, p. 4267-4276, DOI:10.1016/j.trpro.2016.05.398. [29] Houenou, A., Bonnifait, P., Cherfaoui, V., Wen, Y. (2013). Vehicle Trajectory Prediction based on Motion Model and Maneuver Recognition. IEEE/RSJ International Conference on Intelligent Robots and Systems, Tokyo, p. 4363-4369, DOI:10.1109/IROS.2013.6696982. [30] Aykent, B., Merienne, F., Guillet, C., Paillot, D., Kemeny, A. (2014). Motion sickness evaluation and comparison for a static driving simulator and a dynamic driving simulator. Journal of Automobile Engineering, vol. 228, no. 7, p. 818­829, DOI:10.1177/0954407013516101. [31] Michael-Grigoriou, D., Kleanthous, M., Savva, M., Christodoulou, S., Pampaka, M., Gregoriades, A. (2014). Impact of immersion and realism in driving simulator studies. International Journal of Interdisciplinary Telecommunications and Networking, vol. 6, no. 1, p. 10-25, DOI:10.4018/ ijitn.2014010102. [32] Dziuda, L., Biernacki, M., Baran, P., Trszczynski, O. (2014). The effects of simulator fog and motion on simulator sickness in a driving simulator and the duration of after-effects. Applied Ergonomics, vol. 45, no. 3, p. 406-412, DOI:10.1016/j. apergo.2013.05.003. [33] Lefèvre, S., Vasquez, D., Laugier, C. (2014). A survey on motion prediction and risk assessment for intelligent vehicles. Robomech Journal, vol. 1, p. 1, DOI:10.1186/s40648-014­0001-z. [34] Xu, J., Luo, X., Shao, Y.M. (2018). Vehicle trajectory at curved sections of two-lane mountain roads: a field study under natural driving conditions. European Transport Research Review, vol. 10, p. 12, DOI:10.1007/s12544-018-0284-x. [35] Lakkad, S., (2004). Modeling and Simulation of Steering Systems for Autonomous Vehicles. PhD Thesis. Florida State University, USA. [36] Chebly, A. (2017). Trajectory planning and tracking for autonomous vehicles navigation. Automatic Control Engineering. PhD Thesis. Université de Technologie de Compiègne, D2392. [37] Wang, R., Li, G., Zhang, S., Liu, L., Zhang, X., Yin, Y. (2021). Trajectory tracking control in real-time of dual-motor-driven driverless racing car based on optimal control theory and fuzzy logic method. Complexity, vol. 2021, p. 5549776, DOI:10.1155/2021/5549776. [38] Li, B., Du, H., Li, W., Zhang, B. (2018). Dynamically integrated spatiotemporal-based trajectory planning and control for autonomous vehicles. Faculty of Engineering and Information Sciences - Papers: Part B, p. 2124. . [39] Ding, W., Shen, S., (2019). Online vehicle trajectory prediction using policy anticipation network and optimization-based context reasoning. International Conference on Robotics and Automation, p. 9610-9616, DOI:10.1109/ICRA.2019.8793568. [40] Kim, B.D., Kang, C.M., Kim, J., Lee, S.H., Chung, C.C., Choi, J.W. (2018). Probabilistic vehicle trajectory prediction over occupancy grid map via recurrent neural network. IEEE 20th International Conference on Intelligent Transportation Systems, p. 399-404, DOI:10.1109/ITSC.2017.8317943. [41] Li, B., Du, H., Li, W. (2016). Trajectory control for autonomous electric vehicles with in-wheel motors based on a dynamics model approach. Institution of Engineering and Technology (IET) Intelligent Transport Systems, vol. 10, p. 318-330. https://doi.org/10.1049/iet-its.2015.0159. [42] Zheng, H., Yang, S. (2019). A trajectory tracking control strategy of 4WIS/4WID electric vehicle with adaptation of driving conditions. Applied Science, vol. 9, p. 168, DOI:10.3390/app9010168. [43] Jiang, H., Zhou, J., Li, A., Zhou, X., Ma, S. (2019). Human-like trapezoidal steering angle model on two-lane urban curves. International Journal of Advanced Robotic Systems, vol. July 2019, DOI:10.1177/1729881419867614. [44] Ossa, E.A., Palacio, C.C., Paniagua M. (2011) Failure analysis of a car suspension system ball joint. Engineering Failure Analysis, vol. 18, no. 5, p. 1388-1394, DOI:10.1016/j. engfailanal.2011.03.013. [45] Bayrak, M., Ozturk, F., Demirezen, M., Evis, Z. (2007). Analysis of tempering treatment on material properties of DIN 41Cr4 and DIN 42CrMo4 steels. Journal of Materials Engineering and Performance, vol. 16, no. 5, p. 597-600, DOI:10.1007/ s11665-007-9043-1. [46] Bol’shakov, V.I., Laukhin, D.V., Sukhomlin, G.D., Kuksenko, V.I. (2004). Effect of heat treatment on formation of acicular ferrite and on the properties of low carbon microalloyed steels 10G2FB and 09G2S. Metal Science and Heat Treatment, vol. 46, no. 11-12, p. 545-50, DOI:10.1007/s11041-005-0016-4. [47] Murakami, Y. (2002). Metal fatigue: effects of small defects and non-metallic inclusions. Elsevier, Kidlington, Oxfordshire. [48] Alsaran, A., Karakan, M., Celik, A. (2002). The investigation of mechanical properties of ion-nitrided AISI 5140 low-alloy steel. Materials Characterization, vol. 48, no. 4, p. 323-327, DOI:10.1016/S1044-5803(02)00275-9. [49] Rayu, Y.I., Kang, D.O., Heo, S.J., Yim, H.J., Jeon, J.I. (2010). Development of analytical process to reduce side load in strut-type suspension. Journal of Mechanical Science and Technology, vol. 24, no. 1, p. 351-356, DOI:10.1007/s12206­009-1103-z. [50] Lee, K.-H., Kim, J.K., Kim, Y.J., Yang, W.H., Park, Y.C. (2011). Structural design of an outer tie rod for a passenger car. International Journal of Automotive Technology, vol. 12, 375, DOI:10.1007/s12239-011-0044-6. [51] Miettinen, J., Andersson, P. (2000). Acoustic emission of rolling bearings lubricated with contaminated grease. Tribology International, vol. 33, no. 11, p. 777-787, DOI:10.1016/S0301-679X(00)00124-9. [52] Sharma, R.B., Parey, A. (2019). Modelling of acoustic emission generated in rolling element bearing. Applied Acoustics, vol. 144, p. 96-112, DOI:10.1016/j.apacoust.2017.07.015. [53] Patil, A.P., Mishra, B.K., Harsha. S.P. (2020). Vibration based modelling of acoustic emission of rolling element bearings. Journal of Sound and Vibration, vol. 468, p. 115117, DOI:10.1016/j.jsv.2019.115117. [54] Wang, Y.L., Wang, W.Z., Zhang, S.G. (2018). Effects of raceway surface roughness in an angular contact ball bearing. Mechanism and Machine Theory, vol. 121, p. 198-212, DOI:10.1016/j.mechmachtheory.2017.10.016. [55] Liu, J., Xu, Y., Pan, G. (2021). A combined acoustic and dynamic model of a defective ball bearing. Journal of Sound and Vibration, vol. 501, p. 116029, DOI:10.1016/j. jsv.2021.116029. [56] Masjedi, M., Khonsari, M.M. (2015). An engineering approach for rapid evaluation of traction coefficient and wear in mixed EHL. Tribology International, vol. 92, p. 184-190, DOI:10.1016/j.triboint.2015.05.013. [57] Tong, V.C., Hong, S.W. (2018). Improved formulation for running torque in angular contact ball bearings. International Journal of Precision Engineering and Manufacturing, vol. 19, p. 47-56, DOI:10.1007/s12541-018-0006-2. [58] Marques, F., Flores, P., Claro, J.C., Lankarani, H. (2018). Modeling and analysis of friction including rolling effects in multibody dynamics: A review. Multibody System Dynamics, vol. 45, no. 4, p. 223-244, DOI:10.1007/s11044-018-09640­6. [59] Dykha, A., Sorokatyi, R., Makovkin, O., Babak, O. (2017). Calculation-experimental modeling of wear of cylindrical sliding bearings. Eastern-European Journal of Enterprise Technologies, vol. 5, no. 1, p. 51-59, DOI:10.15587/1729­4061.2017.109638. [60] Dykha, A., Marchenko, D. (2018). Prediction the wear of sliding bearings. International Journal of Engineering & Technology, vol. 7, no. 2.23, p. 4-8, DOI:10.14419/ijet.v7i2.23.11872. [61] Rezaei, A., Ost, W., Van Paepegem, W., De Baets, P., Degrieck, J. (2011). Experimental study and numerical simulation of the large-scale testing of polymeric composite journal bearings: Three-dimensional and dynamic modelling. Wear, vol. 270, no. 7-8, p. 431-438, DOI:10.1016/j.wear.2010.11.005. [62] Rezaei, A., Ost, W., Van Paepegem, W., De Baets, P., Degrieck, J. (2012). A study on the effect of the clearance on the contact stresses and kinematics of polymeric composite journal bearings under reciprocating sliding conditions. Tribology International, vol. 48, p. 8-14, DOI:10.1016/j. triboint.2011.06.031. [63] Rezaei, A., Van Paepegem, W., De Baets, P., Ost, W., Degrieck, J. (2012). Adaptive finite element simulation of wear evolution in radial sliding bearing. Wear, vol. 296, no. 1-2, p. 660-671, DOI:10.1016/j.wear.2012.08.013. [64] Sorokatyi, R.V. (2003). Solution of the problem of wear of a fine elastic layer with a rigid bearing mounted on a rigid shaft using the method of triboelements. Journal of Friction and Wear, vol. 24, no. 1, p. 35-41. [65] Sorokatyi, R.V. (2003). Evaluation of efficiency of sliding bearings during reciprocation. Journal of Friction and Wear, vol. 24, no. 2, p. 136-143. [66] Chernets, M.V., Zhydyk, V.B., Chernets’, Yu.M. (2014). Accuracy of evaluation of the service life of a plain bearing according to the generalized cumulative model of wear. Materials Science, vol. 50, p. 39-45, DOI:10.1007/s11003- 014-9689-4. [67] Chernets, M., Chernets, Ju. (2015). Generalized method for calculating the durability of sliding bearings with technological out-of-roundness of details. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 229, no. 2, p. 216-226, DOI:10.1177/1350650114554242. [68] Chernets, M.V. (2015). Prediction of the life of a sliding bearing based on a cumulative wear model taking into account the lobing of the shaft contour. Journal of Friction and Wear, vol. 36, p. 163-169, DOI:10.3103/S1068366615020038. [69] Sep, J., Galda, L., Oliwa, R., Dudek, K. (2020). Surface layer analysis of helical grooved journal bearings after abrasive tests. Wear, vol. 448-449, p. 203233, DOI:10.1016/j. wear.2020.203233. [70] Knig, F., Marheineke, J., Jacobs, G., Sous, C., Zuo, M.J., Tian, Z. (2021). Data-driven wear monitoring for sliding bearings using acoustic emission signals and long short-term memory neural networks. Wear, vol. 476, p. 203616, DOI:10.1016/j. wear.2021.203616. [71] Ejtehadi, M.H., Klaus, H., Sommer, S., Haensel, H., Scholten, J. (2011). Running-in phase of spherical chassis joints - identification of the main influence parameter and implementation in a wear simulation tool. International Journal of Advanced Manufacturing Technology, vol. 55, no. 9-12, 983-995, DOI:10.1007/s00170-010-3136-y. [72] Kadhim, N.A., Abdullah, S., Ariffin A. (2012). Effective strain damage model associated with finite element modelling and experimental validation. International Journal of Fatigue, vol. 36, no. 1, p. 194-205, DOI:10.1016/j.ijfatigue.2011.07.012. [73] Mao, J., Luo, Y., Liu, J. (2014). Dynamics performance and abrasive wear of the automotive drive shaft. Advances in Mechanical Engineering, vol. 6, p. 713824-713824, DOI:10.1155/2014/713824. [74] Scholten, J., Haensel, H., Krekeler, N., Fuchs, H., Stenke, R., Ejtehadi, M.H. (2010). Modellierung des Einflusses der Verschleissverteilung auf die Beanspruchung von Fachwerksgelenken. ‘Inhalt’ Materials Testing, vol. 52, no. 7-8, p. 463-469, DOI:10.3139/120.010072. [75] Watrin, J., Makich, H., Haddag, B., Nouari, M., Grandjean, X. (2017). Analytical modelling of the ball pin and plastic socket contact in a ball joint. Congrès français de mécanique, Lille. [76] Wu, S., Hung, J., Shu, C. (2003). The computer simulation of wear behavior appearing in total hip prosthesis. Computer Methods and Programs in Biomedicine, vol. 70, no. 1, p. 81­91, DOI:10.1016/s0169-2607(01)00199-7. [77] Pietrabessa, R., Raimondi, M., Di Martino, E. (1998). Wear of polyethylenecups in total hip arthroplasty: a parametric mathematical model. Medical Engineering & Physics, vol. 20, no. 3, p. 199-210, DOI:10.1016/s1350-4533(98)00004-6. [78] Maxian, T., Brown, T.D., Pederson, R.D., Callaghan, J.J. (1996). A sliding distance coupled finite element formulation for polyethylene wear in total hip arthroplasty. Journal of Biomechanics, vol. 29, no. 5, p. 687-692, DOI:10.1016/0021­9290(95)00125-5. [79] Chernets, M., Opielak, M., Kornienko, A., Radko, O. (2021). Predictive estimation of sliding bearing load-carrying capacity and tribological durability. Strojniški vestnik - Journal of Mechanical Engineering, vol. 67, no. 7-8, p. 363-368, DOI:10.5545/sv-jme.2021.7139. [80] Mukras, S.M.S. (2020). Computer Simulation/Prediction of Wear in Mechanical Components. Advances in Tribology, vol. 2020, p. 8867351, DOI:10.1155/2020/8867351. [81] Mukras, S., Kim, N.H., Sawyer, W.G., Jackson, D.B., Bergquist, L.W. (2009). Numerical integration schemes and parallel computation for wear prediction using finite element method. Wear, vol. 266, no. 7, p. 822-831, DOI:10.1016/j. wear.2008.12.016. [82] Mukras, S., Kim, N.H., Mauntler, N.A., Schmitz, T.L., Sawyer, W.G. (2010). Analysis of planar multibody systems with revolute joint wear. Wear, vol. 268, no. 5-6, pp. 643-652, DOI:10.1016/j.wear.2009.10.014. [83] Mukras, S., Kim, N.H., Mauntler, N.A., Schmitz, T.L., Sawyer, W.G. (2010). Comparison between elastic foundation and contact force models in wear analysis of planar multibody system. Journal of Tribology, vol. 132, no. 3, p. 031604, DOI:10.1115/1.4001786. [84] Kim, N.H., Won, D., Burris, D., Holtkamp, B., Gessel, G.R., Swanson, P., Sawyer, W.G. (2005). Finite element analysis and experiments of metal/metal wear in oscillatory contacts. Wear, vol. 258, no. 11-12, p. 1787-1793, DOI:10.1016/j. wear.2004.12.014. [85] Uguz, A., Oka, S. (2004). Finite element modeling of ball joints under dynamic loading. Materials Testing, vol. 46, no. 10, p. 506-512, DOI:10.3139/120.100618. [86] Martins, R.C., Moura, P.S. and Seabra, J.O. (2006). MoS2/Ti Low-Friction Coating for Gears. Tribology International, vol. 39, no. 12, p. 1686-1697, DOI:10.1016/j.triboint.2006.02.065. [87] Bruno, G., Fanara, C., Guglielmetti, F. and Malard, B. (2006). Characterization and Residual Stress Analysis of WearResistant Mo Thermal Spray-Coated Steel Gear Wheels. Surface and Coatings Technology, vol. 200, no. 14, p. 4266­4276, DOI:10.1016/j.surfcoat.2005.01.103. [88] Yonekura, D., Chittenden, R.J. and Dearnley, P.A. (2005). Wear mechanisms of steel roller bearings protected by thin, hard and low friction coatings. Wear, vol. 259, no. 1, p. 779-788, DOI:10.1016/j.wear.2004.12.008. [89] Liu, C., Bi, Q., Matthews, A. (2003). Tribological and electrochemical performance of PVD tin coatings on the femoral head of Ti6Al4V articial hip joints. Surface and Coatings Technology, vol. 163-164, no. 7, p. 597-604, DOI:10.1016/S0257-8972(02)00630-8. [90] Diao, D., Kandori, A. (2006). Finite element analysis of the effect of interfacial roughness and adhesion strength onthe local delamination of hard coating under sliding contact. Tribology International, vol. 39, no. 9, p. 849-855, DOI:10.1016/j.triboint.2005.07.037. [91] Kanber, B., Demirhan, N. (2013). Finite element analysis of frictional contact of elastic solids with thin and moderately thick coatings. Turkish Journal of Engineering and Environmental Sciences, vol. 37, no. 2, p. 162-177, DOI:10.3906/muh-1205-14. [92] Syed, Z., Mohamed, H., Mohammed, M.A.H. (2014). Evaluation of contact stress distribution of hip joint model using finite element method. International Journal of Engineering Research & Technology, vol. 3, no. 6, p. 1603­1610. [93] Cilingir, A.C., Ucar, V., Udofia, I.J., Jinet, Z.M. (2008). Biphasic finite element modelling of contact mechanics of hemi­arthoplasty of human hip joint. Part II: Polycarbonate urethane on cartilage contact. Journal of Trends Biomaterial Artificial Organs, vol. 22, no. 2, p. 65-72. [94] Shinde, J., Kadam, S., Patil, A., Pandit, S. (2016). Design modification and analysis of suspension ball joint using finite element analysis. International Journal of Innovative Research in Science, Engineering and Technology, vol. 5, no. 7, p. 12797-12804, DOI:10.15680/IJIRSET.2016.0507135. [95] Shinde, J., Kadam, S. (2016). Design of suspension ball joint using FEA and experimental method. International Research Journal of Engineering and Technology, vol. 3, no. 7, p. 1853­1858. [96] Jang, B.H., Lee, K.H. (2014). Analysis and design of a ball joint, considering manufacturing process. Proceedings of the Institution of Mechanical Engineers, Part C: The Journal of Mechanical Engineering Science, vol. 228, no. 1, p. 146-151, DOI:10.1177/0954406213497317. [97] Hugnell, A., Andersson, S. (1994) Simulating follower wear in a cam-follower contact. Wear, vol. 179, no. 1-2, p. 101-107, DOI:10.1016/0043-1648(94)90226-7. [98] Hugnell, A., Bjrklund, S., Andersson, S. (1996). Simulation of the mild wear in a cam-follower contact with follower rotation. Wear, vol. 199, no. 2, p. 202-210, DOI:10.1016/0043­1648(96)06920-7. [99] Flodin, A., Andersson, S. (1997). Simulation of mild wear in spur gears. Wear, vol. 207, no. 1-2, p. 16-23, DOI:10.1016/ s0043-1648(96)07467-4. [100] Flodin, A., Andersson, S. (1999). Wear simulation of spur gears. Tribotest, vol. 5, no. 3, p. 225-249, DOI:10.1002/ tt.3020050303. [101] Flodin, A. (2000). Wear of Spur and Helical Gears. Royal Institute of Technology, Stockholm. [102] Flodin, A., Andersson, S. (2000). Simulation of mild wear in helical gears. Wear, vol. 241, no. 2, p. 123-128, DOI:10.1016/ S0043-1648(00)00384-7. [103] Pdra, P., Andersson, S. (1997). Wear simulation with the Winkler surface model. Wear, vol. 207, no. 1, p. 79-85, DOI:10.1016/S0043-1648(96)07468-6. [104] Sfantos, G.K., Aliabadi, M.H. (2007). A boundary element formulation for three-dimensional sliding wear simulation. Wear, vol. 262, no. 5, pp. 672-683, DOI:10.1016/j. wear.2006.08.008. [105] Sfantos, G.K., Aliabadi, M.H. (2006). Wear simulation using an incremental sliding boundary element method. Wear, vol. 260, no. 9-10, p. 1119-1128, DOI:10.1016/j. wear.2005.07.020. [106] Rodríguez-Tembleque, L., Abascal, R., Abascal, R., Aliabadi, M.H. (2010). A boundary element formulation for wear modeling on 3d contact and rolling-contact problems. International Journal of Solids and Structures, vol. 47, no. 18­19, p. 2600-2612, DOI:10.1016/j.ijsolstr.2010.05.021. [107] Wang, C., Schipper, D.J. (2020). A numerical-analytical approach to determining the real contact area of rough surface contact. Tribology - Materials, Surfaces & Interfaces, vol. 14, no. 3, p. 166-176, DOI:10.1080/17515831.2020.172 0382. [108] Liew, H.V., Lim, T.C. (2005). Analysis of time-varying rolling element bearing characteristics. Journal of Sound and Vibration, vol. 283, no. 3-5, p. 1163-1179, DOI:10.1016/j. jsv.2004.06.022. [109] Liu, F., Wu, W., Hu, J., Yuan, S. (2019). Design of multi-range hydro-mechanical transmission using modular method. Mechanical Systems and Signal Processing, vol. 126, p. 1-20. DOI:10.1016/j.ymssp.2019.01.061. [110] Liu J., Tang C., Shao Y. (2019). An innovative dynamic model for vibration analysis of a flexible roller bearing. Mechanism and Machine Theory, vol. 135, p. 27-39, DOI:10.1016/j. mechmachtheory.2019.01.027. [111] Ascubal Ascurol. (2021). Rod ends with roller bearing, from https://www.directindustry.com/industrial-manufacturer/ rod-end-roller-bearing-249195.html, accessed on 2021-12­31. [112] Kipp.com (2021). K0716 Rod ends with ball bearing external thread, from https://www.kipp.com/gb/en/Products/ Operating-parts-standard-elements/Joints/K0716-Rod-ends­with-ball-bearing-and-external-thread.html, accessed on 2021-12-31. [113] Liu, J., Tang, C., Wu, H., Xu, Z., Wang, L. (2019). An analytical calculation method of the load distribution and stiffness of an angular contact ball bearing. Mechanism and Machine Theory, vol. 142, p. 103597, DOI:10.1016/j. mechmachtheory.2019.103597. [114] Palmgren, A. (1959). Ball and Roller Bearing Engineering. SKF Industries Inc., Philadelphia, USA. [115] Houpert, L. (1997). A uniform analytical approach for ball and roller bearings calculations. Journal of Tribology, vol. 119, no. 4, p. 851-858, DOI:10.1115/1.2833896. [116] Hernot, X., Sartor, M., Guillot, J. (2000). Calculation of the stiffness matrix of angular contact ball bearings by using the analytical approach. Journal of Mechanical Design, vol. 122, no. 1, p. 83-90, DOI:10.1115/1.533548. [117] Harris, T.A. (2001). Rolling Bearing Analysis. 4th ed., John Wiley & Sons, New York. [118] Spiewak, S.A., Nickel, T. (2001). Vibration based preload estimation in machine tool spindles. International Journal of Machine Tools and Manufacture, vol. 41, no. 4, p. 567-588, DOI:10.1016/S0890-6955(00)00081-X. [119] Hwang, Y.K., Lee, C.M. (2009). Development of automatic variable preload device for spindle bearing by using centrifugal force. International Journal of Machine Tools and Manufacture, vol. 49, no. 10, p. 781-787, DOI:10.1016/j. ijmachtools.2009.04.002. [120] Chen, J.S., Chen, K.W. (2005). Bearing load analysis and control of a motorized high speed spindle. International Journal of Machine Tools and Manufacture, vol. 45, no. 12­13, p. 1487-1493, DOI:10.1016/j.ijmachtools.2005.01.024. [121] Jones, A.B. (1960). A general theory for elastically constrained ball and radial roller bearings under arbitrary load and speed conditions. Journal of Basic Engineering, vol. 82, no. 2, p. 309-320, DOI:10.1115/1.3662587. [122] Li, X., Yu, K., Ma, H., Cao, L., Luo, Z., Li, H., Che, L. (2018). Analysis of varying contact angles and load distributions in defective angular contact ball bearing. Engineering Failure Analysis, vol. 91, p. 449-464, DOI:10.1016/j. engfailanal.2018.04.050. [123] Fang, B., Zhang, J., Hong, J. (2017). Quick calculation method and contact angle analysis for high-speed angular contact ball bearing under combined loads. Journal of Xi’an Jiaotong University, vol. 51, no. 6, p. 115-121, DOI:10.7652/ xjtuxb201706019. [124] Tang, Y.B., Gao, D.P., Luo, G.H. (2006). Research of aero-engine high-speed ball bearing. Journal of Aerospace Power, vol. 21, p. 354-360. [125] Feng, J., Sun, Z., Li, H. (2016). Investigation of the stiffness and displacement of multiple combinations of angular contact ball bearings. Journal of Vibrational Measurement & Diagnosis, vol. 36, no. 4, p. 784-789, DOI:10.16450/j.cnki. issn.1004-6801.2016.04.027. [126] Noel, D., Ritou, M., Furet, B., Le Loch, S. (2013). Complete analytical expression of the stiffness matrix of angular contact ball bearings. Journal of Tribology, vol. 135, no. 4, p. 041101, DOI:10.1115/1.4024109. [127] Sheng, X., Li, B., Wu, Z., Li, H. (2014). Calculation of ball bearing speed-varying stiffness. Mechanism and Machine Theory, vol. 81, p. 166-180, DOI:10.1016/j.mechmachtheory.2014.07.003. [128] Zhang, J., Fang, B., Hong, J., Wan, S., Zhu, Y. (2017). A general model for preload calculation and stiffness analysis for combined angular contact ball bearings. Journal of Sound and Vibration, vol. 411, p. 435-449, DOI:10.1016/j. jsv.2017.09.019. [129] Zhang, J., Fang, B., Zhu, Y., Hong, J. (2017). A comparative study and stiffness analysis of angular contact ball bearings under different preload mechanisms. Mechanism and Machine Theory, vol. 115, p. 1-17, DOI:10.1016/j. mechmachtheory.2017.03.012. [130] Fang, B., Zhang, J., Yan, K., Hong, J., Wang, M.Y. (2019). A comprehensive study on the speed-varying stiffness of ball bearing under different load conditions. Mechanism and Machine Theory, vol. 136, p. 1-13, handle/1783.1/96320. [131] Cheng, Q., Qi, B., Liu, Z., Zhang, C., Xue, D. (2019). An accuracy degradation analysis of ball screw mechanism considering time-varying motion and loading working conditions. Mechanism and Machine Theory, vol. 134, p. 1-23, DOI:10.1016/j.mechmachtheory.2018.12.024. [132] Jedrzejewski, J., Kwasny, W. (2010). Modelling of angular contact ball bearings and axial displacements for high-speed spindles. CIRP Annals, vol. 59, no. 1, p. 377-382, DOI:10.1016/j.cirp.2010.03.026. [133] Zhang, J., Fang, B., Hong, J., Zhu, Y. (2017). Effect of preload on ball-raceway contact state and fatigue life of angular contact ball bearing. Tribology International, vol. 114, p. 365­372, DOI:10.1016/j.triboint.2017.04.029. [134] Ren, X., Zhai, J., Ren, G. (2017). Calculation of radial load distribution on ball and roller bearings with positive, negative and zero clearance. International Journal of Mechanical Sciences, vol. 131 p. 1-7, handle/181551/8818. [135] Guo, Y., Parker, R.G. (2012). Stiffness matrix calculation of rolling element bearings using a finite element/contact mechanics model. Mechanism and Machine Theory, vol. 51, p. 32-45, DOI:10.1016/j.mechmachtheory.2011.12.006. [136] Daidié, A., Chaib, Z., Ghosn, A. (2008). 3D simplified finite elements analysis of load and contact angle in a slewing ball bearing. Journal of Mechanical Design, vol. 130, no. 8, p. 082601, DOI:10.1115/1.2918915. [137] Wang, H., Han, Q., Zhou, D. (2017). Nonlinear dynamic modeling of rotor system supported by angular contact ball bearings. Mechanical Systems and Signal Processing, vol. 85, p. 16-40, DOI:10.1016/j.ymssp.2016.07.049. [138] Ansys Contact Technology Guide. Ansys Release 2020 R1, from https://d.shikey.com/down/Ansys.Products.2020. R1.x64/install_docs/Ansys.Products.PDF.Docs.2020R1/ readme.html, accessed on 2021-12-31. [139] Shinohara, K. (2006). Thermoplastic Resin Composition Containing Mesoporous Powders Absorbed with Lubricating Oils, Patent WO 2006057974 A1. [140] Rabinowicz, E. (1995). Friction and Wear of Materials. Wiley-Interscience. [141] Zhang, B., Quinez, A., Venner, C. (2020). Effect of material anisotropy on rolling contact fatigue life under dry and lubricated point contact conditions: A numerical study. Tribology International, vol. 152, p. 106584, DOI:10.1016/j. triboint.2020.106584. [142] Otero, J., De la Guerra, E., Tanarro, E.C., Río, B. (2017). Friction coefficient in mixed lubrication: A simplified analytical approach for highly loaded non-conformal contacts. Advances in Mechanical Engineering, vol. 9, no. 7, DOI:10.1177/1687814017706266. [143] Zhu, D., Hu, Y. (2001). A computer program package for the prediction of EHL and mixed lubrication characteristics, friction, subsurface stresses, and flash temperatures based on measured 3-D surface roughness. Tribology Transactions, vol. 44, p. 383-390, DOI:10.1080/10402000108982471. [144] Castro, J., Seabra, J. (2007) Coefficient of friction in mixed film lubrication: gears versus twin-discs. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 221, no. 3, p. 399-411, DOI:10.1243/13506501JET257. [145] van Beek, A. (2012). Advanced engineering design: Lifetime performance and reliability. TU Delft, The Netherlands. [146] Thompson, M.K. (2007). A multi-scale iterative approach for finite element modeling of thermal contact resistance. PhD Thesis. Massachusetts Institute of Technology, Massachusetts, DOI:10.1115/MNHT2008-52385. [147] Nayak, P.R. (1971). Random process model of rough surfaces. Journal of Lubrication Technology, vol. 93, no. 3, p. 398-407, DOI:10.1115/1.3451608. [148] Kartini, Saputra, E., Ismail, R., Jamari, J., Bayuseno, A.P. (2016). Analysis of the contact area of smooth and rough surfaces in contact with sphere indenter using finite element method. MATEC Web of Conferences, vol. 58, p. 04007, DOI:10.1051/matecconf/20165804007. [149] Wu, J.-J. (2000). Simulation of rough surfaces with FFT. Tribology International, vol. 33, no. 1, p. 47-58, DOI:10.1016/ S0301-679X(00)00016-5. [150] Hu, Y.Z., Tonder, K. (1992). Simulation of 3-D random surface by 2-D digital filter and Fourier analysis. International Journal of Machine Tools and Manufacture, vol. 32, no. 1-2, p. 82-90. [151] Perez-Rafols, F., Almqvist, A. (2019). Generating randomly rough surfaces with given height probability distribution and power spectrum. Tribology International, vol. 131, p. 591-604, DOI:10.1016/j.triboint.2018.11.020. [152] Pawlus, P.; Reizer, R.; Wieczorowski, M. (2020). A review of methods of random surface topography modeling. Tribology International, vol. 152, p. 106530, DOI:10.1016/j. triboint.2020.106530. [153] You, S.J., Ehmann, K.F. (1991). Computer synthesis of three- dimensional surfaces. Wear, vol. 145, no. 1, p. 29-42, DOI:10.1016/0043-1648(91)90237-O. [154] Chang, W.R., Etsion, I., Bogy, D.B. (1987). An elasticplastic model for the contact of rough surfaces. Journal of Tribology, vol. 109, p. 257-263. [155] Zhao, Y., Maietta, D.M., Chang, L. (2000). An asperity microcontact model incorporating the transition from elastic deformation to fully plastic flow. Journal of Tribology, vol. 122, no. 1, p. 86-93, DOI:10.1115/1.555332. [156] Kogut, L., Etsion, I. (2002). Elastic-plastic contact analysis of a sphere and a rigid flat. Journal of Applied Mechanics, vol. 69, no. 5, p. 657-662, DOI:10.1115/1.1490373. [157] Jackson, R.L., Green, I. (2005). A finite element study of elasto-plastic hemispherical contact against a rigid flat. Journal of Tribology, vol. 127, no. 2, p. 343-354, DOI:10.1115/1.1866166. [158] Mishra, T., de Rooij, M., Schipper, D.J. (2021). The effect of asperity geometry on the wear behaviour in sliding of an elliptical asperity. Wear, vol. 470-471, p. 203615, DOI:10.1016/j.wear.2021.203615. [159] Zhou, C., Wang, H., Wang, H., Hu, B. (2021). Three-dimensional asperity model of rough surfaces based on valley-peak ratio of the maximum peak. Meccanica, vol. 56, p. 711-730, DOI:10.1007/s11012-021-01309-3. [160] Zhu, L., Chen, J., Sun, Y. (2021). On the approximating criteria of parabolic asperities for measured surface profiles. American Institute of Physics Advances, vol. 11, p. 035134, DOI:10.1063/5.0046579. [161] Tsao, Y.H., Tong, K.N. (1975). A model for mixed lubrication. American Society of Lubrication Engineers Transactions, vol. 18, no. 2, p. 90-96, DOI:1080/05698197508982750. [162] Szczerek M., Wisniewski M. (2001). Theories of mixed friction. Szczerek M., Wisniewski M. (eds). Tribolgy and Tribotechnics. ITeE. Radom. [163] Polyacetal (POM) Plastem, (2021). from http://www.plastem. pl/en/offer/plastics/polyacetal-pom/, accessed on 2021-12­31. [164] Bowden, F.P., Tabor, D. (1950) & (1964). The Friction and Lubrication of Solids, Part I (1950) & Part II (1964). Oxford Clarendon Press, Oxford. [165] Patir, N., Cheng, H.S. (1978). An average flow model for determining effects of three-dimensional roughness on partial hydrodynamic lubrication. Journal of Lubrication Technology, vol. 100, no. 1, p. 12-17, DOI:10.1115/1.3453103. [166] Wilson, W., Marsault, N. (1998). Partial hydrodynamic lubrication with large fractional contact areas. Journal of Tribology, vol. 120, no. 1, p. 16-20, DOI:10.1115/1.2834180. [167] Kakoi, K. (2021). Formulation to calculate isothermal, non-newtonian elastohydrodynamic lubrication problems using a pressure gradient coordinate system, and its verification by an experimental grease. Lubricants, vol. 9, no. 5, p. 56, DOI:10.3390/lubricants9050056. [168] Radulescu, A.V., Radulescu, I. (2006). Rheological models for lithium and calcium greases. Mechanica, vol. 59, no. 3, p. 67­70. [169] Clarke, S., Bell, M.A. (2008). Modelling of friction in a ball joint in mixed lubrication. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 222, p. 267-275, DOI:10.1243/09544062JMES437. [170] Walicki, E., Walicka, A. (2000). Mathematical modelling of some biological bearings. Smart Materials and Structures, vol. 9, no. 3, p. 280-283, DOI:10.1088/0964-1726/9/3/305. [171] Pacejka, H.B., Bakker, E. (1993). The Magic Formula Tyre Model. Pacejka H.B. (eds). Tyre Models For Vehicle Dynamics Analysis, 1st ed., CRC Press, p. 1-18, DOI:10.1080/00423119208969994. [172] Rega, G., Wang, X., Shi, S., (2015). Analysis of vehicle steering and driving bifurcation characteristics. Mathematical Problems in Engineering, vol. 2015, p. 847258, DOI:10.1155/2015/847258. [173] Nguyen, V. (2005). Vehicle handling, stability, and bifurcation analysis for nonlinear vehicle models. PhD thesis. College Park, University of Maryland, Maryland. [174] Sahin, Y., Pauw, J., Sukumaran, J., de Baets, P. (2015). Sliding Friction and Wear of Polyoxymethylene Polymer. Synergy International Conferences - Engineering, Agriculture and Green Industry Innovation. [175] Mergler, Y.J., Schaake, R.P., Huis in’t Veld, A.J. (2004). Material transfer of POM in sliding contact. Wear, vol. 256, p. 294-301, DOI:10.1016/S0043-1648(03)00410-1. [176] Zunda, A., Padgurskas, J., Jankauskas, V., Levinskas, R. (2010). Wear Resistance of Industrial Polymers Under Lubrication with Oils. Scientific Journal of Riga Technical University. Material Science and Applied Chemistry, vol. 21, p. 21-25. [177] Golchin, A., Simmons, G., Glavatskih, S., Prakash, B. (2013). Tribological behaviour of polymeric materials in water-lubricated contacts. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 227, no. 8, p. 811-825, DOI:10.1177/1350650113476441. Strojniški vestnik - Journal of Mechanical Engineering 68(2022)2, 126-140 Received for review: 2021-10-25 © 2022 The Authors. CC BY-NC 4.0 Int. Licensee: SV-JME Received revised form: 2021-12-30 DOI:10.5545/sv-jme.2021.7452 Original Scientific Paper Accepted for publication: 2022-01-04 Simulation and Exp erimental R esearch on Electrical Control Anti-backlash Based on a Novel Typ e of Variable Tooth Thickness Involute Gear Pair Guangjian Wang1,2 – D ongdong Z hu1 ,2 – S huaidong Z ou1,2,* – Yujiang Jiang1,2 – X in Tian1,2 1 Chongqing University, State Key Laboratory of Mechanical Transmissions, China 2 Chongqing University, College of Mechanical and Vehicle Engineering, China This paper presents an electrical backlash control method based on a novel variable tooth thickness involute gear pair (VTTIGP) through acceleration judgment. First, the principle, characteristics and theoretical tooth surface equation of the novel type of variable tooth thickness involute gear (VTTIG) are proposed. Compared with the traditional variable tooth thickness gear (VTTG), each transverse of the novel type of VTTIG has the same tooth root circle and the addendum circle. Therefore, there will be no interference caused by the excessive axial relative movement distance between the gear pairs when eliminating the backlash. The theoretical calculation equation of no-load transmission error (TE) and time-varying backlash for a dual-eccentric novel type of VTTIGP are deduced. Second, using the theoretical backlash equation and the eccentricities and the initial phases obtained through the experimental TE curve obtained with angular displacement sensors and the developed TE equation, the continuous backlash curve of the gear pair is obtained. Based on the RecurDyn software and its Colink control module, the co-simulation model of the novel type of VTTIGP and its backlash control system is used to explore the electrical backlash control method through acceleration judgment in the case of time-varying backlash. Finally, using angular displacement sensors and linear actuators, the experimental study on the anti-backlash control with different reversed positions under different input speeds is carried out, demonstrating that this novel electrical anti-backlash control method is effective. Keywords: anti-backlash, electrical control, time-varying backlash, transmission error, variable tooth thickness gear pair Highlights • An electrical backlash control method based on the new type of VTTIGP through acceleration judgment is proposed. • The principle, characteristics, and theoretical tooth surface equation of the novel type of VTTIG are proposed. • The formula of the non-load TE and time-varying backlash based on the eccentricity error model of the novel type of VTTIG is derived. • The simulation and experiment further verify the effectiveness of the electrical anti-backlash control method based on acceleration judgment. 0 INTRODUCTION Gear transmission is widely used in precision servo transmission systems, such as robots, machine tools and tracking and positioning systems, which require fast responses, high positioning accuracy, and good stability [1] and [2]. However, gears will cause errors in the process of machining and assembly, and these errors often affect the backlash of the gear transmission process [1] and [3]. Gear backlash is a nonlinear factor affecting gear performance [4]. For a reciprocating system, the non-linearity of the gear tooth backlash will directly affect the stability, accuracy, and sensitivity of the servo control, resulting in delay and discontinuity of the output motion response [5] to [8]. Therefore, adjusting and controlling the gear backlash is of great significance to improve the accuracy of the transmission system. In practical applications, it is necessary to adopt special anti-backlash mechanisms and technical measures to eliminate backlash when the machining and assembly accuracy reaches a certain level [9] and [10]. R esearch on gear anti-backlash has always been a topic of great interest, and various methods of anti-backlash exist. Generally speaking, the commonly used anti-backlash methods include adjustable centre distance, spring-loaded gears, fixed split gears, and adjustable tooth thickness gear [10] and [11]. The variable tooth thickness gear, also known as the beveloid gear, which was first proposed by Beam in 1954 [12], has become an important gear drive for power and motion transmission due to its advantages, such as abilities to backlash control, misalignments absorption, and easy manufacturing. It is widely used for precision transmission, such as rotate vector (R V) reducer, microscope gear, and others [13] and [14]. In recent years, people have conducted extensive research on the geometric design, meshing characteristics and applications of VTTG. Do et al. [15] modelled and simulated the variable-thickness gear transmission. Based on the hob-machining theory, Yang et al. [16] established a mathematical model of the tooth surface *Corr. Author’s Address: Chongqing University, State Key Laboratory of Mechanical Transmissions, Chongqing, 400044, China, zoushuaidong@cqu.edu.cn of a VTTG considering the machine tool adjustment errors. Cao et al. [17] proposed a plunge-shaving model for variable tooth thickness based on ordinary gear-shaving machine tools’ structure. Brauer [18] used a global–l ocal FE approach to study the transmission error of VTTG. Song et al. studied the meshing stiffness of VTTG pairs [14], [19] and [20]. Chen et al. [13] studied the contact characteristics of variable-thickness internal gear transmission and the machining principle of internal variable-thickness gears. Li et al. [21] used the VTTIG as the transmission part of the robot R V reducer. Yan et al. [22] used the VTTIG as the transmission part of the 3K reducer. Z heng et al. [23] and Ye et al. [24] combined VTTIG with a double-lead involute cylindrical worm and proposed a new type of worm gear with adjustable backlash. However, the above-mentioned VTTG realiz es the axial change of the tooth thickness by adopting different radial modification coefficients, which makes the tooth root circle and the addendum circle of each transverse different. Therefore, interference is easy to occur during the anti-backlash process of the axial movement of the gear pair. Among the commonly used anti-backlash methods, adjustable centre distance, fixed split gears, and adjustable tooth thickness gear only eliminate constant backlash, while spring-loaded gears are limited to light load applications. With the development of control technology, it has become increasingly common to eliminate the backlash through algorithms control or the combination of mechanical structures and algorithms control [25] to [28]. Shi and Z uo [29] proposed a “ soft degree” concept and presented a backstepping algorithm to overcome the backlash non-linearity in gear transmission servo systems. For low operating speed systems with relative cyclic motion, Warnecke and Jouaneh [7] presented an approach for backlash compensation via modification of the input motion. However, this method has a high requirement for motor response, and it also has low anti-interference capability. Sato et al. [30] adopted the piez oelectric actuators to create a preload in a lead screw and nut mechanism to control and eliminate the backlash. However, it is limited to low torque and force applications. Presently, in robots, tracking and positioning systems, and computer numerical control (CNC) machine tools, the backlash is adjusted and controlled by the method of dual-motor driven electrical anti-backlash based on bias torque or speed difference [1] and [31] to [33]. However, because of the need for two sets of transmission chains, two sets of servo systems, and a backlash control card, the volume and the cost of the system are high. In a precision transmission system, eccentricity is the main source of large-period time-varying TE and periodic time-varying backlash [1]. If the total backlash law curve, including constant backlash and periodic time-varying backlash, can be determined, the backlash can be controlled and adjusted according to the characteristics that the novel type of VTTIGP can adjust the backlash through relative axial movement. Based on this, this paper presents a method of electrical control anti-backlash based on the novel type of VTTIGP, which has a much relative smaller volume and lower cost: based on determining the continuous backlash curve of the novel type of VTTIGP, according to the acceleration change of the novel type of VTTIGP, and by adjusting the real-time position relationship between the driving gear and the driven gear of the novel type of VTTIGP, the backlash can be controlled. The structure of this paper is as follows. First, the principle, characteristics and theoretical tooth surface equation of the novel type of variable tooth thickness gear are proposed. Then, calculation models with a theoretical no-load TE and backlash of the novel type of VTTIGP with dual-eccentricities were constructed. Next, the eccentricity and its initial phase of the experimental device were obtained by fitting the theoretical no-load TE equation and then the backlash curve of the gear pair was acquired through the backlash equation and the resultant eccentricity. Finally, according to the acquired backlash curve, and through simulation and experiment, the electrical anti-backlash control method based on acceleration judgment was verified. 1 METHODS 1.1 The Novel Type of VTTIG Based on Compound Modification 1.1.1 The Principle and Characteristics of A New Type of VTTIG As shown in Figs. 1 and 2, the two tooth surfaces of the novel type of VTTIG form a certain cone angle along the gear axis. By adjusting the relative position of the meshing gears along the axis, the meshing backlash can be adjusted to realiz e a slight backlash or no backlash transmission. In order to facilitate the meshing backlash control of the gear pair, the novel type of VTTIG adopts a constant radial modification and a gradual tangential modification at the same time, and the tangential modification of the left and right tooth surfaces changes linearly at the same slope. Compared with the traditional VTTG, each transverse of the novel type of VTTIG has the same tooth root circle and the addendum circle. Therefore, no interference will be caused by the excessive axial relative movement distance between the gear pairs while eliminating the backlash. Fig. 2. The 3D model of the novel type of VTTIGP As shown in Fig. 3a , the indexing tooth thickness of the novel type of VTTIG consists of standard tooth thickness (1) , radial modification coefficient increments (2), and gradually changing tangential modification coefficient increment (3 ). As shown in Fig. 3b, the tooth thickness of the intermediate transverse is p, and through tangential modification, the tooth thickness of the other transverse is p plus or minus the tangential modification coefficient increment x . When the tangential modification coefficient changes linearly along the axial direction, the generated gear tooth surface is an involute helical surface, and the left and right tooth surfaces of the gear teeth are double helical surfaces with different helix angles or directions. At the same time, it can also be modified according to specific regular changes to form drum-shaped teeth and other teeth shapes, thereby improving the load-bearing capacity and reducing vibration, impact and noise during transmission. However, it is the most convenient to use the linearly changing tangential modification coefficient in the actual application process. Fig. 3. The comparison of indexing circle expansion and end face; a) index circle expansion diagram, and b) tangentially modified tooth profile 1.1.2 The Tooth Surface Equation of the New Type of VTTIG Litvin [34] proposed a method of generating helical gears with a pair of conjugate female surfaces. This method can be applied to the derivation of the tooth surface equation and the normal vector of any point on the tooth surface [35]. Assuming that the gear is processed by the rack cutter through the method of generating an envelope, the rigid body motion relationship is mainly considered [36], and the influence of elastic deformation is neglected. The parameters of the rack cutter are as follows, the left and right helix angles of the rack are ß and ß, RL and the left and right helix angles of the gear can be obtained as ßL and ßR ; normal pressure angle a= 20° , n transverse pressure angle: AC transverse pressure angle is aat1 L , DF transverse pressure angle is aat1, and their siz es may be different; the modulus is m, the bottom clearance cp = 0.25 m, the addendum hap = m, and dedendum hf p = 1.25 m. As shown in Fig. 4, the blue contour is the rack cutter of a novel type of VTTIG with compound modification. The helix angle of the orange part is ß. After tangential modification, ß is R RL obtained by keeping the transverse pressure angle of the indexing circle unchanged. According to the basic theorem of the Willis principle, the gear tooth profile curve and the gear cutter curve are conjugate tooth shapes. At any instant of transmission, their common normal line at the contact point must pass through the instantaneous node. Based on this, the discrete equation of the contact point K can be obtained. Coordinate systems [36] are established, as illustrated in Fig. 5 . S 1(O 1 x 1 y1 1 ), S 2(O 2x 2y2 2) and S 0(O 0 –x 0y0 0) are the coordinate system fixedly connected to the rack, the fixedly connected gear rotating coordinate system and auxiliary static coordinate system, respectively. r2 and f are the reference pitch circle radius and the rotation angle for the processing, respectively. In the S 1 coordinate system, the coordinate equation R 1 AC of the point K on the line segment AC on any transverse can be expressed as follows: R AC dl ˜xyz . (,) ° 1 111 d sin ˆ.l tan ...m /4 .h tan ˆT at L 1 L ap 1 at L  ˜d coos ˆat L .hap, (1) 1  l  R DF (,) dl ˜°xyz . 1 111 .d sin ˆ.l tan ...m /4 .h tan ˆT at R 1 R ap 1 at R  ˜d ccos ˆat R .hap, (2) 1  l  where d is the distance from any point K to point F on the line segment DF on any transverse, and l is the Z coordinate value of the studied transverse. When the rack pitch plane and the gear pitch circle roll with each other, the gear tooth surface is processed using the rack tooth surface envelope. By deriving the variables in the coordinate equation R 1 AC of any point on the line AC segment respectively, the following equations can be obtained as: . sin ˆ . tan  ˆ ˜ ˜˜ ˜ RR cos dl Thus, the normal vector of the tooth surface of the AC segment of the basic rack can be obtained as: 1 at L L AC AC .... .... .... . ... (3) ° ° 1 1 0 , . 1 at L 0 1 . ˆ cos  at L 1 .... .... . (4) nAC ˜ ° sin at 1 L at L 1   ° tan cos L Similarly, the tooth surface equation and normal vector of the other parts of the rack can be obtained. According to the basic principle of gear meshing, Fig. 4. The parameters of rack cutter the condition for any point K (x 1 i, y1 j, 1 k) on the tooth surface to become a contact point is that the normal vector direction is perpendicular to the sliding speed direction [35]. If the sliding speed of the gear contact point K relative to the rack is represented by v121 , the speed of the rack in the S 1 coordinate system is v11 , and the speed of the gear at point K in the S 1 coordinate system is v2 , the following equation can be obtained 1 as: ˜ ..° v 11 v 12 . ˆˆ. . ˜.°  . ri 2 OK 2 .ˆ. . (5) In the S 1 coordinate system, the angular velocity vector and vector O 2K can be expressed as follows: where d is the distance from any point K to point A on ˜˜°.˜k , (6) the line segment AC on any transverse, and l is the Z OK OO °OO OK ˜° coordinate value of the studied transverse. 2 20 011 ri ) °. (7) In the S 1 coordinate system, the coordinate ˜( x 1 ..2) °( y 1 °rj zk 21 equation R DF of the point K on the line segment DF 1 According to Eqs. (5) to (7) , the equation of on any transverse can be expressed as follows: relative sliding speed can be obtained as: v12 ˜.(( x °.) °(8) rjyi ). 1 121 The relationship between the normal vector and relative sliding speed can be represented as follows: n AC ˜v 112 °0. (9) According to Eqs. (4) , (8) , and (9) , the following equation can be obtained as: y 1 ˜°x tan ...r tan .1. (10) 1 at 1 L 2 atL According to Eqs. (1) and (10) , the general meshing equation of rack AC segment and gear tooth surface can be represented as follows: tan ˜( d sin ˜.l tan °ˆ.m /4 ˆh tan ˜1) at L 1 atL 1 L ap atL .d cos ˜at L ˆ.r 2ta1 ˆhap.0. (1 1) an 1 ˜at L The conversion relationship between the contact point K in the gear coordinate system and the rack coordinate system can be determined by the following equation [35]: ˜T , ,, °T xyzt ,,, °.M ˜xyzt , (12) 2 22 21111 where the matric M21 is the transformation matrix as: .cos .sin .0 r (sin °cos) . ....2  °sin .cos .0 r 2 (cos .sin) ... M21 ˜.. (13) .001 0  . 000 1 ˆ After finishing Eq. (1 1) , the expression of variable d can be obtained as follows: d ˜sin .( .m /4 °h tan .1 at 1 L ap atL hap °°.rl ˆ 2 .tan L ). (14) tan .at L 1 According to Eqs. (1) , (12) , and (13) , the coordinate conversion relationship can be obtained as follows: (sin .... .x ˜x cos .°y sin .°r cos) 21 12 . y ˜.x sin .°y cos .°r (cos .°.s.(15) ˆsin). 21 12 . z ˜l .2 According to Eqs. (1) , (14) , and (15) , the final tooth surface equation without the intermediate parameter d can be obtained. Similarly, the tooth surface equations of other parts of the gear surface can be obtained. 1.2 The Principle and Strategy of Backlash Control for the Novel Type of VTTIGP The anti-backlash principle for the novel type of VTTIGP is based on the change of tooth thickness along the axial direction of gear tooth thickness. As shown in Fig. 6, backlash is controlled by adjusting the relative axial displacement between the driving and the driven gear. Fig. 6. Anti-backlash principle of the composite modified VTTIGP The composite modified VTTG with left and right tooth profiles of identical tangential modification is introduced in [37], and the relationship between the axial displacement of the driven gear and the backlash can be expressed as follows: j ˜2tans °, (16 ) where j is the backlash of the gear pair, s is the axial displacement of driven gear to eliminate backlash, ß is the helix angle of the novel type of VTTIG. Therefore, the axial displacement of the driven gear required to eliminate backlash can be obtained with the known backlash. Fig. 7 shows a flow chart of backlash adjustment based on acceleration judgment. When the gear pair is detected as beginning to decelerate, the driven gear will be driven by a linear actuator to move axially based on Eq. (16) , which eliminates the backlash when the rotational direction of the gear pair reverses. Subsequently, when the gear pair is detected speeding up, the driven gear will return to the initial axial displacement to prepare for the next backlash-eliminating sequence. 1.3 Calculation of Backlash of Gear Pair with Dual-eccentricities Eccentricity error is the main factor leading to large-period time-varying TE and periodic time-varying backlash. Therefore, this paper first theoretically derives the TE and time-varying backlash equations of the novel type of VTTIGP under eccentricity Fig. 7. Backlash control mode based on acceleration judgment error. Then through the experimental TE curve and the derived TE equation, the eccentricity error of the novel type of VTTIGP is solved in reverse, and then the time-varying backlash curve under the gear pair eccentricity error is obtained. Finally, through the above-mentioned anti-backlash strategy, the active anti-backlash of the novel type of VTTIGP is realiz ed. 1.3.1 Theoretical No-Load TE Calculation Since each transverse of the novel type of VTTIG is an involute gear pair, the novel type of VTTIGP can be analysed according to the relevant theories of involute cylindrical gear pair. If the driving gear has an eccentricity, the meshing point displacement along the line of action will change after rotating each unit angle. Compared with the ideal gear transmission, when the driving gear rotates counter-clockwise or clockwise, as shown in Fig. 8, the displacement increment of the meshing point along the line of action can be expressed as follows [1] and [38]: 1 °) . ˜F .°sin(  11 cc 1 .2 ..e 1 .., (17) ˜F 1 .ˆ.sin( 1 .c ) .1 ˆ where .F 1 1 and .F 12 represent the displacement increment along the line of action caused by eccentricity when the driving gear rotates counter­clockwise and clockwise, respectively; f1 cc , f1 c , e1, and .1 are the angular displacement along the counter­clockwise direction, the angular displacement along the clockwise direction, the eccentricity, and initial phase of eccentricity of driving gear, respectively. Similarly, if the driven gear has an eccentricity and when the driving rotates counter-clockwise or clockwise, as shown in Fig. 9, the displacement increment along the line of action can be calculated as follows: °˜F 1 .°sin(  ) . 22 cc 2 .2 ..e 2 .., (18) .˜F .sin( ) .2 ˆ.2 c 2 ˆ where .F 21 and .F 22 represent the displacement increments along the line of action caused by the eccentricity when the driving gear rotates counter­clockwise and clockwise, respectively; f2, f2, e2, cc c and .2 are the angular displacement along the counter­clockwise direction, the angular displacement along the clockwise direction, the eccentricity, and initial phase of the driven gear, respectively. Fig. 8. Displacement increment along the line of action calculation model with the eccentricity of the driving gear at any transverse section; a) driving gear rotates counter-clockwise, and b) driving gear rotates clockwise Fig. 9. Displacement increment along the line of action calculation model with the eccentricity of driven gear at any transverse section; a) driving gear rotates counter-clockwise, and b) driving gear rotates clockwise To determine that TE is caused only by the eccentricity of no-load TE, the no-load TE of counter­clockwise and clockwise can be acquired when the displacement increment along the line of action translates into the rotational angle of the driven gear. 180 e .. ˜°ˆ.[sin( °..) cc 11 cc 1 .R 2 cos . °...) e sin( ..) .e sin( . 22 cc 2 11 ..(19) e 2 sin( 2 )], ˜°ˆ180 .e . [sin( °..) c 11 c 1 .R 2 cos . °..) e sin( ..) e sin( .... 22 c 2 11 .e sin( .. ..)], (20) 22 here f is no-load TE of counter-clockwise cc rotation, f is no-load TE of clockwise rotation, and c the unit is degree, and R 2 is the reference radius of driven gear. 1.3.2 Theoretical Backlash Calculation The total backlash [39] of the gear pair can be divided into two components: the constant backlash (which is caused by factors including tooth thickness tolerance and allowance, and centre distance tolerance and allowance) and the variable backlash that is due mainly to the eccentricity. Therefore, the variable backlash can be calculated as follows [40]: 180 °2 tan . j ˜°v ˆR 2 [cos( e ...) .e cos( ...)], (21 ) 22 cc 21 1 cc 1 and, following that, the total backlash can be expressed as follows, j ˜°jv .jc , (22) where jf, jv, jv are total backlash, variable backlash, and constant backlash respectively with the units in degrees. Consequently, both no-load TE and backlash can be obtained via the eccentricity error and its initial phase. Therefore, it is first based on the experimental TE curve in order to obtain the eccentricity and its initial phase; then based on Eq. (21) , the acquired eccentricity and initial phase, with the continuous backlash being obtained subsequently. 2 EXPERIMENTAL AND SIMULATION 2.1 Measurement of No-Load Torque TE and Backlash 2.1.1 The Test Rig The experiment used the incremental encoder ER N180 made by Heidenhain in Germany; the system accuracy is 18 arc seconds, the line number is 3600, to output the V signal. Through the subdivision PP circuit matched with the encoder, the signal can be subdivided 20, 25, and 50 times and output TTL signal. In this experiment, 25 times the subdivision is selected, the output pulses per second is 90,000, and the resolution is 14.4 a rc seconds. The experimental test rig (Fig. 10) comprises a novel type of VTTIGP whose parameters are shown in Table 1. The driving gear is fixed on the primary shaft, which is moved by a servo motor. The driven gear has been splined to its shaft and can be moved axially by a linear motor so that the value of the gear backlash can be modified by adjusting the relative position of the driven gear. 1 ˜°°°, (23) .. 21 2 where f1, f2 are the angular displacement of driving gear and driven gear, respectively. To achieve this, two incremental encoders were used as angular displacement sensors, which were fixed to the end of the two shafts, as shown in Fig. 10. Table 1. Parameters of the novel type of VTTIGP Parameters Driving gear Driven gear Number of teeth 45 80 Tip diameter [mm] 94 164 Root diameter [mm] 85 155 Module [mm] 2 2 Tooth width [mm] 20 20 Pressure angle [°] 20 20 Helical angle [°] 3 3 Centre distance [mm] 125.12 Tooth thickness of the thinnest end [mm] 2.093 2.093 Tooth thickness of the thickest end [mm] 4.195 4.195 Transmission ratio 1.778 2.1.2 Measurement of TE The experimental TE curve is acquired when the driving gear rotates counter-clockwise at 60 rpm at an angular displacement of 7560 degrees, then stops two seconds and finally backs to the initial position as shown in Fig. 1 1a . Fig. 1 1b shows the fitting TE curve based on the equation of no-load TE derived above and the experimental TE curve; the eccentricities and the initial phases of the novel type of VTTIGP can be obtained by the fitting method [38]. The acquired eccentricities and the initial phases of the driving and driven gear can be seen in Table 2. The obtained initial phases and the rotation displacements were recorded by the two recorders at the next start of the initial phases of the driving and driven gear so that it is not required to seek the initial phase of the gear pair again at each of the novel type of VTTIGP transmission start. Table 2. Eccentricity and its initial phase of the novel type of VTTIGP Parameters Driving gear Driven gear Eccentricity [mm] 0.0226 0.0273 Initial phase of eccentricity [°] 49.5 79.3 Despite the case of gear eccentric error values and phase determination, TE is a periodic function of the input gear corner and the output gear corner. However, when the input gear rotates at a uniform speed, ignoring the slight influence of the TE, the output gear can also be considered as rotating at a uniform speed, and the TE can be expressed as a function of gear speed and time. If the frequency spectrum analysis is performed on a TE curve obtained under a uniform rotation, the amplitudes of the input frequency and output frequency respectively reflect the eccentricity error of the two gears. In reference [41], it is based on this method to identify the gear eccentricity error. the figure that the amplitudes corresponding to the rotation frequencies of the driving and driven gears are more consistent with the eccentricity values of the two gears obtained by the previous fitting, which proves the correctness of the eccentricity error obtained by the fitting. a) b) c) Fig. 11. TE measurement; a) experimental TE curve, b) theoretical TE equation fitting curve, and c) experimental transmission error FFT diagram Fig. 1 1c is an FFT diagram of the TE in the uniform speed section intercepted when the counter­clockwise input speed is 60 rpm. It can be seen from a) b) Fig. 12. Experimental TE at 60 rpm; a) TE curve of driving gear rotating counter-clockwise, and b) TE curve of driving gear rotating clockwise Since the two gears have 45 and 80 teeth, the driving gear’ s rotation of 5760 degrees is an angular period of the TE. Figs. 12a and b are the experimental TE curves when the driving gear rotates counter­clockwise and clockwise at 60 rpm, respectively. It can be seen from the figures that when the driving gear rotates 5760 degrees, the new angular period of the TE starts, which proves the fundamental periodicity of Eqs. (19) and (20). 2.1.3 Calculation and Measurement of Backlash The theoretical variable backlash of each position of the gear pair is obtained (as shown in Fig. 13) using Eq. (21) and the obtained eccentricities and the initial phases (as shown in Table 1) . Additionally, the constant backlash of gear pair can be obtained using Eq. (22), which is based on the total backlash value when the driving wheel is reversed at the position of 7560 degrees as shown in Fig. 1 1a and the variable backlash at the position of 7560 degrees as calculated by Eq. (21) . Therefore, the theoretical total backlash curve of the gear can be obtained by adding the constant backlash curve and the variable backlash curve. In order to verify the obtained theoretical total backlash curve, the mechanical return backlash method (MR BM), which is obtained by fixing the driving gear at a specified position and then rotating the driven gear clockwise and counter-clockwise (as shown in Fig. 14) is used to calculate the angle difference that is equal to the backlash of the specified position of the driving gear. The results show that the backlash curve obtained by the theoretical equation agrees with the backlash value by MR BM. In summary, the axial displacement to eliminate backlash can be calculated when the gear pair has a reverse motion based on the acquired total backlash and Eq. (16) . 2.2 Backlash Control of Simulation Verification The mechanical and control co-simulation model of the novel type of VTTIGP was established in R ecurDyn (as shown in Fig. 15) , whose gear parameters are consistent with those in Table 1. Simultaneously, a revolute joint is constructed at the eccentric position of the driving gear, and a rotating and axial moving constraint is set up at the eccentric position of the driven gear. The solid contact of the gear pair is also established. Moreover, the input speed of the driving gear is set to 100 rpm, and the driving gear will stop at 7524 degrees and then back to the initial position. Additionally, the external loads include the friction of the revolution joints and the inertia force of the gear pair. Finally, according to the control mode shown in Fig. 7, a simulation control program was established in Colink, as shown in Fig. 15. In order to better show the effect of backlash, the gear pair will stop some seconds when the driving gear is reversed. Fig. 16 shows the results of the simulation TE curve with and without backlash control. It can be seen that the TE curve without backlash control has a jump when the driving gear is reversed, and the jump value is equal to the backlash at the reverse position. However, in the TE curve with backlash control, the a) b) jump phenomenon is not obvious, which demonstrates that the backlash has been controlled effectively. 3 RESULTS AND DISCUSSION 3.1 Comparative Study on Backlash Control through Experimental and Simulation Method 3.1.1 Experimental Anti-Backlash Method As Fig. 7 shows, during gear transmission, the backlash controller receives the angular displacement signal of the encoder of the driving gear in real time and determines whether the speed is reduced or not. When the deceleration signal is detected, according to the current angular displacement signal, the controller calculates the corresponding backlash value and the axial displacement of driven gear required for eliminating backlash. The linear actuator will then receive a signal from the controller to move the driven gear to the specified axial position for anti-backlash. Finally, when the controller detects the acceleration signal of the driving gear, the linear actuator returns to its original axial position. 3.1.2 Backlash Control Comparison and Validation Due to the output part in the application needed to eliminate backlash always having a relatively low speed, the speed of the experiment in this paper is carried out at a relatively low level. Figs. 17 to 20 show the simulation and experimental TE curves of the novel type of VTTIGP with and without backlash control at different reverse positions under different input speeds. According to the encoder signal, the driving gear will return to the initial rotational position after each experimental measurement and then rotate a large period cycle (a large period cycle is defined as that where the driving gear rotates 576 0 degrees), which ensures the initial phases of both the driving and driven gear are the same as the parameters in Table 1. In the figures below, the simulation and experimental TE curves agree well with each other demonstrating that the co-simulation model of the novel type of VTTIGP is consistent with the experimental rig. Due to the eccentricity of the gear pair, the backlash values will not be the same at the different reverse positions when the novel type of VTTIGP is without backlash control, so that the jump values in the TE curve vary, as shown in Figs. 17 a, 18a , 19a , and 20a. However, after backlash control, the jump phenomenon in the TE curve improves in each experiment (as shown in Figs. 17b, 18b, 19b, and 20b), which verifies the validation of the eliminating backlash device and the anti-backlash control mode. a) without backlash control, and b) with backlash control a) without backlash control, and b) with backlash control Wang, G. – Zhu, D. – Zou, S. – Jiang, Y. – Tian, X. a) without backlash control, and b) with backlash control a) without backlash control, and b) with backlash control 3.2 Comparison of Transmission Error with and without maximum value, minimum value, and backlash are Backlash Control compared, as shown in Table 3. It is easy to know that the average, maximum, According to the TE curves measured by all the above and minimum values of the TE with and without the simulations and experiments, the average value, backlash control are relatively consistent, while the instantaneous backlash is significantly reduced after the backlash control. Therefore, it can be seen that the instantaneous backlash control can reduce the TE at a certain time after the backlash control, and the transmission process is more stable, which proves the effectiveness of the electronically controlled anti-backlash mechanism. Since this article aims to realiz e and verify the anti-backlash function of the electronically controlled anti-backlash mechanism, the machining and assembly accuracy of gears and other parts adopts ordinary accuracy (gear accuracy is 7 levels). Therefore, in the entire transmission process with and without the backlash control, the main factor that affects the large-period TE is the eccentricity error. Table 3. Comparison of TE with and without backlash control Rotating speed Backlash control Avg TE [°] Max Min Backlash [°] 10 rpm With With Experiment Simulation Experiment 0.0412 0.0427 0.0425 0.1117 0.1118 0.1183 0.016 -0.0135 -0.0136 --0.0342 out Simulation 0.0419 0.1121 -0.0135 0.0262 30 rpm With With Experiment Simulation Experiment 0.0437 0.0430 0.0408 0.1316 0.1253 0.1316 -0.0145 -0.0125 -0.0154 --0.0471 out Simulation 0.0407 0.1253 -0.0125 0.0443 60 rpm With With Experiment Simulation Experiment 0.0345 0.0385 0.0347 0.1228 0.1278 0.1302 -0.0175 -0.0120 -0.0194 --0.0412 out Simulation 0.0362 0.1279 -0.0126 0.0443 100 rpm With With Experiment Simulation Experiment 0.0368 0.0387 0.0331 0.1224 0.1275 0.1381 -0.018 -0.0125 -0.0213 --0.0437 out Simulation 0.0356 0.1276 -0.0126 0.0436 4 CONCLUSIONS (1) This paper has presented an electrical backlash control method based on the novel type of VTTIGP. According to the angular acceleration and angular displacement of the driving gear, backlash elimination can be realiz ed at the time of reversal by adjusting the relative axial displacement between the driving and driven gear. (2) The continuous theoretical backlash curve was attained through the eccentricity and its initial phase of the novel type of VTTIGP, which can be obtained by fitting the experimental TE curve through the developed no-load TE equation. Comparing the theoretical backlash to the experiment one, it can be seen that the backlash curve obtained by the theoretical equation agrees with the backlash value determined by MR BM. (3) An experimental backlash control investigation was also carried out, which showed that under different input speeds and at various reverse positions, the TE curves without backlash control both have a jump phenomenon at the time of reverse; however, after backlash control, the jump phenomenon in TE curves improves. This demonstrates that this backlash control method can eliminate backlash effectively at the time of reversal. Additionally, the TE curves of the simulation and experiment agreed well with each other, which indicates that the constructed co-simulation model of the novel type of VTTIGP is consistent with the experimental rig. The conclusion shows that under no-load conditions, the theory and experiment have achieved good agreement. However, the load has a certain effect on the backlash. Therefore, it is necessary to establish the backlash control model of the novel type of VTTIGP under load and carry out related experimental research. More related topics will be discussed in the follow-up. 5 ACKNOWLEDGEMENTS This work is supported by the National Natural Science Foundation of China (Grant No. 92048201 ) and the Fundamental R esearch Funds for the Central Universities (Project No.2021C DJJMR H-008) . The authors thank the reviewer for his/her valuable comments on the manuscript. 6 REFERENCES [1] Wang, G., Chen, L., Yu, L., Zou, S. (2017). Research on the dynamic transmission error of a spur gear pair with eccentricities by finite element method. Mechanism and Machine Theory, vol. 109, p. 1-13, DOI:10.1016/j. mechmachtheory.2016.11.006. [2] Wang, G., Su, L., Zou, S. (2020). Uneven load contact dynamic modelling and transmission error analysis of a 2K-V reducer with eccentricity excitation. Strojniški vestnik - Journal of Mechanical Engineering, vol. 66, no. 2, p. 91-104, DOI:10.5545/sv-jme.2019.6298. [3] Yang, Q., Liu, T., Wu, X., Deng, Y., Chen, Q. (2021). A planetary gear reducer backlash identification based on servo motor current signal and optimized fisher discriminant analysis. ISA Transaction, vol. 112, p. 350-362, DOI:10.1016/j. isatra.2020.12.016. [4] Liu, W., Zhao, H., Lin, T., Gao, B., Yang, Y. (2020). Vibration characteristic analysis of gearbox based on dynamic excitation with eccentricity error. Journal of Mechanical Science and Technology, vol. 34, p. 4545-4562, DOI:10.1007/s12206-020­1014-6. [5] Ren, X., Li, D., Sun, G., Zhao, W. (2016). Eso-based adaptive robust control of dual motor driving servo system. Asian Journal of Control, vol. 18, n. 6, p. 2358-2365, DOI:10.1002/ asjc.1325. [6] Dong, R., Tan, Y., Chen, H. (2009). Recursive identification for dynamic systems with backlash. Asian Journal of Control, vol. 12, no. 1, p. 26-38, DOI:10.1002/asjc.157. [7] Warnecke, M., Jouaneh, M. (2003). Backlash compensation in gear trains by means of open-loop modification of the input trajectory. Journal of Mechanical Design, vol. 125, no. 3, p. 620-624, DOI:10.1115/1.1596241. [8] Yang, Q., Liu, T., Wu, X., Deng, Y. (2020). Gear backlash detection and evaluation based on current characteristic extraction and selection. IEEE Access, vol. 8, p. 107161­107176, DOI:10.1109/access.2020.2999478. [9] Yu, C. (2016). Analysis and study on eliminating clearance structure of mechanical precision transmission. Mechanical Research & Application, vol. 29, no. 1, p. 8-11. (in Chinese) [10] Nordin, M., Gutman, P.-O. (2002). Controlling mechanical systems with backlash-a survey. Automatica, vol. 38, no. 10, p. 1633-1649, DOI:10.1016/s0005-1098(02)00047-x. [11] Michaelec, G.W. (1966). Precision Gearing-Theory and Practice. John Wiley & Sons, New York. [12] Beam, A.S. (1954). Beveloid gear. Machine Design, vol. 26, no. 12, p. 220-238. [13] Chen, Q., Song, C., Zhu, C., Du, X., Ni, G. (2017). Manufacturing and contact characteristics analysis of internal straight beveloid gear pair. Mechanism and Machine Theory, vol. 114, p. 60-73, DOI:10.1016/j.mechmachtheory.2017.04.002. [14] Sun, R., Song, C., Zhu, C., Yang, X., Li, X. (2021). Computational study of pitting defect influence on mesh stiffness for straight beveloid gear. Engineering Failure Analysis, vol. 119, p. 1-15, DOI:10.1016/j.engfailanal.2020.104971. [15] Do, T.P., Ziegler, P., Eberhard, P. (2015). Review on contact simulation of beveloid and cycloid gears and application of a modern approach to treat deformations. Mathematical and Computer Modelling of Dynamical Systems, vol. 21, no. 4, p. 359-388, DOI:10.1080/13873954.2015.1012838. [16] Yang, X., Song, C., Zhu, C., Liu, S. (2018). Tooth surface deviation and mesh analysis of beveloid gears with parallel axis considering machine tool adjustment errors. Journal of Advanced Mechanical Design, Systems, and Manufacturing, vol. 12, no. 4, p. 1-15, DOI:10.1299/jamdsm.2018 jamdsm0082. [17] Cao, B., Li, G. (2021). Computerized design of plunge shaving tool for beveloid gears and plunge shaving characteristic analysis. Mechanism and Machine Theory, vol. 161, p. 1-16, DOI:10.1016/j. mechmachtheory.2021.104325. [18] Brauer, J. (2005). Transmission error in anti-backlash conical involute gear transmissions: A global-local fe approach. Finite Elements in Analysis and Design, vol. 41, no. 5, p. 431-457, DOI:10.1016/j.finel.2004.04.007. [19] Song, C., Zhou, S., Zhu, C., Yang, X., Li, Z., Sun, R. (2018). Modeling and analysis of mesh stiffness for straight beveloid gear with parallel axes based on potential energy method. Journal of Advanced Mechanical Design, Systems, and Manufacturing, vol. 12, no. 7, p. 1-14, DOI:10.1299/ jamdsm.2018jamdsm0122. [20] Sun, R., Song, C., Zhu, C., Wang, Y., Yang, X. (2020). Computational studies on mesh stiffness of paralleled helical beveloid gear pair. International Journal of Precision Engineering and Manufacturing, vol. 22, no. 1, p. 123-137, DOI:10.1007/s12541-020-00452-3. [21] Li, G., Wu J., Li H. (2000). Design and calculation of parallel shaft internal meshing beveloid gear. Chinese Mechanical Engineering, vol. 8, no.11, p. 886-889. (in Chinese) [22] Yan, H., Huang, Q., Xing, J., Zhou, F. (2017). The analysis of the efficiency of 3K type conical involute planetary gears transmission reducer. Machinery Design & Manufacture, vol. 5, p. 158-161. (in Chinese) [23] Zheng, Z., Chen, Y., Chen, B., Du, X., Li, C. (2020). Meshing performance investigations on a novel point-contact hourglass worm drive with backlash-adjustable. Mechanism and Machine Theory, vol. 149, p. 1-17, DOI:10.1016/j. mechmachtheory.2020.103841. [24] Ye, X., Chen, Y., Lu, B., Luo, W., Chen, B. (2022). Study on a novel backlash-adjustable worm drive via the involute helical beveloid gear meshing with dual-lead involute cylindrical worm. Mechanism and Machine Theory, vol. 167, p. 1-21, DOI:10.1016/j.mechmachtheory.2021.104466. [25] Tsai, L. W., Chang, S.L. (1994). Backlash control via redundant drives: An experimental verification. Journal of Mechanical Design, vol. 116, no. 3, p. 963-967, DOI:10.1115/1.2919477. [26] dos SantosCoelho, L., Cunha, M.A.B. (2011). Adaptive cascade control of a hydraulic actuator with an adaptive dead-zone compensation and optimization based on evolutionary algorithms. Expert Systems with Applications, vol. 38, n. 10, p. 12262-12269, DOI:10.1016/j.eswa.2011.04.004. [27] Rostalski, P., Besselmann, T., Baric, M., Belzen, F. V., Morari, M. (2007). A hybrid approach to modelling, control and state estimation of mechanical systems with backlash. International Journal of Control, vol. 80, n. 11, p. 1729-1740, DOI:10.1080/00207170701493985. [28] Na, J., Ren, X., Herrmann, G., Qiao, Z. (2011). Adaptive neural dynamic surface control for servo systems with unknown dead-zone. Control Engineering Practice, vol. 19, no. 11, p. 1328-1343, DOI:10.1016/j.conengprac.2011.07.005. [29] Shi, Z., Zuo, Z. (2015). Backstepping control for gear transmission servo systems with backlash non-linearity. IEEE Transactions on Automation Science and Engineering, vol. 12, no. 2, p. 752-757, DOI:10.1109/tase.2014.2369430. [30] Sato, K., Murayama, Y., Imada, S., Shimokohbe, A. (1995). Control and elimination of lead screw backlash for ultra-precision positioning. JSME International Journal. Ser. C, Dynamics, Control, Robotics, Design and Manufacturing, vol. 38, no. 1, p. 36-41, DOI:10.1299/jsmec1993.38.36. [31] Chen, W., Wu, Y., Du, R., Wu, X. (2014). Fault diagnosis and fault tolerant control for the servo system driven by two motors synchronously. Control Theory & Application, vol. 31, no. 1, p. 27-34. (in Chinese) [32] Zhang, Z., Peng, C., Shao, S. (2017). A Rotation speed offset regulation and model predictive control based anti-backlash control approach for dual motor servo system, MICROMOTORS, vol. 50, no. 1, p. 54-58. (in Chinese) [33] Jiang, Y., Zhang, W., Liu, X., Jin, B. (2020). simulation and application of anti-backlash technique with dual motor in tandem manipulators. China Mechanical Engineering, vol. 31, no. 16, p. 1991-1997. (in Chinese) [34] Litvin, F.L. (2008). Gear Geometry and Applied Theory. Shanghai Science and Technology Press, Shanghai. [35] Wu, X.T. (1982). Gear Meshing Principle. Machinery Industry Press, China. [36] Tao, L. (2005). Performance comparison analysis of stepped double involute gear and ordinary involute gear. Chongqing University, Chongqing . [37] Yu, L., Wang, G., Zou, S. (2017). The calculation of meshing efficiency of a new type of conical involute gear. Strojniški vestnik - Journal of Mechanical Engineering, vol. 63, no. 5, p. 320-330, DOI:10.5545/sv-jme.2016.3843. [38] Yu, L., Wang, G., Zou, S. (2018). The experimental research on gear eccentricity error of backlash-compensation gear device based on transmission error. International Journal of Precision Engineering and Manufacturing, vol. 19, p. 5-12, DOI:10.1007/s12541-018-0001-7. [39] Michalec, G.W. (1973). Correlation of probabilistic backlash with measurements. Mechanism and Machine Theory, vol. 8, no. 2, p. 161-173, DOI:10.1016/0094-114x(73)90050-5. [40] Wu, C. (1982). Research on transmission error by eccentric gear. Journal of Nanjing Institute of Technology, vol. 4, p. 133­ 145. (in Chinese) [41] Rocca, E., Russo, R. (2011). Theoretical and experimental investigation into the influence of the periodic backlash fluctuations on the gear rattle. Journal of Sound and Vibration, vol. 330, no. 20, p. 4738-4752, DOI:10.1016/j. jsv.2011.04.008. Strojniški vestnik - Journal of Mechanical Engineering 68(2022)2 Vsebina Vsebina Strojniški vestnik - Journal of Mechanical Engineering letnik 68, (2022), številka 2 L jubljana, februar 2022 ISSN 0039-2480 aa eeo R az širjeni povz etki (extended abstracts) Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš: Modeliranje kavitacijske eroz ije na radialno divergentnem testnem odseku z uporabo R ANS SI 1 1 Anz e Sitar: Nova nestacionarna metoda merjenja prepustnosti preverjena na tkaninah SI 12 Leilei hao, ueei u, Jianhu Cao, eiei hou: elinearen sklopljeni dinamicni model sistema voznik-sedež-kabina in napoved biomehanskega vedenja SI 13 Marek ozniak, Adam Rylski, Krzysztof Siczek: Meritve obrabe komponent koncnikov SI 14 Guangjian Wang, Dongdong Z hu, Shuaidong Z ou, Yujiang Jiang, X in Tian: Simulacija in eksperimentalna raziskava elektricne regulacije razstopa pri novi evolventni zobniški dvojici z variabilno debelino z ob SI 15 Strojniški vestnik - Journal of Mechanical Engineering 68(2022)2, SI 11 Prejeto v recenzijo: 2021-08-10 © 2022 Avtorji. CC BY 4.0 Int. Prejeto v recenzijo: 2021-12-17 Odobreno za objavo:2022-01-13 Modeliranje kavitacijske eroz ije na radialno divergentnem testnem odseku z uporabo R ANS * Luka Kevorkijan – L uka Lešnik – I gnacijo Biluš Univerz a v Mariboru, Fakulteta z a strojništvo, Slovenija Kavitacija je pojav uparjanja kapljevine v hidravlicnih sistemih zaradi padca tlaka pod tlak uparjanja. Parne strukture, ki so izpostavljene višjemu tlaku kapljevine v okolici, so podvržene kondenzaciji in posledicno kolapsirajo. Pri tem se iz inženirskega vidika pojavi vec neželjenih ucinkov, med drugim kavitacijska erozija. a izboljšave obratovanja hidravlicnih sistemov je želja boljše poznavanje in napovedovanje pojava kavitacijske eroz ije. a numericno modeliranje tega v splošnem nezaželjenega pojava, povezanega s poškodbami materiala, so na voljo razlicni validirani matematicni modeli prenosa parne faze. V preteklosti so se za napoved kavitacijske erozije uporabljali razlicni, predvsem eksperimentalni pristopi. razvojem numericnih simulacij so se v zadnjem casu razvili razlicni numericni modeli za napoved kavitacijske erozije, ki temeljijo na razlicnih hipotezah, kot sta hipoteza tlacnega vala in hipoteza mikrocurka ali kombinaciji obeh. V pricujocem clanku je bila izvedena numericna simulacija turbulentnega toka s pojavom kavitacije z namenom napovedi kavitacijske eroz ije. Z a dvofaz ni tok je uporabljen model homogene z mesi, turbulenca je modelirana po pristopu R ANS. Z a modeliranje kavitacije v toku je uporabljen Schnerr-Sauer prenosni kavitacijski model. umericni model je nadgrajen z modelom erozijske potencialne energije 1 do 3, ki temelji na hipotezi tlacnega vala. V erozijskem modelu je izvedena sprostitev tlacnega vala na dva nacina, s sprotnim sprošcanjem kavitacijske potencialne energije ob poteku kolapsa parne faze in s fokusiranim sprošcanjem kavitacijske potencialne energije na koncu kolapsa. adgrajen numericni model, izražen s tremi razlicnimi matematicnimi forumalcijami eroz ijske potencialne energije, je uporabljen na primeru radialno divergentnega testnega odseka [ 4] . R ez ultati simulacije so validirani s primerjavo z eksperimentalnimi rez ultati drugih avtorjev [ 4] . Vsklopu študijeje ugotovljeno, da se porazdelitev nakopicene energije na površini sklada z eksperimentalnimi rezultati, vendar pa obstajajo razlike med razlicnimi matematicnimi formulacijami erozijskega modela. Primerjava napovedanega položaja poškodba ni pokazala razlik med fokusiranim in sprotnim sprošcanjem erozijske potencialne energije, vendar pa je v primeru fokusiranja napovedana eroz ija bolj koncentrirana, kar predstavlja korak k boljši napovedi kavitacijske erozije. a natancno napoved kavitacijske erozije je kljucen natancen opis kavitacijske dinamike. Schnerr-Sauer kavitacijski model pravilno napove razsežnost in dinamiko kavitacije. Predstavljena metodologija za napoved kavitacijske erozije izkazuje primernost za nadaljnjo uporabo in možnost nadgradnje v smeri modeliranja odz iva materiala površine iz postavljene kavitacijski eroz iji. le eede aitacia eroia eroia oteciala eeria eria ilacia Strojniški vestnik - Journal of Mechanical Engineering 68(2022)2, SI 12 Prejeto v recenzijo: 2021-10-26 © 2022 Avtorji. CC BY 4.0 Int. Prejeto popravljeno: 2021-12-29 Odobreno za objavo: 2022-01-11 Nova nestacionarna metoda merjenja prepustnosti preverjena na tkaninah * Anz e Sitar Fundacija Bruno Kessler, Italija Merjenje prepustnosti materialov je pomembno na številnih podrocjih in nekatera od njih že imajo na voljo ustrezne merilne tehnike. a podrocju liofiliziranih produktov, poroznih umetnih kosti, 3D tiskanih materialov itd. imajo proizvodi lahko povsem edinstveno obliko, kar onemogoca uporabo standardizirane opreme. Iz tega razloga je bila razvita nova metoda za dolocanje prepustnosti materiala z analizo visokofrekvencnega merjenja nestacionarnega tlaka med permeacijo delovne tekocine. Metoda je zaradi hitrosti zajema in analize tlacnega signala primerna tudi za meritve v realnem casu. a fizikalno izhodišce smo na podlagi režima toka delovne tekocine uporabili Darcy-jev zakon, ki smo ga lahko poenostavili z aradi uporabe primerjalne analiz e iz merkov, saj se viskoz nost z raka, površina tkanine in prostornina toka zraka ni spreminjala med razlicnimi meritvami. Primerjalna analiza merilnih rezultatov je omogocila enostavnejši preracun prepustnosti, ki ne zahteva poznavanja niti dodatnih termodinamicnih lastnosti delovne tekocine ali merjenega materiala niti podrobnega poznavanju fizikalnega ozadja toka tekocin. V primeru obravnavanih tkanin smo za primerjalni parameter dolocili integralni parameter dx /PT , ki vkljucuje debelino tkanine dx in normalizirano vrednost produkta tlaka in casa PT . Primerjalni parameter iz haja iz iz merjenega nor m nor m tlacnega signala in bi lahko bil sestavljen iz naklona, najvišje vrednosti, odvodov, integrala ali kombinacije naštetega, ce bi se drug primerjalni parameter v konkretnem primeru izkazal z boljšim determinacijskim koeficientom. R az vita metoda je nova, z ato so bili pridobljeni merilni rez ultati primerjani z iz merki standardne naprave za merjenje zracne prepustnosti tkanin. Tako pri novi nestacionarni metodi merjenja prepustnosti kot tudi pri standardizirani stacionarni metodi smo uporabljali zrak kot delovno tekocino, ki prehaja skozi tkanine. Izmerjene prepustnosti so z našale med 8 Da in 5 0 Da z a vseh pet obravnavanih tkanin. Primerjalna analiz a linearnega regresijskega modela je podala z elo spodbudne rez ultate, saj je imela determinacijski koeficient R 2 = 0.98 z uporabo le enega iz merka vsake iz med petih tkanin, ki je bil pridobljen z novo nestacionarno metodo. Poleg merjenja prepustnosti tkanin je bila nova metoda preiz kušena tudi na liofiliz iranih produktih v steklenih vialah, kar je najpogostejši nacin uporabe procesa zmrzovalnega sušenja farmacevtskih proizvodov. Potrdili smo možnost uporabe razvite nestacionarne metode za merjenje prepustnosti liofiliziranih farmacevtskih produktov, saj smo lahko z az navali raz like v prepustnosti liofiliz irane 5 wt% in 12 wt% vodne raz topine manitola, ki imata razlicne prepustnosti zaradi razlicnih poroznosti, le-ta pa je odvisna od zacetne kolicine manitola v raztopini. Dolocitev prehodnosti liofiliziranih farmacevtskih produktov je kljucnega pomen za optimizacijo procesa z mrz ovalnega sušenja. R az vita metoda je po mojem vedenju prva, ki lahko prepustnost tovrstnih farmacevtskih produktov iz meri neposredno. le eede oa erila etoda ereailot retot etacioara etoda erea tkanine, liofiliz irani proiz vodi Strojniški vestnik - Journal of Mechanical Engineering 68(2022)2, SI 13 Prejeto v recenzijo: 2021-10-07 © 2022 Avtorji. CC BY 4.0 Int. Prejeto popravljeno: 2021-11-22 Prejeto popravljeno: 2021-12-16 elieare lolei diaii odel itea oiede kabina in napoved biomehanskega vedenja Leilei Z hao1 – Yuewei Yu1,* – J ianhu Cao1 – Yuewei Yu2 1 Tehniška univerz a v Shandongu, Šola z a transport in avtomobilsko tehniko, Kitajska 2Poklicni inštitut ibo, ola za avtomobilsko inženirstvo, Kitajska Biomehanski odziv sistema voznik-tovornjak v dinamicnem okolju je eden od glavnih dejavnikov pri projektiranju in upravljanju tovornjakov. a ovrednotenje voznega udobja sistema sedež-kabina je treba napovedati biomehanski odziv razlicnih delov voznikovega telesa v razlicnih smereh. Trenutno ni zanesljivih modelov in metod za ucinkovito napovedovanje karakteristicnega odziva sistema voznik-sedež-kabina. ato je bil na osnovi obstojecegabiodinamicnega modela sedecega cloveka s 7 prostostnimi stopnjami najprej ustvarjen nelinearen dinamicni model sistema voznik-sedež-kabina z 10 prostostnimi stopnjami in izpeljane so bile diferencialne enacbe za popis vibracij. V naslednjem koraku so bili za simulacijo in verifikacijo eksperimentalno pridobljeni signali vibracij tovornjaka med vožnjo po cesti. V tretjem koraku je predstavljen specificen postopek reševanja modela na osnovi integracijske metode emark-. elinearni koeficienti dušenja prednjih in zadnjih blažilnikov kabine so bili izmerjeni na preizkuševališcu. pravljene so bile simulacije na podlagi izmerjenih parametrov modela z z branimi signali vibracij šasije kot vhodnimi podatki. R ez ultati simulacije se ujemajo z rezultati preskusov, kar dokazuje, da lahko dinamicni model ucinkovito napove voznikove biomehanske odzive. Analiza s simulacijami je dala vec uporabnih zakljuckov. (1) Primerjava je pokaz ala, da se rez ultati simulacije ujemajo z rez ultati preskusov. Model je uporaben in njegova tocnost je sprejemljiva. Predlagani model lahko uspešno napove voznikov biomehanski odziv in osnovne dinamicne lastnosti tovornjakov, ki so med vožnjo po cesti podvrženi vplivu nakljucnih vzbujanj. (2) asovne vrste vertikalnih sil in pospeškov spodnjega dela trupa, drobovja, rok, zgornjega dela trupa in glave pod vplivom nakljucnih vzbujanj med vožnjo s tovornjakom so si podobne. (3) Energija vertikalnih vibracij voznikovega telesa med vožnjo po asfaltirani cesti je koncentrirana v nizkofrekvencnem obmocju (0 Hz do 10 Hz). Energija vibracij voznikovega telesa v visokofrekvencnem obmocju (10 Hz do 100 Hz) je prakticno enaka nic. Predlagani model odpravlja pomanjkljivosti obstojecega modela in lahko bolje simulira pospeške voznikovega hrbtnega naslona in kolkov z aradi vibracij. Novi model je z ato uporabnejši z a projektiranje sistemov vz metenja kabin. Primeren je tudi za analizo vibracij voznikovega celotnega telesa med vožnjo s tovornjakom kot teoreticno osnovo za dolocanje ucinkovitih ukrepov za preprecevanje poškodb voznikovega drobovja. Model ni primeren za analizo utrujanja cloveškega telesa zaradi vibracij. jegova pomanjkljivost je, da ne more tocno reproducirati realnih vzdolžnih pospeškov. Vnadaljnjih študijah bo zato treba upoštevati še vzdolžne prostostne stopnje kabine ter vdelati vodilni mehanizem sistema vzmetenja kabine v dinamicni model. le eede atooila teia aoot oe ite oiedeaia oioo drae diaio odelirae ioeai odiY Strojniški vestnik - Journal of Mechanical Engineering 68(2022)2, SI 14 Prejeto v recenzijo: 2021-09-06 © 2022 Avtorji. CC BY 4.0 Int. Prejeto popravljeno: 2022-01-10 Odobreno za objavo: 2022-02-02 erite orae ooet oioY Marek Woz niak1 – Adam R ylski2 – K rz ysz tof Sicz ek1,* 1 Tehniška univerza v Lodžu, ddelek za vozila in osnove konstruiranja strojev, Poljska 2Tehniška univerza v Lodžu, Inštitut za materiale in inženiring, Poljska amen raziskave je bila dolocitev intenzivnosti obrabe ležišc koncnikov iz materiala polioksimetilen (PM) v z gibu z jekleno kroglo. Intenzivnost obrabe je bila dolocena z naslednjimi predpostavkami: režim vožnje med popravili krmilnega sestava, realen model vožnje v ovinek. Opredeljeni so bili: model krmilnega mehanizma, ki omogoca dolocitev obremenitve cepa koncnika v odvisnosti od krmilnega kota, model obrabe kroglicnega cepa v odvisnosti od kontaktnega tlaka v analizirani coni zgiba, model za dolocitev deleža obremenitev, ki se prenašajo prek deformiranih grobih neravnin in s hidrodinamicnim maz anjem v pogojih mešanega trenja v analiz irani kontaktni coni, model za dolocitev tlaka v stisnjeni plasti masti v analizirani kontaktni coni. braba krogle je bila dolocena kot razlika v prostorninah nove krogle in obrabljene krogle v mejah maksimalne tangencialne zunanje površine. braba ležišca je bila dolocena kot razlika v prostorninah obrabljenega in novega ležišca v mejah minimalne tangencialne notranje površine. Meritve so bile opravljene z mikrometrom. a podlagi razvitih modelov in izmerjenih vrednosti volumetricne obrabe je bila numericno dolocena intenzivnost obrabe ležišca iz PM, temu pa je sledila primerjava z vrednostmi za izbrane kroglicne zgibe. Kontaktni tlak v analizirani kontaktni coni je bil dolocen po metodi koncnih elementov. Eksperimenti so pokazali, da je obraba krogle zgiba bistveno manjša od obrabe ležišca iz PM.a ležišce iz PM je bila ocenjena vrednost faktorja intenzivnosti obrabe k v višini 1,3810-19 mN-2. Ta 4 ustreza intenzivnosti linearne obrabe Ih v višini 5,10-9 , ki lahko mocno variira v odvisnosti od kontaktnega tlaka, drsne hitrosti, temperature, vlažnosti, delovnih pogojev in mazanja sestava koncnika. Ugotovljene vrednosti faktorja intenz ivnosti obrabe k in intenz ivnosti linearne obrabe Ih lahko odstopajo od dejanskih vrednosti, saj so mocno odvisne od metode dolocanja povprecnega kontaktnega tlaka v coni med kroglico in ležišcem. Povprecne vrednosti kontaktnega tlaka v analizirani kontaktni coni narašcajo prakticno linearno z aksialno obremenitvijo koncnika F . Povprecna vrednost kontaktnega tlaka p pri vrednosti sile 300 je znašala 0,21 MPa. Vrednosti kontaktnega tlaka pri obremenitvi kroglicnega zgiba z enako silo F so podobne tistim pri stiskanju plasti masti v reži med kroglo in ležišcem. Podobne so tudi vrednostim tlaka v obmocju hidrodinamicnega mazanja z mikroklini masti, ko se je podobno obremenjen kroglicni cep vrtel v ležišcu z najvecjo drsno hitrostjo 0,03 ms. Kontaktni tlak v kontaktnih conah med plasticno deformiranimi neravninami je dosegel vrednost 150 MPa. Pri izvedbi študije je bila uporabljena predpostavka, da so vrednosti vsebnosti vlage, tujkov in onesnaževal v masti ter temperature v analiz irani kontaktni coni konstantne ter da ne vplivajo na reološke lastnosti masti oz . na obrabo sticnih površin in same masti. Vpliv sprememb temperature, vrste masti in omenjenih parametrov na obrabo v analiz irani kontaktni coni bo predmet prihodnjih raz iskav. Avtorji so eksperimentalno ocenili vrednost faktorja intenzivnosti obrabe k za ležišce iz PM, ki je del kroglicnega zgiba. Teh podatkov ni bilo mogoce najti v obstojeci literaturi. Faktor ustreza intenzivnosti linearne obrabe Ih, ki je v literaturi le redko uporabljena z a kvantifikacijo procesa obrabljanja z gibov. Ocenjen je bil vpliv maz anja v kontaktni coni na poraz delitev obremenitev med neposrednimi stiki grobih neravnin in hidrodinamicnim mazanjem, ki ga zagotavljajo mikroklini masti v prostoru med neravninami. Analiziran je bil tudi vpliv stiskanja mazalnega filma na tlak v sloju masti med kroglo in ležišcem. le eede oi toa orae etoda oi eleeto ite rilea eaia kontaktni tlak, mast, stisnjen film Stefanowskiego Str. 1/5, 90-537 Lodz, Poljska, ks670907@p.lodz.pl Strojniški vestnik - Journal of Mechanical Engineering 68(2022)2, SI 15 Prejeto v recenzijo: 2021-10-25 © 2022 Avtorji. CC BY-NC 4.0 Int. Prejeto popravljeno: 2021-12-30 Odobreno za objavo: 2022-01-04 Simulacija in eksperimentalna raz iskava eletrie relacie ratoa ri oi eoleti z obniški dvojici z variabilno debelino z ob Guangjian Wang1,2 – D ongdong Z hu1,2 – S huaidong Z ou1,2,* – Yujiang Jiang1,2 – X in Tian1,2 1 Univerza v ongcingu, Državni laboratorij za mehanske prenose, Kitajska 2 Univerza v ongcingu, Kolidž za strojništvo in avtomobilsko tehniko, Kitajska Razstop pri preciznih prenosnikih neposredno vpliva na stabilnost, tocnost in obcutljivost servokrmiljenja, saj povzroca casovni zamik in nezveznosti v izhodnem gibanju. V clanku je predstavljena metoda za elektricno regulacijo razstopa v realnem casu pri novi evolventni zobniški dvojici z variabilno debelino zob (VTTIGP). Glavni vir nizkofrekvencne komponente napak prenosa in variabilnega razstopa pri preciznih prenosnikih je ekscentricnost zobnikov. Po dolocitvi krivulje skupnega razstopa za nove VTTIGP, ki jo sestavljata konstantna in periodicna casovno variabilna komponenta, je razstop mogoce regulirati in ga prilagajati z relativnimi aksialnimi premiki zobnikov. Regulacija poteka na osnovi zvezne krivulje razstopa in sprememb pospeška z realnocasovnim uravnavanjem medsebojnega položaja pogonskega in gnanega zobnika. lanek predstavlja princip delovanja, lastnosti in teoreticno enacbo površine zob z variabilno debelino (VTTIG). Vprimerjavi s tradicionalnimi zobniki z variabilno debelino zob ima novi tip VTTIG v vsakem precnem prerezu enak vznožni in temenski krog. To pomeni, da pri odpravljanju razstopa ne more priti do prekrivanja zaradi prevelike relativne dolžine premika med zobnikoma. Podana je teoreticna izpeljava enacbe napake prenosa v neobremenjenem stanju in casovno spremenljivega razstopa za VTTIGP. vezna krivulja razstopa zobniške dvojice je bila dolocena s pomocjo teoreticne enacbe za razstop, ekscentricnosti ter zacetnih faz, pridobljenih iz krivulje napake razstopa, ki je bila dolocena eksperimentalno s senzorji kotnega pomika in enacbo za napako prenosa. S programsko opremo RecurDyn in pripadajocim modulom za regulacijo Colink je bil ustvarjen kosimulacijski model VTTIGP in sistema za regulacijo razstopa. Preucena je bila metoda elektricne regulacije razstopa z oceno pospeškov za primer casovno spremenljivega razstopa. Koncno je bila opravljena še eksperimentalna študija sistema za regulacijo razstopa z razlicnimi položaji spremembe smeri prenosa in vhodnimi hitrostmi. V ta namen so bili uporabljeni senz orji kotnega pomika in linearni aktuatorji. Rezultati za razlicne vhodne hitrosti in položaje spremembe smeri prenosa za primer brez regulacije razstopa so pokaz ali skok na krivuljah napake prenosa. Ta skok na krivulji se z manjša z uvedbo regulacije raz stopa in s tem je bilo dokazano, da regulacija ucinkovito odpravlja razstop v trenutku spremembe smeri vrtenja. Primerjava krivulj razstopa, ki so bile dolocene teoreticno in eksperimentalno po metodi mehanskega povratka, je pokazala dobro ujemanje. Dobro se ujemajo tudi krivulje napake prenosa, dolocene s simulacijo in eksperimentalno. Iz tega sledi sklep, da je ustvarjeni kosimulacijski model VTTIGP skladen s sistemom na preizkuševališcu. lanek se dotika nekaterih vidikov, ki še niso bili obravnavani v strokovni literaturi, kot so princip delovanja, lastnosti in teoreticna enacba površine zob novega VTTIG, ter metoda elektricne regulacije razstopa na osnovi ocene sprememb pospeška. R az iskava je bila opravljena v pogojih brez obremenitve, teorija pa se dobro ujema z rezultati eksperimentov. bremenitve vplivajo na razstop in zato bo treba v prihodnjih raziskavah dolociti model z a regulacijo raz stopa pri VTTIGP pod obremenitvijo in opraviti ustrez ne eksperimentalne raz iskave. Predstavljena metoda za elektricno regulacijo razstopa je primerna za nadaljnje raziskave in razvoj z namenom izboljšanja tocnosti prenosa pri preciznih servopogonskih sistemih. le eede odraa ratoa eletria relacia aoo ariaile rato aaa reoa z obniška dvojica z variabilno debelino z ob Guide for Authors All manuscripts must be in English. Pages should be numbered sequentially. The manuscript should be composed in accordance with the Article Template given above. The suggested length of contributions is 10 to 20 pages. Longer contributions will only be accepted if authors provide justification in a cover letter. For full instructions see the Information for Authors section on the journal’s website: http://en.sv-jme.eu . SUBMISSION: Submission to SV-JME is made with the implicit understanding that neither the manuscript nor the essence of its content has been published previously either in whole or in part and that it is not being considered for publication elsewhere. All the listed authors should have agreed on the content and the corresponding (submitting) author is responsible for having ensured that this agreement has been reached. The acceptance of an article is based entirely on its scientific merit, as judged by peer review. Scientific articles comprising simulations only will not be accepted for publication; simulations must be accompanied by experimental results carried out to confirm or deny the accuracy of the simulation. Every manuscript submitted to the SV-JME undergoes a peer-review process. The authors are kindly invited to submit the paper through our web site: http://ojs.sv­jme.eu. The Author is able to track the submission through the editorial process - as well as participate in the copyediting and proofreading of submissions accepted for publication - by logging in, and using the username and password provided. SUBMISSION CONTENT: The typical submission material consists of: - A manuscript (A PDF file, with title, all authors with affiliations, abstract, keywords, highlights, inserted figures and tables and references), - Supplementary files: • a manuscript in a WORD file format • a cover letter (please see instructions for composing the cover letter) • a ZIP file containing figures in high resolution in one of the graphical formats (please see instructions for preparing the figure files) • possible appendicies (optional), cover materials, video materials, etc. Incomplete or improperly prepared submissions will be rejected with explanatory comments provided. In this case we will kindly ask the authors to carefully read the Information for Authors and to resubmit their manuscripts taking into consideration our comments. COVER LETTER INSTRUCTIONS: Please add a cover letter stating the following information about the submitted paper: 1. Paper title, list of authors and their affiliations. One corresponding author should be provided. 2. Type of paper: original scientific paper (1.01), review scientific paper (1.02) or short scientific paper (1.03). 3. A declaration that neither the manuscript nor the essence of its content has been pub­lished in whole or in part previously and that it is not being considered for publication elsewhere. 4. State the value of the paper or its practical, theoretical and scientific implications. What is new in the paper with respect to the state-of-the-art in the published papers? Do not repeat the content of your abstract for this purpose. 5. We kindly ask you to suggest at least two reviewers for your paper and give us their names, their full affiliation and contact information, and their scientific research inter­est. The suggested reviewers should have at least two relevant references (with an impact factor) to the scientific field concerned; they should not be from the same country as the authors and should have no close connection with the authors. FORMAT OF THE MANUSCRIPT: The manuscript should be composed in accordance with the Article Template. The manuscript should be written in the following format: - A Title that adequately describes the content of the manuscript. - A list of Authors and their affiliations. - An Abstract that should not exceed 250 words. The Abstract should state the principal objectives and the scope of the investigation, as well as the methodology employed. It should summarize the results and state the principal conclusions. - 4 to 6 significant key words should follow the abstract to aid indexing. - 4 to 6 highlights; a short collection of bullet points that convey the core findings and provide readers with a quick textual overview of the article. These four to six bullet points should describe the essence of the research (e.g. results or conclusions) and high­light what is distinctive about it. - An Introduction that should provide a review of recent literature and sufficient back­ground information to allow the results of the article to be understood and evaluated. - A Methods section detailing the theoretical or experimental methods used. - An Experimental section that should provide details of the experimental set-up and the methods used to obtain the results. - A Results section that should clearly and concisely present the data, using figures and tables where appropriate. - A Discussion section that should describe the relationships and generalizations shown by the results and discuss the significance of the results, making comparisons with pre­viously published work. (It may be appropriate to combine the Results and Discussion sections into a single section to improve clarity.) - A Conclusions section that should present one or more conclusions drawn from the results and subsequent discussion and should not duplicate the Abstract. - Acknowledgement (optional) of collaboration or preparation assistance may be includ­ed. Please note the source of funding for the research. - Nomenclature (optional). Papers with many symbols should have a nomenclature that defines all symbols with units, inserted above the references. If one is used, it must con­tain all the symbols used in the manuscript and the definitions should not be repeated in the text. In all cases, identify the symbols used if they are not widely recognized in the profession. Define acronyms in the text, not in the nomenclature. - References must be cited consecutively in the text using square brackets [1] and col­lected together in a reference list at the end of the manuscript. - Appendix(-icies) if any. v, T, n, etc.). Symbols for units that consist of letters should be in plain text (e.g. ms-1, K, min, mm, etc.). Please also see: http://physics.nist.gov/cuu/pdf/sp811.pdf . Abbreviations should be spelt out in full on first appearance followed by the abbreviation in parentheses, e.g. variable time geometry (VTG). The meaning of symbols and units belonging to symbols should be explained in each case or cited in a nomenclature section at the end of the manuscript before the References. Figures (figures, graphs, illustrations digital images, photographs) must be cited in consecutive numerical order in the text and referred to in both the text and the captions as Fig. 1, Fig. 2, etc. Figures should be prepared without borders and on white grounding and should be sent separately in their original formats. If a figure is composed of several parts, please mark each part with a), b), c), etc. and provide an explanation for each part in Figure caption. The caption should be self-explanatory. Letters and numbers should be readable (Arial or Times New Roman, min 6 pt with equal sizes and fonts in all figures). Graphics (submitted as supplementary files) may be exported in resolution good enough for printing (min. 300 dpi) in any common format, e.g. TIFF, BMP or JPG, PDF and should be named Fig1.jpg, Fig2.tif, etc. However, graphs and line drawings should be prepared as vector images, e.g. CDR, AI. Multi-curve graphs should have individual curves marked with a symbol or otherwise provide distinguishing differences using, for example, different thicknesses or dashing. Tables should carry separate titles and must be numbered in consecutive numerical order in the text and referred to in both the text and the captions as Table 1, Table 2, etc. In addition to the physical quantities, such as t (in italics), the units [s] (normal text) should be added in square brackets. Tables should not duplicate data found elsewhere in the manuscript. Tables should be prepared using a table editor and not inserted as a graphic. REFERENCES: A reference list must be included using the following information as a guide. Only cited text references are to be included. Each reference is to be referred to in the text by a number enclosed in a square bracket (i.e. [3] or [2] to [4] for more references; do not combine more than 3 references, explain each). No reference to the author is necessary. References must be numbered and ordered according to where they are first mentioned in the paper, not alphabetically. All references must be complete and accurate. Please add DOI code when available. Examples follow. Journal Papers: Surname 1, Initials, Surname 2, Initials (year). Title. Journal, volume, number, pages, DOI code. [1] Hackenschmidt, R., Alber-Laukant, B., Rieg, F. (2010). Simulating nonlinear materials under centrifugal forces by using intelligent cross-linked simulations. Strojniški vest-nik - Journal of Mechanical Engineering, vol. 57, no. 7-8, p. 531-538, DOI:10.5545/sv­jme.2011.013. Journal titles should not be abbreviated. Note that journal title is set in italics. Books: Surname 1, Initials, Surname 2, Initials (year). Title. Publisher, place of publication. [2] Groover, M.P. (2007). Fundamentals of Modern Manufacturing. John Wiley & Sons, Hoboken. Note that the title of the book is italicized. Chapters in Books: Surname 1, Initials, Surname 2, Initials (year). Chapter title. Editor(s) of book, book title. Publisher, place of publication, pages. [3] Carbone, G., Ceccarelli, M. (2005). Legged robotic systems. Kordi˜, V., Lazinica, A., Merdan, M. (Eds.), Cutting Edge Robotics. Pro literatur Verlag, Mammendorf, p. 553­576. Proceedings Papers: Surname 1, Initials, Surname 2, Initials (year). Paper title. Proceedings title, pages. [4] Štefani˜, N., Martin°evi˜-Miki˜, S., Tošanovi˜, N. (2009). Applied lean system in process industry. MOTSP Conference Proceedings, p. 422-427. Standards: Standard-Code (year). Title. Organisation. Place. [5] ISO/DIS 16000-6.2:2002. Indoor Air – Part 6: Determination of Volatile Organic Com­pounds in Indoor and Chamber Air by Active Sampling on TENAX TA Sorbent, Thermal Desorption and Gas Chromatography using MSD/FID. International Organization for Standardization. Geneva. WWW pages: Surname, Initials or Company name. Title, from http://address, date of access. [6] Rockwell Automation. Arena, from http://www.arenasimulation.com, accessed on 2009­09-07. EXTENDED ABSTRACT: When the paper is accepted for publishing, the authors will be requested to send an extended abstract (approx. one A4 page or 3500 to 4000 characters or approx. 600 words). The instruction for composing the extended abstract are published on-line: http://www.sv-jme.eu/information-for-authors/ . COPYRIGHT: Authors submitting a manuscript do so on the understanding that the work has not been published before, is not being considered for publication elsewhere and has been read and approved by all authors. The submission of the manuscript by the authors means that the authors automatically agree to publish the paper uder CC-BY 4.0 Int. or CC-BY-NC 4.0 Int. when the manuscript is accepted for publication. All accepted manuscripts must be accompanied by a Copyright Agreement, which should be sent to the editor. The work should be original work by the authors and not be published elsewhere in any language without the written consent of the publisher. The proof will be sent to the author showing the final layout of the article. Proof correction must be minimal and executed quickly. Thus it is essential that manuscripts are accurate when submitted. Authors can track the status of their accepted articles on https://en.sv-jme.eu/. PUBLICATION FEE: Authors will be asked to pay a publication fee for each article prior to the article appearing in the journal. However, this fee only needs to be paid after the article has been accepted for publishing. The fee is 380 EUR (for articles with maximum of 6 pages), 470 EUR (for articles with maximum of 10 pages), plus 50 EUR for each additional page. The additional cost for a color page is 90.00 EUR (only for a journal hard copy; optional upon author’s request). These fees do not include tax. SPECIAL NOTES Units: The SI system of units for nomenclature, symbols and abbreviations should be Strojniški vestnik -Journal of Mechanical Engineering followed closely. Symbols for physical quantities in the text should be written in italics (e.g. Ašker°eva 6, 1000 Ljubljana, Slovenia, e-mail: info@sv-jme.eu http://www.sv-jme.eu Contents Papers 71 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš: Cavitation erosion modelling on a radial divergent test section using RANS 82 Anze Sitar: Novel Unsteady State Method of Measuring Permeability Tested on Fabric Materials 90 Leilei Zhao, Yuewei Yu, Jianhu Cao, Weiwei Zhou: Nonlinear Coupled Dynamic Modelling of Driver-seat-cab System and Biomechanical Behaviour Prediction 101 Marek Wozniak, Adam Rylski, Krzysztof Siczek: The Measurement of the Wear of Tie Rod End Components 126 Guangjian Wang, Dongdong Zhu, Shuaidong Zou, Yujiang Jiang, Xin Tian: Simulation and Experimental Research on Electrical Control Anti-backlash Based on a Novel Type of Variable Tooth Thickness Involute Gear Pair