Since 1955 „ > Strojniški vestnik Journal of Mechanical Engineering Strojniški vestnik - Journal of Mechanical Engineering (SV-JME) Aim and Scope The international journal publishes original and (mini)review articles covering the concepts of materials science, mechanics, kinematics, thermodynamics, energy and environment, mechatronics and robotics, fluid mechanics, tribology, cybernetics, industrial engineering and structural analysis. The journal follows new trends and progress proven practice in the mechanical engineering and also in the closely related sciences as are electrical, civil and process engineering, medicine, microbiology, ecology, agriculture, transport systems, aviation, and others, thus creating a unique forum for interdisciplinary or multidisciplinary dialogue. The international conferences selected papers are welcome for publishing as a special issue of SV-JME with invited co-editor(s). Editor in Chief Vincenc Butala University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Technical Editor Pika Škraba University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Founding Editor Bojan Kraut University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Editorial Office University of Ljubljana, Faculty of Mechanical Engineering SV-JME, Aškerčeva 6, SI-1000 Ljubljana, Slovenia Phone: 386 (0)1 4771 137 Fax: 386 (0)1 2518 567 info@sv-jme.eu, http://www. sv-jme.eu Print: Koštomaj printing office, printed in 275 copies Founders and Publishers University of Ljubljana, Faculty of Mechanical Engineering, Slovenia University of Maribor, Faculty of Mechanical Engineering, Slovenia Association of Mechanical Engineers of Slovenia Chamber of Commerce and Industry of Slovenia, Metal Processing Industry Association President of Publishing Council Mitjan Kalin University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Vice-President of Publishing Council Bojan Dolšak University of Maribor, Faculty of Mechanical Engineering, Slovenia Strojniški vestnik Journal of Mechanical ''■s^^'' Engineering Cover: A sequence of images, taken in IR and visible spectrum, showing sludge layer during drying. An IR image and the next visible image were taken at the same time. Time between each image pair was 20 minutes. Image courtesy: Faculty of Mechanical Engineering, Ljubljana, Slovenia ISSN 0039-2480, ISSN 2536-2948 (online) International Editorial Board Kamil Arslan, Karabuk University, Turkey Hafiz Muhammad Ali, King Fahd U. of Petroleum & Minerals, Saudi Arabia Josep M. Bergada, Politechnical University of Catalonia, Spain Anton Bergant, Litostroj Power, Slovenia Miha Boltežar, University of Ljubljana, Slovenia Filippo Cianetti, University of Perugia, Italy Janez Diaci, University of Ljubljana, Slovenia Anselmo Eduardo Diniz, State University of Campinas, Brazil Jožef Duhovnik, University of Ljubljana, Slovenia Igor Emri, University of Ljubljana, Slovenia Imre Felde, Obuda University, Faculty of Informatics, Hungary Janez Grum, University of Ljubljana, Slovenia Imre Horvath, Delft University of Technology, The Netherlands Aleš Hribernik, University of Maribor, Slovenia Soichi Ibaraki, Kyoto University, Department of Micro Eng., Japan Julius Kaplunov, Brunel University, West London, UK Iyas Khader, Fraunhofer Institute for Mechanics of Materials, Germany Jernej Klemenc, University of Ljubljana, Slovenia Milan Kljajin, J.J. Strossmayer University of Osijek, Croatia Peter Krajnik, Chalmers University of Technology, Sweden Janez Kušar, University of Ljubljana, Slovenia Gorazd Lojen, University of Maribor, Slovenia Darko Lovrec, University of Maribor, Slovenia Thomas Lubben, University of Bremen, Germany Jure Marn, University of Maribor, Slovenia George K. Nikas, KADMOS Engineering, UK Tomaž Pepelnjak, University of Ljubljana, Slovenia Vladimir Popovič, University of Belgrade, Serbia Franci Pušavec, University of Ljubljana, Slovenia Mohammad Reza Safaei, Florida International University, USA Marco Sortino, University of Udine, Italy Branko Vasič, University of Belgrade, Serbia Arkady Voloshin, Lehigh University, Bethlehem, USA General information Strojniški vestnik - Journal of Mechanical Engineering is published in 11 issues per year (July and August is a double issue). Institutional prices include print & online access: institutional subscription price and foreign subscription €100,00 (the price of a single issue is €10,00); general public subscription and student subscription €50,00 (the price of a single issue is €5,00). Prices are exclusive of tax. Delivery is included in the price. The recipient is responsible for paying any import duties or taxes. Legal title passes to the customer on dispatch by our distributor. Single issues from current and recent volumes are available at the current single-issue price. To order the journal, please complete the form on our website. For submissions, subscriptions and all other information please visit: http:// www.sv-jme.eu. You can advertise on the inner and outer side of the back cover of the journal. The authors of the published papers are invited to send photos or pictures with short explanation for cover content. We would like to thank the reviewers who have taken part in the peer-review process. © 2020 Strojniški vestnik - Journal of Mechanical Engineering. All rights reserved. SV-JME is indexed / abstracted in: SCI-Expanded, Compendex, Inspec, ProQuest-CSA, SCOPUS, TEMA. The list of the remaining bases, in which SV-JME is indexed, is available on the website. The journal is subsidized by Slovenian Research Agency. Strojniški vestnik - Journal of Mechanical Engineering is available on https://www.sv-jme.eu. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9 Contents Contents Strojniški vestnik - Journal of Mechanical Engineering volume 66, (2020), number 9 Ljubljana, September 2020 ISSN 0039-2480 Published monthly Papers Andraž Lipolt, Brane Širok, Marko Hočevar, Lovrenc Novak: Convective Drying of Sewage Sludge Layer in Through-flow 481 Anna Trelka, Wojciech Zórawski, Anna Góral: Micro structure and Property Modification of Cold Sprayed Coatings Using Different Grain Sizes of Cr3C2-25(Ni20Cr) Composite Powder 494 Frank Goldschmidtboeing, Uwe Pelz, Karen Claire-Zimmet, Michael Wolf, Ralf Goerlach, Peter Woias: Bristle Motion, Forces, and Related Vertical Translation for a Novel Electric Toothbrush Design 505 Xiangyang Xu, Xupeng Fan, Peitang Wei, Baojun Yang: Research on the Lubrication Characteristics of Harmonic Gear Transmission Meshing Areas 513 Andres Lopez-Lopez, Jose Billennan Robles-Ocampo, Perla Yazmin Sevilla-Camacho, Orlando Lastres-Danguillecourt, Jesús Muniz, Bianca Yadira Perez-Sariñana, Sergio de la Cruz: Dynamic Instability of a Wind Turbine Blade Due to Large Deflections: An Experimental Validation 523 Yan Li, Lianhui Li: Enhancing the Optimization of the Selection of a Product Service System Scheme: A Digital Twin-Driven Framework 534 Urban Močnik, Bogdan Blagojevič, Simon Muliič: Numerical Analysis with Experimental Validation of Single-Phase Fluid Flow in a Dimple Pattern Heat Exchanger Channel 544 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 481-493 © 2020 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2020.6717 Original Scientific Paper Received for review: 2020-04-12 Received revised form: 2020-07-21 Accepted for publication: 2020-08-12 Convective Dry ng of Sewage Sludge Lay r in Through-flow Andraž Lipolt1 - Brane Širok2 - Marko Hočevar2 - Lovrenc Novak2 * 1 Petrol d.d., Slovenia 2 University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Drying of the sewage sludge layer was Investigated In a convective laboratory dryer at air temperatures of 65 ° C and 80 °C and air speeds of 0.53 m/s and 0.83 m/s. The sludge layer was formed by loading cylindrical extradâtes on a grate of 0.5 m x 0.5 m size. The drying air was directed through the layer, as typically encountered in Industrial belt dryers. Under such setup, the sludge layer structure and porosity significantly affect the air flow conditions and thus the drying rates. Shrinkage and cracking of the material during drying caused changes in the layer's porous structure, that affected the pressure drop and the drag force due to passing of air through the layer. The decreasing of drag force over time was modeled by a simple function that showed excellent agreement to the selected measured data. The sludge layer drying kinetics was determined by fitting the measured data to the most common drying models. Two models, the modified Nadhari and the Wang Singh model, were determined as most suitable for modeling of drying curves. The total drying time per kilogram of sludge was modeled as a function of drying air temperature, drying air velocity and initial sludge d/y matter content. The coefficient of determination (R2) of the model is 0.944. Total drying times between 43 minutes per kilogram and 76 minutes per kilogram of sludge were obtained for the investigated range of drying air conditions. Keywords: wastewater sludge, thin-layer drying, porous layer, drying kinetics, drag force Highlights • The sewage sludge d tying experiments were conduced by blowing hot air through the porous sludge layer of 0.25 m2 area. • The through-flow setup enabled representation of conditions in belt dryers and capturing of the effects of time-varying layer porosity. • The modified Nadhari and the Wang Singh model were determined as most suitable for modeling of drying curves. • The drag force due to passing of air through the layer was modeled by a simple decreasing function. • The total drying time per kilogram of sludge was modeled as a function of the sludge dry matter content and the drying air conditions. 0 INTRODUCTION The sludge generated in the wastewater treatment process is a very complex material with many of its properties depending on the wastewater characteristics and on the treatment technology. Wastewater characteristics vary by source, which can be any combination of domestic, municipal, also known as sewage, and industrial activities. Wastewater may contain physical, chemical and biological contaminants, therefore handling and disposal of sludge is of great concern. The most widely available options for sludge end-use in the EU are the agriculture utilization, the waste disposal sites, the land reclamation and restoration, the incineration and other novel uses [1]. Legal limitations concerning landfilling and agricultural reuse due to contaminants are increasing the role of sludge incineration in combination with recovery of nutrients, especially phosphorus [2] and [3], The sludge produced by wastewater treatment plants must generally be dehydrated before incineration, as it contains a high percentage of water. By using mechanical water extraction, it is typically possible to achieve sludge solid content up to 30 wt%. Further reduction of water content can be achieved by drying. Typical sludge drying methods used in the industry are convective drying, conductive drying and solar drying [4], The most common method is convective drying; compared to conductive drying it generally involves a simpler design of the drying system and shorter drying times [5], The convective drying systems are based on belt dryers, flash dryers, fluidized bed dryers and rotary dryers. The present study is focused on convective drying with specifics that are relevant to conditions in belt dryers. Fundamental aspects of sludge drying by the convective process have been examined by several authors. They presented the influential parameters of convective drying and concluded that the drying air temperature, its velocity and humidity are the most important parameters that affect the drying rates [6], Sludge composition has also been found to influence experimental results [7], which indicates the need to study drying behavior for sludges of different characteristics. Determination of optimum drying conditions prior to drying can improve energy efficiency of the process and avoid problems such as *Corr. Author's Address: University of Ljubljana, Faculty of Mechanical Engineering, Slovenia, Aškerčeva 6, Ljubljana, Slovenia, lovro.novak@fs.uni-lj.si 481 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 481-493 inhomogeneously dried material with the formation of a skin [8] and [9]. Most reported sludge drying studies are performed on small sludge samples with masses of a few grams in form of single cylindrical extrudates, balls or cakes [6] to [15]. Some studies also deal with larger sludge samples in form of a porous bed [16] to [20]. Study [5] compared drying of a small sludge sample (2.5 g) and sludge bed (250 g) and found significant differences between the two drying processes. Drying of sludge bed exhibited three typical drying phases: the adaptation phase, the constant drying rate phase and the falling drying rate phase, while no constant rate drying phase was detected during the drying of small sample. Size of sample has further implications on the drying process as some effects become evident only when drying larger quantities of sludge, typically in shape of porous layers. The layer characteristics and its transformation during drying, mainly due to shrinkage and cracking, influence the drying process. Moreover, in case of drying air being forced through the porous sludge layer, as typically encountered in industrial large-scale belt dryers, these effects are even more significant. Several researchers investigated the influence of sludge structural characteristics on heat and mass transfer phenomena [10], [14], and [19] to [22], with emphasis on changes in porous structure that are caused by shrinking and cracking of the material. Relationships between the sludge structure and the functional characteristics of drying, such as the energy efficiency of sludge drying, were presented. In most cases, the sludge obtained after the initial mechanical water removal processes is a soft, pasty and sticky wet material that is difficult to handle and dry. Due to the low stiffness of the resulting porous sludge structure, the layer deforms, which results in the reduction of free pores in the sludge and, consequently, the reduction of the transfer surface on which convective transfer of heat and water vapor takes place. The phenomenon was discussed by the authors in study [22], where adding of dry matter (sawdust) increased the stiffness of a porous substance on a laboratory scale, resulted in shorter drying times and higher energy efficiency of the process. However, it should be emphasized that technologies with the introduction of secondary raw materials such as sawdust [17] to [19] and [22] or pulverized coal [23] and [24] are associated with economic limitations arising from the costs of these materials and requirements for increased capacity of drying facilities. The aim of this work is to analyze the sewage sludge drying, where sludge is shaped as a porous layer and the drying air is blown through the layer. For this purpose, experiments were conducted on a laboratory drying facility, where sludge samples of several kilograms were dried in a layer of 0.5 m x 0.5 m size. The selected sample size enabled effects related to time-varying characteristics of the porous layer to manifest in large scale and thus affect the drying process. Contrary to most reported research, current study involved airflow being forced through the sludge layer. Under such setup the layer porosity on both local and integral scale significantly affects the air flow paths and thus the drying rates. The analysis of measured results was focused at evaluation of drying kinetics and at development of models for description of the sludge drying process at integral level. Since conditions in the presented study are closely related to the industrial belt dryers, the presented results could be scaled-up and applied to real large-scale dryers. 1 MATERIALS AND METHODS 2.1 Experimental Setup The sludge drying experiments were performed on a laboratory scale convective drying system as shown in Fig. 1. The central part of the system is the drying chamber with perforated grate, which holds the porous sludge layer. The grate area is 0.25 m2 and covers the entire chamber cross section. The circulation of air through the chamber is provided by a circulating fan located outside the drying chamber. The circulating drying air enters the drying chamber at the top, passes through the sludge layer and is extracted out of the chamber at the bottom. To maintain the desired drying air temperature, an electric heater is installed in the circulating air system. The heater has a nominal electrical power of 2 kW and is controlled by a temperature regulator. The extraction fan at the top of the chamber is used for partial extraction of moist air from the system and thus for supply of fresh air into the chamber. By adjusting the extraction flow, the relative humidity of the drying air in the chamber is controlled. Both the circulating and the extraction fan are controlled by variable frequency drives. By independently adjusting the circulating and extraction fan speeds and by controlling the heater, a wide range of drying air parameters above the layer can be achieved: air velocity up to 1 m/s, air temperature up to 90 °C (limited by the sludge weighing system) and air relative humidity up to 90 % (lower limit depends on ambient conditions and other parameters of the drying air). 482 Lipolt, A. - Širok, B. - Hočevar, M. - Novak, L. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 481-493 IR camera extraction fan circulation udg lave grate nr. ml e ectrica heater circu ation ar extraction Fig. 1. Laboratory convective dryer The convective dryer is equipped with a measurement system for monitoring of process variables, including the air and the sludge temperatures, the air relative humidity, the air static pressure, the air flow rate and the sludge mass. All temperatures were measured by the PtlOO resistive temperature detectors (RTD), accuracy class A. The air relative humidity (RH) was measured by the calibrated Honeywell HIH-4000 hygrometers (±3.5 % RH accuracy). The air static pressures were measured by the Endress+Hauser Deltabar PMD235 differential pressure transducers (±0.1 % accuracy). The drying air conditions were measured above the sludge layer (4 RTDs), below the sludge layer (4 RTDs) and in the extraction pipe at the top of the drying chamber (1 RTD and 1 hygrometer). Additionally, the ambient air conditions at the chamber inlet were also measured (1 RTD and 1 hygrometer). Downstream of the extraction fan an orifice plate was installed into the pipe to enable determination of the extracted air flow rate. The orifice plate and the pressure taps were manufactured according to the ISO standard [28] and connected to a differential pressure transducer. Air temperature was also measured at the orifice plate (1 RTD) to enable the air density determination and calculation of volume flow rate (accuracy ±3 %). Static pressure difference between the top and the bottom section of the drying chamber (pressure drop across the sludge and grate) was measured by a differential pressure transducer, connected to pressure taps at the chamber walls. The sludge condition was monitored by measuring its temperature and mass. The sludge surface temperature was measured by the FLIR T425 infra-red camera. The camera was mounted at the top of the drying chamber, where a ZnSe window was installed to enable observation of the sludge in IR spectrum. Online weighing of sludge mass was performed by a high temperature load cell Pavone Systems C2G1 HT with 20 kg capacity (combined error ±0.02 % of full scale output), which was installed into the drying chamber below the grate. Weighing of the sludge before and after drying was additionally performed on a Kern FKB 15K0.5A scale with linearity of 1.5 g. All sensors were connected to the Agilent 34970A data acquisition (DAQ) unit and were sampled in 20-second intervals over the entire drying time. Measured data was monitored and saved for further analysis by using the Benchlink Data Logger software running on a PC. The DAQ unit features high accuracy: ±(0.004 % reading + 0.004 % range) for DC voltage at 100 mV range and ±(0.002 % reading + 0.0007 % range) for DC voltage at 10 V range. Its contribution to combined measurement uncertainty in all measured variables is very low or insignificant. Nevertheless, it was included for calculation of the combined measurement uncertainty for online sludge weighing. According to the uncertainty propagation laws the sludge mass weighing uncertainty was calculated to be ±5.7 g. The drying air velocity was not measured directly but was calculated from the volume flow rate, which was determined on the basis of the fan performance curve, the measured pressure drop, the temperature and the relative humidity of the drying air. Additionally, the layer pressure drop vs. air velocity characteristic was determined prior to experiments. For this purpose, the layer pressure drop was varied by placing mats of different porosity on the grate. For each pressure drop, the air velocity was measured by a hot-film probe at 12 points just above the layer. The pressure drop was plotted as a function of the area-weighted average air velocity to obtain the layer resistance curve, which was then used for determination of the air velocity from the measured pressure drop. It lias to be noted, that the layer pressure drop represented a minor part of the total pressure drop across the air circulation system, therefore its variation had little influence on the fan performance and consequently the air velocity. Based on the typical fan performance curve measurement uncertainty and the presented layer resistance curve measurement method it is estimated that the air velocity was determined at ±5 % accuracy. 2.2 Experimental Procedure The sewage sludge used in the experiments was obtained from the sewage sludge belt dryer located in Ihan, Slovenia. The sludge was essentially a mixture Convectrue Drying of Sewage Sludge Layer in Through-flow 483 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 481-493 of sludges from different municipal wastewater treatment plants in Slovenia. The sludge was collected at the outlet of the die press system, which generates extradâtes of approximately 10 mm diameter. The extruded sludge was daily delivered to the location of the laboratory dryer in a sealed container in amounts that were sufficient for 2 drying tests. The initial moisture content of the sludge was determined from samples taken before each drying experiment and was on average 4.24 kg (water) per kg (dry matter) or 80.9 % water. The moisture content was determined by complete evaporation of water from the samples at 90 °C. The same method was employed to determine the final moisture content of sludge samples that were taken at the end of each drying test. Before each drying test, the grate was removed from the drying chamber and the laboratory dryer was preheated to the selected temperature level. The sludge was applied to the grate manually, where special attention was given to an even distribution of extradates and generation of an uniform layer. The initial sludge layer mass was 4.8 kg on average but could vary by as much as ±0.6 kg from test to test. The initial sludge layer thickness was approximately 40 ±5 mm. The extradate length after formation of the layer was approximately 50 mm, however, during manipulation some formation of larger sludge aggregates could not be avoided due to the soft and sticky nature of the material. Placement of extradates was predominantly parallel to the grate (horizontal) with random orientation. The prepared sludge layer was weighed on the laboratory scale and then inserted into the heated dryer. During drying, the circulation fan was periodically stopped for up to 30 seconds in order to eliminate the drag force from the weighed data according to the procedure, explained later in Eq. (1). The drying process was performed until no significant sludge mass variation could be detected. Then the grate with the dried sludge was removed from the dryer and again weighed on the laboratory scale. The experiment was carried out at different drying air conditions (temperature, velocity, relative humidity) to determine their effect on drying kinetics. Two levels were selected for each of the stated variables by considering the operating conditions, present in the full-scale convective dryer in Ihan, Slovenia. The low and high temperature level were selected to be 65 °C and 80 °C, respectively, and were maintained during drying by the temperature regulator. The high air velocity level was selected as the maximum possible (full fan speed) and was later determined to be on average 0.83 m/s, while the low 484 Lipolt, A. - Širok, B. - velocity level was selected by setting the fan rotational speed to 60 %, which corresponded to 0.53 m/s air velocity. The relative humidity (RH) of the drying air was not regulated during drying but instead two levels of fixed air extraction rates were selected. Under the low and high extraction rates, the average relative humidity of the air above the sludge layer was 20 % to 30 % and 40 % to 60 %, respectively. Table 1. Definition of operating conditions Test v [m/s] TTC] v ' ex RHave [%] 1-HHH 0.82 80 high 22 2-HHL 0.82 80 low 35 3-HLH* 0.84 65 high 32 4-HLL 0.84 65 low 60 5-LLL* 0.53 65 low 39 6-LLH 0.53 65 high 31 7-LHH 0.53 80 high 21 8-LHL 0.53 80 low 39 *complete drying curve not available from measured data In total 8 drying tests were required to cover all possible combinations of the selected process variable levels. In fact, more tests were performed, however, some had to be rejected due to incomplete or invalid data. After final analysis, 6 tests were identified to provide complete drying kinetics data and 2 additional tests were identified to provide limited data, sufficient only for the analysis of the total drying time. Further repetition of tests was not possible due to interruption in operation of the dryer in Ihan and consequent variability in supplied sludge properties. An overview of the performed tests is given in Table 1. Each test was given a designation consisting of a sequential number and a sequence of letters, indicating the set level for air velocity (v), temperature (T) and extraction rate ()t'v). where H and L indicate high and low level, respectively. It is evident from Table 1 that all possible combinations of the chosen process variables were measured, however, sludge mass data for 2 tests (3-HLH and 5-LLL) was insufficient for a complete drying kinetics analysis. It also has to be noted that the method of controlling the drying air RH by setting fixed air extraction rates was quite rough and could not maintain fixed RH levels or provide equal RH differences between the tests. 2.3 Processing of Measured Data The weighing data from the load cell in general included gravity forces due to the mass of sludge, the Hočevar, M. - Novak, L. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 481-493 mass of steel grate and the drag force due to passing of drying air through the porous sludge layer. After subtracting the mass of steel grate from the weighed data, the data still contained contributions from both the sludge mass and the drag force. The drag force was eliminated periodically by stopping the circulation fan for up to 30 seconds. The corresponding weighing at zero air flow provided the actual sludge mass. The drag force was calculated by subtracting the actual sludge mass from the weighed mass, obtained as soon as the circulation fan was back in normal operation. The delay between the actual sludge mass determination and the drag force calculation was typically 40 seconds, which corresponds to 2 data samples. drag,i (1) where md agj is the drag force (expressed as mass), mwe is the total weighed mass (without grate mass), msl is the sludge mass and i is the data sample number. The drying curves typically display either the moisture content (A") or the moisture ratio (MR) as a function of time. The moisture content is defined as the mass of water per mass of dry matter and was calculated according to Eq. (2). A' = - m., m., (2) where mw is the water mass, msl is the sludge mass and md, is the sludge dry matter. The dry matter was determined from sludge samples, taken before each drying test, by complete evaporation of water from the sample. The dimensionless moisture ratio can be calculated using the equation MR = (Xt-Xe)/(X0~Xe). The equilibrium moisture content, A'e. is often assumed to be zero, which simplifies the moisture ratio definition to [29]: MR = —, AV (3) where A' is the moisture content at given time and A"0 the initial moisture content. 2.4 Mathematical Modeling of Drying Curves The sludge mass as a function of time represents the fundamental characteristic of the drying process and is referred to as the drying curve. Characterization of drying kinetics typically involves modeling of drying curves by fitting measured values to various theoretical, semi-theoretical and empirical drying laws. Researchers [27] identified 67 different drying curve equations of various complexity for modeling of thin-layer drying, however, only 2 out of 390 reviewed references were dealing with sludge drying. Generally, modeling of sewage sludge convective drying by semi-theoretical or empirical equations has been attempted by several researchers [10], [15], [19], [21], [24] and [28] to [31], Studies demonstrated various models as the most appropriate for description of the drying curves, most probably due to the fact that sewage sludge composition and properties (hence drying kinetics) are highly variable. Furthermore, drying experiments were performed in different setups and under a wide range of conditions, which also affects the drying kinetics. Table 2. Used drying models Model name Model equation Parameters Newton [32] MR =exp (-kt) k Page [33] MR =exp (-kt») k, n Wang-Singh [34] MR = 1 +at. +b-f- a, b Logarithmic [35] MR =aexp(-kt)+b a, k, b Midilli [26] MR =aexp(-kt»)+bt a, k, n, b Mod. Midilli MR =exp (~kt»)+bt k, n, b Two-term [36] MR =a exp (-kt. )+b exp (-ti ) a, k, b, n Nadhari [37] MR =aexp(-kt»)+b a, k, n, b Modified Nadhari MR =a exp(-kt")+(\-a) a, k, n In the present study, drying curves were modeled by fitting measured data (moisture ratio versus time) to the most important semi-theoretical and empirical models (Table 2) that are widely used in the scientific literature to describe the kinetics of the drying process [27] and [38]. The original Midilli [26] and Nadhari [37] models were modified by considering the initial condition of MR = 1 at t = 0. which made one parameter in both of the models into constant, thus reducing the parameter count from 4 to 3. This simplification caused an insignificant reduction of the models' predictive power while making them more comparable to the other models that include 3 or less parameters. Regression analysis was used to obtain the parameters in each of the selected models. The analysis was performed in MS Excel by using the GRG nonlinear solver to minimize the sum of squared errors. The coefficient of determination (R2). reduced chi-square (■/}) and the root mean square error (RMK ) were calculated for each model in order to test their accuracy in reproducing the experimental data. The higher values of the coefficient of determination and the lower values of the reduced chi-square and RMK Convectrue Drying of Sewage Sludge Layer in Through-flow 485 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 481-493 were chosen for goodness of fit [38], These parameters were calculated as: %' = RAISE = N-n (4) (5) (6) where MR, MR expi is the experimental moisture ratio, -pre} the predicted moisture ratio, MR, expßv e the average experimental moisture ratio, N the number of data points and n the number of constants in the regression model. 3 RESULTS AND DISCUSSION 3.1 Modeling of Drying Curves The most significant measured quantities for a selected representative operating condition (8-LHL) are presented in Fig. 2. Thermal conditions are represented by the measured temperatures of the air below and above the sludge layer and by the calculated temperature difference across the layer ("'i layer" curve). The displayed air temperatures are in fact average values calculated from readings of 4 RTDs above and below the layer, respectively. The temperature of drying air was set to 80 °C in the presented operating condition however, it is evident that the initial air temperature was actually around 70 °C and reached the set value only after 60 minutes. The length of the air warming-up stage was a function of the maximum heater power and the operating conditions, primarily the set temperature. In operating conditions with lower set temperature (65 °C), the warming-up stage typically involved a starting temperature of 60 °C and lasted less than 20 minutes. Generally, the temperature difference across the layer was not affected by the initial air warming-up stage and instead showed a shorter transient period, lasting approximately 15 minutes. This period is characterized by the highest temperature differences and by fast increasing of the drying air relative humidity and reflects the initial sludge warming-up stage. From here on, the temperature difference is decreasing with an almost steady gradient which signifies the falling drying rate period. At around 135 minutes the falling rate of temperature and RH ("J layer" and "RII extraction" in Fig. 2a) is slightly increased, which could indicate a change in drying rate. However, such evident change in falling rate was not typical for other operating conditions. In contrast, the variation of pressure drop across the sludge layer (fi layer) as presented in Fig. 2b was typical for all operating conditions. The pressure drop decreased continuously with time. The highest pressure drop gradient and hence most pressure drop reduction occurred during the first half of drying. In the second half of drying the pressure drop gradually stabilized at an almost constant value. Periodical pressure drops of 0 Pa reflect stopping of the circulation fan for the purpose of sludge weighing. More details on the weighed sludge mass for the selected operating condition 8-LHL are given in Fig. 3. The actual sludge mass was recorded only when the circulation fan was stopped and consequently the drag force was eliminated. In the presented 7~above layer S 70 J" below layer RH ambient RH extraction d J" layer Fig. 2 O 50 100 150 ZOO 250 Time [min] Measured quantities for a representative operating condition; dp layer 100 150 Time [min] a) temperatures and relative humidities, and b) pressure drop 486 Lipolt, A. - Širok, B. - Hočevar, M. - Novak, L. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 481-493 ;>no:i too ) 3Ü00 1000 120 150 180 210 240 270 Time [min] Fig. 3. Weighed masses and sludge surface thermograms for operating condition 8-LHL operating condition, there were 19 instances of fan stopping for the purpose of weighing. The point at t = 0 corresponds to sludge weighing before drying and represents an additional sludge mass value. The drag force was determined by subtracting the sludge mass from the weighed total mass at normal fan operation (Eq. (1)). Fig. 3 also displays 3 thermograms of the sludge surface, taken at 3 different sludge weighing instances. The thermograms were used as an indication of the evenness of the drying rates across the sludge layer. Evenness or uniformity of drying is ensured by an uniform airflow through the sludge layer, which mostly depends on the layer structure and porosity distribution. High drying uniformity is desired as it is associated with higher drying efficiency, lower drying times and uniformly dried end-product. In the presented research, uniformity of sludge layer structure was monitored and assessed from the recorded thermograms as satisfactory, however, a more detailed analysis and objective quantification of these effects will be performed in scope of future research. The operating conditions with the most instances of fan stopping and hence the most measured values of actual sludge mass (4-HHL, 8-LHL) were expected to provide the definition of the best drying model equation. However, statistical criteria (RAIE , R2. X2) led to selection of a different equation for each operating condition. In the case of operating condition 4-HHL the best fit was obtained by the Wang-Singh model while in the case of operating condition 8-LHL the best fit was obtained by the modified Nadhari model. Summary of parameters related to fitting the two selected operating conditions is given in Table 3. All remaining measured operating conditions with less frequent sludge weighing were then fitted to the same models. The best fit of data as indicated by the statistical parameters was again obtained either for the Wang Singh or the Nadhari model, as evident in Table 4. In general, tests with higher drying rates were better modeled by the Nadhari model while tests with lower drying rates (drying at lower temperatures and flow rates) by the Wang Singh model, with the exception of operating condition 8-LHL, which was also better modeled by the Nadhari model. The modified Midilli model was statistically always very close to the modified Nadhari model, therefore its parameters were included in Table 4. Table 3. Parameters and statistics for drying models fitted to tests 4-HLL and 8-LHL Model param. Newton Page Wang Singh Logarithmic Two-term Nadhari a - - -0.00745 1.27274 - -14.024 1.1399 b_-_-_1.403E-05 0.006785 0.004348 15.04186_- k_0.009982 0.002569_-_-0.25278 1.144849 0.00399 0.0043 4-HLL n_-_1.290009_-_-_-0.00039 0.004282 1.1239 RAIE_0.04549 0.01688 0.00291 0.00846 0.00516 0.00719 0.00461 R2_0.9787_0.9971_0.9999_0.9993_0.9997_0.9995_0.9998 _X2_2.18E-03 3.18E-04 9.46E-06 8.49E-05 3.17E-05 6.54E-05 2.53E-05 a_-_-_-0.00688 1.278724_-_1.127077 1.0696 b_-_-_1.169E-05 0.007368 0.002077 -0.14182_- k_0.009215 0.001293_-_-0.19177 1.304418 0.010671 0.0023 8-LHL n_-_1.415025_-_-_-0.00016 1.210392 1.2720 RAIE_0.05692 0.01533 0.01694 0.01726 0.00817 0.03896 0.00713 R2_0.9724_0.9980_0.9976_0.9975_0.9994_0.9871_0.9994 _X2_3.42E-03 2.63E-04 3.21E-04 3.54E-04 7.92E-05 1.70E-03 8.09E-05 Convectrue Drying of Sewage Sludge Layer in Through-flow 487 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 481-493 Table 4. Parameters and statistics for drying models fitted to tests that are not included in Table 3 Param. Wang Singh Modified Midilli Modified Nadhari a -0.01087 - 1.0099 b 2.899E-05 0.002483 - k - 1.405839 0.0026 1-HHH n - -3.91 E-05 1.3967 RAIE 0.01745 0.00442 0.00417 R^ 0.99785 0.99986 0.99988 X2 4.26E-04 3.42E-05 3.04E-05 a -0.01043 - 0.9978 b 2.718E-05 0.001389 - k - 1.555674 0.0014 2-HHL n - 7.3E-06 1.5611 RAIE 0.02817 0.00729 0.00725 R^ 0.9956 0.99971 0.99971 X2 1.11 E—03 9.29E-05 9.20E-05 a —0.00686 - 1.0767 b 1.213E-05 0.003927 - k - 1.153867 0.0040 6-LLH n - -0.00019 1.1346 RAIE 0.00374 0.00434 0.00396 R^ 0.99989 0.99986 0.99988 X2 2.10E-05 3.76E-05 3.14E-05 a -0.00751 - 1.0443 b 1.447E-05 0.003548 - k - 1.208275 0.0036 7-LHH n - -0.00013 1.1940 RAIE 0.00407 0.00617 0.00575 R^ 0.99989 0.99976 0.99979 X2 2.76E-05 9.50E-05 8.28E-05 respectively. The measured conditons include weighed sludge mass before, during and after drying. Both the Wang Singh and the modified Nadhari models are plotted. Fig. 4. Weighed and fitted moisture ratio for operating condition 4-HLL Figs. 4 and 5 show measured and fitted moisture ratios for operating conditions 4-HLL and 8-LHL, Fig. 5. Weighed and fitted moisture ratio for operating condition 8-LHL In the case of operating condition 4-HLL both models produce very similar drying curves while in the case of operating condition 8-LHL, differences between the models are more evident. In the latter case the drying curve takes a different shape with initial adaptation stage clearly visible. Both models show some weaknesses at the end of drying, where they could predict negative drying rates (Nadhari) or change from falling to rising drying rates (going past the minimum of the Wang Singh curve, which is a parabola). Fig. 6 presents modeled moisture ratio for all operating conditions. The best model according to the statistical parameters (Tables 3 and 4) was selected for each operating condition. In total, 3 operating conditions were best fitted by the Wang Singh model and 3 by the modified Nahdari model. It is evident that combining high temperature and high air velocity resulted in significantly higher drying rates with drying kinetics, best described by the modified Nadhari model. Other combinations of velocity and temperature resulted in similar drying curves, that were mostly fitted by the Wang Singh model. Similar observations can be made on the basis of Krischer curves (Fig. 7) that present the drying rate as a function of moisture content. The initial drying rates (curves at the right-hand side of the diagram) are similar for all operating conditions, however, the curve shapes differ significantly. Operating conditions that were modeled by the Wang Singh model clearly show a continuously falling drying rate for the whole range of moisture rates. On the contrary, operating conditions that were 488 Lipolt, A. - Širok, B. - Hočevar, M. - Novak, L. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 481-493 best modeled by the modified Nadhari model exhibit a significant adaptation stage with rising drying rates that eventually level off and start to fall as the sludge gets dryer. The different drying behaviors, as reflected by the two types of drying curves, could be explained by the limiting mechanism for moisture evaporation. Léonard et al. [6] dried two sludges of different origin and concluded that the drying curve without adaptation phase reflects the presence of internal (diffusional) limitations from the beginning of drying, while the more typical drying curve reflects control of drying by external (sludge to air) limitations. It has to be noted though, that this explanation can not be proven in relation to the presented work and therefore has to be treated as hypothetical. - 1-HHH - 2-HHL -4-HLL 6-LLH V fr — •6- - \\ -- 7-LHH % V 11 8-LHL V 0 30 60 90 120 150 180 210 240 270 Time [min] Fig. 6. Modeled moisture ratio for all operating conditions: solid lines/markers according to modified Nadhari model, dashed lines/ no-fill markers according to Wang Singh model 0.05 0.0 0,5 1.0 1.5 2,0 2.5 3,0 3,5 4.0 4,S *[kg/kg] Fig. 7. Krischer curves for all operating conditions 3.2 Layer Porosity and Drag Force The presented drying curves enable basic qualitative assessment of relations between operating parameters and drying curves, but do not provide enough data for modeling of these relations. Furthermore, some uncertainty can be observed in the results which causes some of the drying curves to deviate from expectations. The unexpected variation of results could in large be attributed to non-controlled experiment parameters, where the sludge layer properties were identified as the most significant. The sludge layer was formed by loading pre-formed sludge extrudates on the drying grate by hand. Despite effort to produce sludge layers with consistent properties for all tests, there was inevitable variation of layer structure from test to test. In some cases, the extrudates were more easily deformed and formed aggregates that could not be manually separated while in some cases the extrudates were more stable and formed a more consistent porous layer. Consequent variation in layer structure and homogeneity could result in different courses of layer shrinkage and cracking, thus different evolution of layer porosity, air flow distribution and drying rate on a local scale. Current analysis could not be performed by considering local anomalies, but rather focused on measured integral (space-averaged) parameters. The most important integral layer properties that varied from test to test were identified to be its mass, initial moisture content and porosity. All these parameters were measured either directly (mass, initial moisture) or indirectly (porosity). Variation of porosity during drying was assessed by examining variation of pressure drop across the layer and the associated drag force. The two parameters were falling continuously during drying, as indicated in Figs. 2 and 3 for operating condition 8-LHL. The layer porosity was increasing with time due to shrinkage and cracking of the material. This phenomena is of great importance for the presented through-flow convective drying setup as it can significantly affect the contact area and the heat and mass transfer rates between the drying air and the dried material. Current results allow identification and qualitative assessment of the effects of variable porosity on drying, whereas a more detailed experimental treatment and analysis is foreseen in further studies. Fig. 8 presents drag force as a function of time for two operating conditions that have the most weighed data and therefore the best defined drying curves. The operating condition 4-HLL actually developed the highest drag force of all the tests whereas the operating Convectrue Drying of Sewage Sludge Layer in Through-flow 489 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 481-493 condition 8-LHL was in this respect a typical test at low air velocity. The measured data, represented by markers in Fig. 8, was fitted by a model that was determined by using the curve fitting tool ZunZun [39], which provides fitting of data with a wide range of functions and selection of best function according to the selected statistical criteria. In the presented case the best function was limited to 3-parameter equations and was defined as the function with the lowest RME value. The resulting model inform of Eq. (7) provided an excellent fit to data for the two measured operating conditions, as evident from the R2 values in Table 5. Fa = a-t +c. (7) The model parameters are clearly reflected in the curve shape: a+c represents the initial value (at i=0), b represents the curve gradient and c represents the final value (at / = inf.). The fitted parameter values are given in Table 5. 6000 5000 g 4000 G -o xt tU ¿z -aj 5 3000 2000 1000 K * \ i V A 4-HLL measured ■ 8-LHL measured ---4-HLL fitted -8-LHL fitted X A, L V. 300 0 100 200 Time [mtnj Fig. 8. Weighed and modeled drag force for operating conditions 4-HLL and 8-LHL Table 5. Drag force model parameters a b C R2 4-HLL 5351 -0.00372 464 0.9997 8-LHL 908 -0.00366 102 0.9963 The decreasing of drag force (or pressure drop through layer) during the course of drying could be associated with increasing of air flow rate. However, since changes in layer pressure drop had a minor effect on the fan performance, as explained in section Experimental setup, the increase in flow rate (or average air velocity) during drying was calculated to be maximum 3.5 %. 3.2 Total Drying Time In industrial sludge drying applications, it is typically required to remove most of water from the sludge, therefore modeling of total drying time under different drying air conditions was performed. For this purpose, the end of drying was defined to occur once the 3-minute averaged sludge mass gradient reached 0.5 g/min. In this instance, all drying curves practically leveled off and the sludge moisture ratios reached values of around 0.01. The selected definition of total drying time enabled evaluation of two additional tests (3-HLH, 5-LLL) that were not included in the preceding analysis due to insufficient data for full definition of their drying curves. Therefore, 8 operating conditions in total were available for the drying time modeling by the multiple regression analysis. The measured drying time data was fitted to the exponential and linear models by trying several combinations of independent variables and finally the exponential model was selected since it provided a better fit to data. The exponential model with n independent variables (x) and one dependent variable (v) is in general written as: v = A,fp (8) The most influential variables that defined the total drying time were identified to be the initial sludge mass (ms). the initial sludge moisture content or dry matter (DM), the drying air temperature (T) and velocity (v). Since influence of the initial sludge mass on the drying time was almost linear it was decided to define the drying time as specific time per mass of sludge, t„, = t/ms [min/kg]. By considering air velocity in m/s, air temperature in °C and dry matter in wt% and fitting coefficients according to the least squares method the following model was obtained: t =1043v-°-51ST-LU2DM-0S (9) Coefficient of determination (R2) for the model is 0.944. Plot of modeled vs. measured total drying time (Fig. 9) shows that the model adequately predicts the drying time within the investigated range of drying conditions. Exponents in the model are all negative which reflects negative correlation between the drying time and the drying air velocity and temperature and the initial sludge dry matter. Interestingly, the exponents values are very close to values of-0.5 and -1.0. By rounding the exponents to the mentioned values and adjusting the leading constant a simplified model is obtained: 490 Lipolt, A. - Širok, B. - Hočevar, M. - Novak, L. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 481-493 t =6i7v"°-Y"1mr1 = 617 s/v T DM' (10) Coefficient of determination (R2) for the above model is 0.935, which is comparable to the initial model (Eq. (9)). However, due to its mathematical simplicity, the Eq. (10) provides a clearer representation of relations between the involved variables. Similarities between the model and the physical laws, related to drying and flows through porous media, could not be determined. The coincidental simplicity of the model is linked to the selected system of units. This is especially evident in case of temperature, where changing units from degrees Celsius to Kelvin causes the exponent for temperature to take value far from -1.0. 80 70 s 60 o 50 40 /* * A * s r' * ✓ * / * $ A a <5 / * <£_ o best fit model * f x simple model -------equality line 40 50 60 70 Vkg] SO Measured time Fig. 9. Modeled vs. measured specific drying time As evident from Fig. 9, two points are located above the line of equality with a slightly higher deviation compared to other points. These two points correspond to tests 4-HLL and 8-LHL, which were reference tests with the most frequent sludge weighing. Incidentally, these two tests also developed higher initial drag forces compared to other tests. These specifics could act as accelerators of drying, since the total drying time in these two tests was actually shorter than predicted by the model. Definition of relations between the drag force (or pressure drop) and drying parameters was not possible based on current experimental data, therefore further experimental work and analysis is planned to investigate these effects in more detail. 4 CONCLUSIONS The sewage sludge drying experiments were conducted by forming the sludge sample as a porous layer and blowing the drying air in through-flow. Weighing of the sludge sample during drying was affected by the drag force, caused by passing of the air through the layer. The drag force variability between the tests was significant and was a function of sludge layer characteristics that could not be fully controlled. Nevertheless, the changing of drag force over time was modeled by a simple monotonic decreasing function that showed excellent agreement to the selected measured data. The sludge layer drying kinetics was determined by fitting the measured data to the most common drying models. Based on statistical criteria it was not possible to determine a single model, instead, two models, the modified Nadhari and the Wang Singh model, were determined as most suitable. The total drying time per kilogram of sludge was modeled as a function of the drying air temperature, the drying air velocity and the initial sludge dry matter content. Total drying times between 43 minutes and 76 minutes per kilogram of sludge were obtained for the investigated range of drying air conditions. The presented experiments were representative of the conditions in belt dryers, where sludge is dried in layers in through-flow. The results demonstrated the importance of time-dependent sludge layer characteristics, most specifically its porosity, that is reflected in the pressure drop through the layer and consequently the drag force. However, exact relations between the drying air parameters, the layer porosity and the drying kinetics could not be defined considering current data. Further experiments with upgraded measurement systems will be conducted in order to analyze changes of layer structure over time due to shrinkage and cracking and to relate these effects to the drying process. 5 ACKNOWLEDGEMENTS The authors acknowledge the financial support from the Slovenian Research Agency (research core Funding No. P2-0401). 6 REFERENCES [1] Fytili, D., Zabaniotou, A. (2008). Utilization of sewage sludge in EU application of old and new methods-A review. Renewable and Sustainable Energy Reviews, vol. 12, no. 1, p. 116-140, D0l:10.1016/j.rser.2006.05.014 [2] Yang, F„ Chen, J., Yang M„ Wang, X., Sun, Y., Xu, Y„ Qian, G. (2019). Phosphorus recovery from sewage sludge via incineration with chlorine-based additives. Waste Management, vol. 95, p. 644-651, D0l:10.1016/j. wasman.2019.06.029 Convectrue Drying of Sewage Sludge Layer in Through-flow 491 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 481-493 [3] Donatello, S., Tong, D., Cheeseman, C. R. (2010). Production of technical grade phosphoric acid from incinerator sewage sludge ash (ISSA). Waste Management, vol. 30, nos. 8-9, p. 1634-1642, D0l:10.1016/j.wasman.2010.04.009. [4] Bennamoun, L., Arlabosse, P., Léonard, A. (2013). Review on fundamental aspect of application of drying process to wastewater sludge. Renewable and Sustainable Energy Reviews, vol. 28, p. 29-43, D0I:10.1016/j.rser.2013.07.043. [5] Bennamoun, L., Fraikin, L., Li, J., Léonard, A. (2016). Forced Convective Drying of Wastewater Sludge with the Presentation of Exergy Analysis of the Dryer. Chemical Engineering Communications, vol. 203, no. 7, p. 855-860, D0I:10.1080/0 0986445.2015.1114475. [6] Léonard, A., Blacher, S., Marchot, P., Pirard, J. P., Crine, M. (2005). Convective drying of wastewater sludges: Influence of air temperature, superficial velocity, and humidity on the kinetics. Drying Technology, vol. 23, no. 8, p. 1667-1679, D0I:10.1081/DRT-200065082. [7] Léonard, A., Vandevenne, P., Salmon, T., Marchot, P., Crine, M. (2004). Wastewater sludge convective drying: Influence of sludge origin. Environmental Technology, vol. 25, no. 9, p. 1051-1057, D0I:10.1080/09593330.2004.9619398. [8] Tao, T., Peng, X. F., Lee, D. J. (2005). Structure of crack in thermally dried sludge cake. Drying Technology, vol. 23, no. 7, p. 1555-1568, D0I:10.1081/DRT-200063547. [9] Font, R., Gomez-Rico, M. F., Fullana, A. (2011). Skin effect in the heat and mass transfer model for sewage sludge drying. Separation and Purification Technology, vol. 77, no. 1, p. 146161, D0I:10.1016/j.seppur.2010.12.001. [10] Bennamoun, L., Fraikin, L., Léonard, A. (2014). Modeling and simulation of heat and mass transfer during convective drying of wastewater sludge with introduction of shrinkage phenomena. Drying Technology, vol. 32, no. 1, p. 13-22, D0I:10.1080/07373937.2013.807281. [11] Leonard, A., Blacher, S., Marchot, P., Crine, M. (2002). Use of X-ray microtomography to follow the convective heat drying of wastewater sludges. Drying Technology, vol. 20, nos. 4-5, p. 1053-1069, D0I:10.1081/DRT-120004013. [12] Fraikin, L., Herbreteau, B., Salmon, T., Nicol, F., Crine, M., Léonard, A. (2015). Use of an experimental design to characterize the convective drying behavior of different sludges. Drying Technology, vol. 33, no. 11, p. 1302-1308, D0I:10.1080/07373937.2015.1026979. [13] Kobayashi, N., Okada, K., Tachibana, Y., Kamiya, K., Ito, T., Ooki, H., Zhang, B., Suami, A., Itaya, Y. (2020). Drying behavior of sludge with drying accelerator. Drying Technology, vol. 38, nos. 1-2, p. 38-47, D0I:10.1080/07373937.2019.1605611. [14] Tao, T., Peng, X. F., Lee, D. J. (2005). Thermal drying of wastewater sludge: Change in drying area owing to volume shrinkage and crack development. Drying Technology, vol. 23, no. 3, p. 669-682, D0I:10.1081/DRT-200054164. [15] Danish, M., Jing, H., Pin, Z., Ziyang, L., Pansheng, Q. (2016). A new drying kinetic model for sewage sludge drying in presence of CaO and NaClO. Applied Thermal Engineering, vol. 106, no. 2016, p. 141-152, D0I:10.1016/j. applthermaleng.2016.05.191. [16] Léonard, A., Meneses, E., Le Trong, E., Salmon, T., Marchot, P., Toye, D., Crine, M. (2008). Influence of back mixing on the convective drying of residual sludges in a fixed bed. Water Research, vol. 42, nos. 10-11, p. 2671-2677, D0l:10.1016/j. watres.2008.01.020. [17] Li, J., Plougonven, E., Fraikin, L., Salmon, T., Toye, D., Nistajakis, E., Léonard, A. (2017). Influence of sawdust addition on drying of wastewater sludges: Comparison of structural characteristics. Drying Technology, vol. 35, no. 8, p. 925-932, D0I:10.1080/07373937.2016.1253020. [18] Li, J., Fraikin, L., Salmon, T., Bennamoun, L., Toye, D., Schreinemachers, R., Léonard, A. (2015). Investigation on convective drying of mixtures of sewage sludge and sawdust in a fixed bed. Drying Technology, vol. 33, no. 6, p. 704-712, D0I:10.1080/07373937.2014.982254. [19] Li, J., Bennamoun, L., Fraikin, L., Salmon, T., Toye, D., Schreinemachers, R., Léonard, A. (2014). Analysis of the shrinkage effect on mass transfer during convective drying of sawdust/sludge mixtures. Drying Technology, vol. 32, no. 14, p. 1706-1717, D0I:10.1080/07373937.2014.924136. [20] Li, J., Fraikin, L., Salmon, T., Plougonven, E., Toye, D., Léonard, A. (2016). Convective drying behavior of sawdust-sludge mixtures in a fixed bed. Drying Technology, vol. 34, no. 4, p. 395-402, D0I:10.1080/07373937.2015.1076835. [21] Bennamoun, L., Crine, M., Léonard, A. (2013). Convective drying of wastewater sludge: Introduction of shrinkage effect in mathematical modeling. Drying Technology, vol. 31, no. 6, p. 643-654, D0I:10.1080/07373937.2012.752743. [22] Li, J., Wu, C. W., Fraikin, L., Salmon, T., Toye, D., Nistajakis, E., Léonard, A. (2019). Convective drying of sawdust-sludge mixtures: The effect of the sludge matrix. Drying Technology, vol. 37, no. 7, p. 920-927, D0I:10.1080/07373937.2018.14 80027. [23] Lijuan, Z., Junhong, Y., Shanshan, W., Zhonghua, W. (2019). CO-drying characteristics of sticky sewage sludge preconditioned with biomass and coal. Drying Technology, vol. 37, p. 111, D0I:10.1080/07373937.2019.1692861. [24] Zhang, X. Y., Chen, M. Q., Huang, Y. W., Xue, F. (2016). Isothermal hot air drying behavior of municipal sewage sludge briquettes coupled with lignite additive. Fuel, vol. 171, p. 108115, D0I:10.1016/j.fuel.2015.12.052. [25] ISO 5167:2003. Measurement of fluid flow by means of pressure differential devices inserted in circular-cross section conduits running full: Orifice Plates. International Organization for Standardization, Geneva. [26] Midilli, A., Kucuk, H., Yapar, Z. (2002). A new model for single-layer drying. Drying Technology, vol. 20, no. 7, p. 1503-1513, D0I:10.1081/DRT-120005864. [27] Kucuk, H., Midilli, A., Kilic, A., Dincer, I. (2014). A Review on thin-layer drying-curve equations. Drying Technology, vol. 32, no. 7, p. 757-773, D0I:10.1080/07373937.2013.873047. [28] Zhang, X. Y., Chen, M. Q. (2016). A comparison of isothermal with nonisothermal drying kinetics of municipal sewage sludge. Journal of Thermal Analysis and Calorimetry, vol. 123, no. 1, p. 665-673, D0I:10.1007/s10973-015-4933-1. [29] Reyes, A., Eckholt, M., Troncoso, F., Efremov, G. (2004). Drying kinetics of sludge from a wastewater treatment plant. Drying Technology, vol. 22, no. 9, p. 2135-2150, D0I:10.1081/DRT-200034218. 492 Lipolt, A. - Širok, B. - Hočevar, M. - Novak, L. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 481-493 [30] Celma, A. R., Cuadros, F., Lopez-Rodriguez, F. (2012). Convective drying characteristics of sludge from treatment plants in tomato processing industries. Food and Bioproducts Processing, vol. 90, no. 2, p. 224-234, D0l:10.1016/j. fbp.2011.04.003. [31] Zhou, J., Zhang, R., Wang, X., Chen, S., Luo, A., Niu, D., Chai, X., Zhao, Y. (2017). NaHC03-enhanced sewage sludge thin-layer drying: Drying characteristics and kinetics. Drying Technology, vol. 35, no. 10, p. 1276-1287, D0I:10.1080/07373937.2016. 1249878. [32] Bruce, D. M. (1985). Exposed-layer barley drying: Three models fitted to new data up to 150°C. Journal of Agricultural Engineering Research, vol. 32, no. 4, p. 337-348, D0I:10.1016/0021-8634(85)90098-8. [33] Page, G.E. (1949). Factors Influencing the Maximum Rates of Air Drying Shelled Corn in Thin Layers. MSc thesis, Purdue University, West Lafayette. [34] Wang, C. Y., Singh, R. P. (1978). A single layer drying equation for rough rice, ASAE Paper, No. 78-3001. [35] Togrul, I. T., Pehlivan, D. (2002). Mathematical modelling of solar drying of apricots in thin layers. Journal of Food Engineering, vol. 55, no. 3, p. 209-216, D0I:10.1016/S0260-8774(02)00065-1. [36] Henderson, S.M. (1974). Progress in developing the thin layer drying equation. Transactions of the ASAE, vol. 17, no. 6, p. 1167-1168, D0I:10.13031/2013.37052. [37] Nadhari, W. N. A. W., Hashim, R., Sulaiman, 0., Jumhuri, N. (2014). Drying kinetics of oil palm trunk waste in control atmosphere and open air convection drying. International Journal of Heat and Mass Transfer, vol. 68, p. 14-20, D0I:10.1016/j.ijheatmasstransfer.2013.09.009. [38] Ertekin, C., Firat, M. Z. (2017). A comprehensive review of thin-layer drying models used in agricultural products. Critical Reviews in Food Science and Nutrition, vol. 57, no. 4, p. 701717, D0I:10.1080/10408398.2014.910493. [39] Philips, J. R. Zunzun, from http://zunzun.com, accessed on 2019-01-22. Convectrue Drying of Sewage Sludge Layer in Through-flow 493 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9,494-504 © 2020 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2020.6688 Original Scientific Paper Received for review: 2020-03-22 Received revised form: 2020-06-12 Accepted for publication: 2020-07-01 Microstructure and PropertjM odiflcation of Cold Spray d Coatings Using Different Grain Sizes of Cr3C2-25(Ni20Cr) Composite Powder Anna Trelka1 * - Wojciech Zorawski2 - Anna Goral1 1 Institute of Metallurgy and Materials Science, Polish Academy of Sciences, Poland 2Laser Processing Research Centre, Kielce University of Technology, Poland This article describes the microstructure and mechanical properties of Cr3C2-25(NI20Cr) cermet coatings cold-sprayed on Al 7075 alloy substrates. The study aimed to determine the effect of three different powder grain sizes and three different standoff distances on the microstructure (phase composition, surface topography volume fraction and distribution of the Cr3C2 phase in the NI20Cr matrix and porosity) and the mechanical properties (hardness and abrasive wear resistance) of the Cr3C2-25(NI20Cr) coatings. The examinations revealed that the thickness of the coatings obtained with the same duration and feed rate decreased with increasing powder grain size and standoff distance. The coatings produced from powders with smaller particles had higher hardness and more compact structures than the coatings obtained from powders where the particles were larger (more than 40 pm). The coatings produced from finer powders were also more resistant to abrasive wear. Keywords: Cr3C2-25(Ni20Cr) composite coating; cold spraying; microstructure; powder grain size; standoff distance; mechanical properties Highlights • The grain size and standoff distance significantly affect the microstructure and properties of Cr3C2-25(Ni20Cr) coatings, e.g., the content of the ceramic phase increases by over 40 % as the grain size Increases. • Cr3C2-25(NI20Cr) cold-sprayed coatings revealed a compact and uniform microstructure. • The thickness of the cold sprayed Cr3C2-25(Ni20Cr) coatings (obtained with the same duration and feed rate) decreased with increasing particle size of the powder and Increasing standoff distance. • The coatings made from the powder with the grain size of 9.5 pm to 55.3 pm had the highest hardness, while those made from the powder with the grain size of 40 pm to 55.3 pm the lowest. 0 INTRODUCTION Thermally sprayed composite coatings are common in industrial applications because of the wide variety of materials from which they can be made as well as numerous methods by which they can be deposited. In recent years, there has been a demand for innovative coatings to improve the properties of coated materials, which has led to a high development of research in the field of composite deposits [1] and [2], Coatings are designed mainly to improve the durability and enviromnental resistance of machine parts; they may also impart catalytic properties to them. Surface layers serve to enhance the mechanical properties of the underlying material, protecting it against deterioration due to corrosion, friction, abrasive wear, unfavourable temperatures and aggressive enviromnents [3] to [5], Cr3C2-NiCr coatings have been widely discussed in the literature because of their excellent resistance to abrasive wear [6] to [12], Composite coatings are used in enviromnents in which severe erosion and corrosion occur [13] to [15], Cr3C2-25(Ni20Cr) coatings are often produced using thermal spraying methods, especially high-velocity oxygen fuel [7], [8], [10] and plasma spraying [11], but, recently, also by cold spraying [6], [9] and [12], In some applications, cold spraying is superior to the other thermal spraying methods [16] because the gas temperature is lower than the melting point of the feedstock material, the powder particles do not undergo severe oxidation, and there is no phase change or growth of particles. In other words, problems typical of thermal spraying do not occur. The low porosity of cold-sprayed coatings and the absence of oxides in them may result in very high electrical conductivity as well as good corrosion resistance [17] and [18], Such coatings are formed as a consequence of severe plastic deformation of powder particles, adiabatic shear instability and mechanical interlocking [19], Powder particles are in the solidstate when they reach the substrate surface. Particles carried by the first stream remove all impurities together with the outer layer of oxides, which activates the substrate, making the next layers of particles adhere easily to it [20], A vital cold spray parameter is critical velocity, dependent on the type, morphology and size of the feedstock particles [21] and [22], As 494 *Corr. Author's Address: Institute of Metallurgy and Materials Science, Polish Academy of Sciences, ul. Reymonta 25,30-059 Krakow, Poland, a.trelka@imim.pl Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 494-504 pointed out by Sevillano et al. [23], the presence of small ceramic particles uniformly distributed in the metallic matrix of Al203-Ni composite coatings is beneficial, because the ceramic phase can act as a barrier against the propagation of cracks, while the metallic particles harden and strengthen under the impact of the subsequent layers, which may result in higher hardness and rigidity of the coatings. Celotto et al. [24] write that to deposit coatings successfully, various metallic particles 5 (im to 50 |im in size should be used. Powders with metallic particles smaller than 5 |im have poor flowability, and the reason is their lower weight and therefore insufficient kinetic energy to deform and adhere to the substrate. Particles larger than 50 |im have too much inertia and as such cannot be accelerated to a velocity higher than the critical velocity [24], Sova et al. [15], who analyse cold sprayed A1203-A1, SiC-Al, Al203-Cu and SiC-Cu composite coatings with different sizes of ceramic phase particles (19 |im to 25 |im and 135 |im to 141 |.un). indicate that, regardless of the type of material used, smaller particles ensured higher deposition efficiency of the feedstock and higher thickness of the coatings. Composite powders with larger particles contributed to coating erosion during deposition, and consequently lower efficiency of the deposition process. Another parameter that significantly affects the efficiency of the cold-spraying process and the quality of the coatings is the standoff distance between the nozzle and the substrate. As recommended in [25], the standoff distance should be less than or equal to the length of the non-scattered portion of the initial jet. The current challenge facing researchers is to develop new cold sprayed composite coatings with properties desirable in the automotive and aviation industries. Although some research has been carried out on cold sprayed composite coatings [5], [6], [9], [12] to [15], [23] and [26] to [29], no studies analysing the combined effect of the standoff distance and the powder grain size have been found. This study was undertaken to address the knowledge gap by investigating how the two parameters influenced the microstructure and the mechanical properties of the Cr3C2-25(Ni20Cr) coatings. 1 MATERIALS AND METHODS The cermet coatings were deposited on Al 7075 alloy substrates using Cr3C2-25(Ni20Cr) Diamalloy 3004 powder (Oerlikon Metco Europe GmbH, the Polish Division, Poznan, Poland) composed of balance Cr, 18.75 % Ni, 9.75 % C, 2.25 % other (max). The cold-spraying process was performed with an Impact Innovations 5/8 Cold Spray System, (Impact Innovations GmbH, Rattenkirchen, Germany) mounted on a Fanuc M-20iA robot arm (ZAP Robotyka, Oströw Wielkopolski, Poland). Table 1 shows the cold spray process parameters. Table 1. Parameters of the cold sprayed Cr3C2-25(NI20Cr) coatings Parameter Value Pressure [MPa] 4 Temperature [°C] 800 Standoff distance [mm] 10; 30; 50 Process gas n2 (DL: 9.5 to 40); Powder grain size [pm] (DO: 9.5 to 55.3); (DT: 40 to 55.3) Powder feeder rate [g/min] 95 ± 5 Number of layers 40 Speed of robot arm [m/s] 0.3 type convergent-divergent Nozzle length [mm] 161 designation INJECTOR-OUT 1-SIC (10202.00.0.01) The experiments were carried out for three different powder grain sizes and three different nozzle-substrate standoff distances: 10 mm, 30 mm, and 50 mm. The original commercially available Cr3C2-25(Ni20Cr) powder, denoted by DO, had a grain size ranging from 9.5 |im to 55.3 |im (Fig. 1), which was established using a Malvern Mastersizer analyser (HELOS Particle Size Analysis, Sympatec GmbH System-Partikel-Technik, Clausthal-Zellerfeld, Germany). The other powders, DL and DT, were obtained by dry sieving the DO powder through a 40 |im mesh screen. The DL powder had a grain size in the range of 9.5 |im to 40 |im. while the DT powder was composed of particles in the range of 40 (im to 55.3 |.un. The content of the ceramic phase in the powders determined using X-ray diffraction and image analysis (ImageJ) based on the powder cross-section microstructures (10 images in the case of each powder) obtained using SEM was: DL (77.6 ± 3.9 % vol.), DO (78.4 ± 2.7 % vol.), DT (80.3 ± 3.1 % vol.). Nine samples were analysed: DL1, DL3, DL5, DOl, D03, D05, DTI, DT3 and DT5. The symbols correspond to the powder grain size and the standoff distance. For example, DL1, DL3 and DL5 represent coatings produced from the DL powder at standoff distances of 10 mm, 30 mm and 50 mm, respectively. The coating thickness analysis was carried out for the optical microscope images using ImageJ software. Prior to the microscopic examinations, the specimens were cut, mounted in synthetic resin, ground with Microstructure and Property Modification of Cold Sprayed Coatings Using Different Grain Sizes of Cr3C2-25(Ni20Cr) Composite Powder 495 Strojnlskl vestnik - Journal of Mechanical Engineering 66(2020)9, 494-504 400 to 7000 grit paper, polished to 3 |im. 1 |im and 0.25 (rm diamond suspension finishes using a Presi MecaPol P320 (PRESI France, Eybens, France) setup, washed and dried. Fig. 1. Grain size distribution for the Cr3C2-25(Ni20Cr) (DO) powder The cross-sectional observations with a Leica DM IRM optical microscope (Leica Microsystems Wetzlar GmbH, Wetzlar, Germany) aimed to determine the thickness of the coatings and visually inspect their quality. The micro structure of the cold sprayed Cr3C2-25(Ni20Cr) coatings was characterised by scanning electron microscopy (FEI /Philips XL30, FEI Company, Oregon. United States). A Bruker D8 Discover diffractometer (Bruker AXS GmbH, Karlsruhe, Germany) with CoKa radiation (a wavelength of 1.7903 A) was used to study the phase compositions of the powders and the coatings. The volume fraction of the ceramic phase and the porosity in the coatings were measured using the ImageJ software. The surface topography analysis was performed with a Keyence VHX-7000 microscope (Keyence Corporation, Osaka, Japan). Three surfaces were selected and characterised on the basis of the respective 600 |im x 3000 |im topographical maps with five roughness parameters each. A CSM Instruments (Anton Paar GmbH, Graz, Austria) low-load Vickers hardness tester (HV0.3) was used to determine the coating hardness in the cross-section. In each case, the hardness result was an average of 12 measurements taken at different points. The coating hardness was measured in accordance with the standard ISO 6507-l:2018(en) [30], Hardness was measured according to the diagonal of the imprint. The abrasive wear tests were carried out using an ITEE T-07 tester (The Institute for Sustainable Technologies, Radom, Poland) (dry sand-rubber wheel) and loose abrasive A1203 particles 250 |im to 300 |im in size, at a flow rate of 250 g/min. a wheel (050 x 20) rotation speed of 200 rpm, and a load of 50 N. 2 RESULTS 2.1 Characterisation of the Cr3C2-25(Ni20Cr) powders Fig. 2 shows the morphologies of the Cr3C2-25(Ni20Cr) Diamalloy 3004 powders differing in the grain size (DL, DO and DT). All the powders consist of a metallic phase and a ceramic phase. The particles of the former are almost spherical or elongate in shape, while those of the latter are irregularly shaped with sharp and uneven edges. As can be seen from Figs. 2a and c, there is a significant difference in grain size between the DL powder (less than 40 |im) and the DT powder (more than 40 |im). Fig. 3 shows the phase analysis for the Cr3C2-25(Ni20Cr) powders. From the diffractogram, it is evident that there are two phases: the Cr3C2 ceramic phase (cell parameters: a = 5.533 A, b = 2.829 A and c = 11.472 A) and the metallic Cr0 25Ni0 75 phase (cell parameter: a = 3.552 A). The Cr0 25Ni0 75 phase contains 22.8 wt% Cr and 77.2 wt% Ni, according to the PDF 04-003-7001 card. 2.2 Characterisation of the Cr3C2-25(Ni20Cr) coatings 2.2.1 Thickness, Morphology and Microstructure The thickness measurement results obtained for the DL, DO, and DT coatings are compared in Fig. 4. The data indicate that the coating thickness decreases with Fig. 2. SEM-BSE (backscattered electrons) images depicting the morphologies of the Cr3C2-25(Ni20Cr) powders a) DL b) DO cj DT 496 Trelka, A. - Zörawsto, W. - G oral, A. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 494-504 increasing standoff distance. From Fig. 4, it is also clear that the greater the powder grain size, the thinner the coating. Fig. 3. X-ray diffraction patterns for the DL3, D03 and DT3 powders Fig. 4. Effect of the grain size and standoff distance on the thicknesses of the Cr3C2-25(Ni20Cr) coatings: DL (9.51urn to 40 pm), DO (9.5 pm to 55.3 pm), DT (40 pm to 55.3 pm) Fig. 5 shows the phase compositions of the DL3, D03, and DT3 coatings. The microstructural analysis of the coatings revealed that they also consisted of a ceramic phase (chromium carbide (Cr3C2)) and a metallic phase (Ni0 75Cr0 25). 70OO ^ MOO 2000 t 4 A AA > A - CrjC, 0 - PNZOCr dl3 ° ..._______ ...J A 003 A /l a dt3 a „ » 30 tB SO SO 70 SO 00 100 110 1JO 2 Theta [°] Fig. 5. X-ray diffraction patterns for the DL3, D03 and DT3 coatings A significant advantage of the cold-spray process, confirmed through the X-ray diffraction analysis, is that no change in the phase composition occurs during deposition. Cold-spraying does not cause oxidation of the material or formation of new phases, as is the case with other thermal spraying methods. The surface morphologies of the coatings tested are shown in Fig. 6. The arrows indicate the ceramic phase. The analysis of the DL coatings revealed the presence of the Ni20Cr phase with embedded ceramic phase Cr3C2 particles. Visual observation revealed that the area of the metallic phase was the smallest when the standoff distance was 10 mm (Fig. 6a). In the coatings produced at 30 mm (DL3) or 50 mm (DL5), the area of the metallic phase was more extensive, while the ceramic phase was highly fragmented (Figs. 6b and c, respectively). The examination of the DO coatings made from the original powder shows that the largest zone of the metallic phase was obtained for D03, sprayed at a distance of 30 mm (Fig. 6e) and that it was a result of high plastic deformation. As can be seen from the morphologies of the coatings presented in Fig. 6, the highest volume fraction of the Cr3C2 ceramic phase at the surface is observed for DT deposits. During the deposition process, the brittle ceramic particles were crushed and unevenly distributed in the metallic phase, giving the impression that powder with smaller particles was used. The higher amount of ceramic particles at the surface was responsible for the higher roughness of the DT coatings. The surface topography parameters of the coatings deposited at the same standoff distance of 30 mm but differing in initial powder grain size are presented in Fig. 7 and Table 2. The surface roughness parameter, i.e., the arithmetic mean height (Sj. increased with increasing powder grain size. The difference in roughness between the DL3 and D03 coatings (Figs. 7a and b, respectively) was small, with the difference in maximum height (S:) being 10.1 |im. In the case of the DT3 coating (Fig. 7c), the parameter was more than twice as high as those obtained for DL3 orD03. All the parameters indicate that the grain size is the factor affecting the surface topography. The texture aspect ratio of the surface (Str) decreased with increasing powder grain size. The arithmetic mean peak curvature (Spc) represents the arithmetic mean of the principal curvature of the peaks on the surface. The lower its value, the more rounded the shapes of the points of contact. The developed interfacial area ratio (Sd) is expressed as the percentage of the definition area's additional surface area contributed by the texture as compared to the planar definition area. Microstructure and Property Modification of Cold Sprayed Coatings Using Different Grain Sizes of Cr3C2-25(Ni20Cr) Composite Powder 497 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 494-504 Fig. 6. SEM-BSE images showing the surface morphologies of the Cr3C2-25(Ni20Cr) coatings from powders differing in grain size cold sprayed at three different standoff distances: a) DL1, b) DL3, c) DL5, d) D01, e) D03, f) D05, g) DTI, h) DT3 i) DT5; arrows indicating Cr3C2 Table 2. Area surface topography parameters determined for the Cr3C2-25(Ni20Cr) coatings Coating surface Parameters * Height Spatial Feature Hybrid Sa [um] S.tum] Str Svc [1/mm] DL3 5.1 55.8 1.0 958.5 0.5 D 03 5.7 65.9 0.9 999.3 0.6 DT3 6.7 138.4 0.7 1143.8 0.9 * Sa arithmetic mean height, S: maximum height, Str texture aspect ratio, Spc arithmetic mean peak curvature, Sd developed interfacial area ratio dark grey, and the metallic Ni20Cr phase, visible as light grey. The substrate is black. The DL coatings have small ceramic phase particles evenly distributed in the Ni20Cr alloy matrix, which is particularly visible for the DL3 coating. The Cr3C2 phase in the DL3 and DL5 coatings is in the form of strands. The most uniform micro structure is that of the DL3 coating. The microstructures of the DO coatings seem more compact than those of the DL and DT coatings. The ceramic particles are distributed in the metallic matrix in the form of strands. The analysis of The parameters Spc and Sd increase with increasing grain size (Table 2). Fig. 8 shows cross-sectional views of the different cold sprayed coatings. The deposits are composed of the ceramic Cr3C2 phase, marked in Fig. 7. Surface topographies of the Cr3C2-25(Ni20Cr) coatings from powders differing in grain size produced at a standoff distance of 30 mm a) DL3, b) D03, c) DT3 um LI m mm 1000 [j m ÎOOD L J I um 1000 um um ■ um 498 Trelka, A. - Zörawsto, W. - Gora/, A. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 494-504 Fig. 8. SEM-BSE Images depicting the cross-sectional mlcrostructures of the Cr3C2-25(NI20Cr) coatings from powders differing In grain size cold sprayed at three different standoff distances: a) DL1, b) DL3, c) DL5, d) D01, e) D03, f) D05, g) DTI, h) DT3, i) DT5 the average porosity of all the coatings (DL, DO and DT) shows that the DO coatings exhibit the lowest porosity, ranging from 1.7 ± 0.4 % vol. to 1.9 ± 0.4 % vol. (Fig. 9). When single coatings are considered, the DL3 coating has the lowest porosity (1.6 % vol.). The volume fraction of the ceramic phase in the DT coatings, varying from 38.8 ± 2.8 % vol. to 39.7 ± 0.9 % vol., is much higher than those in the DL, being in the range of 26.4 ± 0.8 % vol. to 28.8 ± 1.9 % vol. or DO, ranging between 31.8 ± 1.8 % vol. and 33.8 ± 1.6 % vol. However, high fragmentation of ceramic particles is visible in the region closest to the substrate. In the case of the DT5 coating, the fragmentation affected porosity, which was about 1 % vol. higher than in the other coatings. The microstructures of the coatings are less compact, which is associated with the cracking and fragmentation of the ceramic phase particles (Figs. 10b and c). The porosity observed between the crushed Cr3C2 particles could have been the result of particle crumbling during sample preparation for scanning electron microscopy (Fig. 10). Fig. 9. Porosities (% vol.) of the Cr3C2-25(NI20Cr) coatings The results provided in Fig. 11 indicate that the volume fraction of the ceramic phase decreased with increasing standoff distance; this, however, is not true for the D05 coating. Higher volume fractions of the ceramic phase were observed when the feedstock powder particles were larger. The highest volume fraction of the ceramic phase was reported for the DT coatings produced from the powder with particles greater than 40 |im. while the lowest for the DL Microstructure and Property Modification of Cold Sprayed Coatings Using Different Grain Sizes of Cr3C2-25(Ni20Cr) Composite Powder 499 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 494-504 Fig. 10. SEM-BSE cross-sectional views of the coating microstructures: a) DL3, b) DO3, c) DT3; CP - porosity caused by the crumbling of finer ceramic particles coatings. The difference between them was 10 % vol., which may have resulted from the bouncing of small ceramic particles off the surface during the cold spray process. The effect of the standoff distance on the volume fraction of the ceramic phase in the coatings was negligible (about 2 % vol.). Fig. 9 shows that the highest porosity occurred at a standoff distance of 50 mm when powder with a grain size greater than 40 ^m was used. Fig. 11. Cr3C2 particles embedded in the Cr3C2-25(Ni20Cr) coatings (% vol.) 2.2.2 Mechanical and Tribological Properties The hardness measurement results show that the DO coatings, obtained from the original powder, have the highest hardness, while the DT coatings the lowest (Fig. 12). The hardness of the DL coatings was slightly lower than that of the DO coatings. The low hardness of the DT coatings was a consequence of their high microstructural porosity related to the crumbling of the high amount of ceramic phase embedded in the Ni20Cr matrix. The high values of the standard deviation obtained for the DT specimens indicate a non-uniform structure resulting from the high porosity of the coatings. No significant relationship was observed between the coating hardness and the standoff distance. The difference between the highest and lowest values of hardness obtained for each standoff distance was the same, irrespective of the grain size; it was about 20 HV0.3. Fig. 13 shows the load-displacement curves and the hardness indentations obtained for the DO3 and DT3 specimens. Fig. 12. Hardness of the Cr3C2-25(Ni20Cr) coatings Fig. 13. a) A load-displacement curve and b) hardness indentations for the DO3 and c) DT3 Cr3C2-25(Ni20Cr) coatings The indentation in the DO3 coating was smaller than that in the DT3 (Figs. 13b and c, respectively), which suggests its higher hardness. The porosity observed around the irregularly shaped indentation in the DT3 specimen (Fig. 13 c) had a considerable effect on the coating hardness. From the load-displacement curves in Fig. 13 a, it is also clear that the hardness of 500 Trelka, A. - Zörawsto, W. - Gora/, A. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 494-504 the DT3 coating is lower. At a load of 2.942 N, the indenter displacement was higher than 4.5 |im for the D03 specimen, while for DT3 it was higher than 5 (im. Fig. 14 shows the mass loss for the DL3, D03 and DT3 coatings after the abrasive wear tests. The DL3 coating containing the lowest volume fraction of the ceramic phase in comparison with the D03 and DT3 deposits revealed the most considerable mass loss. The D03 and DT3 coatings had an identical mass loss after 30 min of the abrasive wear tests (125 mg), which was 29 % lower than that reported for DL3 (161 mg). It was also found that the wear time for DT3 was shorter (30 min) than for the other coatings (40 min). Compared to DL3 or D03, this coating was about 37.7 % thinner, which was the reason for its shorter wear time. The microstructures of the DL3 and D03 coatings studied after the abrasive wear tests by optical and scanning electron microscopes are shown in Fig. 15. The differences in the appearance of wear track between the coatings after the tests with loose abrasive corundum (A1203) were insignificant. The wear track reported for D03 (Fig. 15b) was more uneven than that observed in DL3 (Fig. 15a), which was probably due to the greater grain size and the higher content of ceramic particles, which broke off during friction. 250 O 10 20 30 40 SO Time [min] Fig. 14. Mass loss during the abrasive wear test 3 DISCUSSION This study provides new insights into cold-sprayed composite coatings, especially the effects of the powder grain size and the standoff distance on the microstructure and properties of the Cr3C2-25(Ni20Cr) cermet coatings. The comparison of the diffraction patterns of the feedstock powders with those of the cold sprayed coatings revealed several differences. The relative intensity of the chromium carbide peaks was higher for the powders than for the coatings (Figs. 3 and Fig. 15. Microstructures of the Cr3C2-25(Ni20Cr) coatings after the abrasive weartests: a) and c) DL3, b) and d) D03 performed using optical (a, b; magnification lOOx) and scanning electron microscope (c, d; magnification 2000x) Microstructure and Property Modification of Cold Sprayed Coatings Using Different Grain Sizes of Cr3C2-25(Ni20Cr) Composite Powder 501 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 494-504 5, respectively). This suggests that the content of chromium carbide decreased during the deposition process, which is confirmed by the results in Fig. 11. Much fewer Cr3C2 particles were found in the coating structure, and this was because they did not deform plastically during the cold spraying process; instead, they bounced off the substrate surface. The broader diffraction peaks for the chromium-nickel phase indicate a decrease in the grain size. Similar relationships are reported by Wolfe and Eden [12] for Cr3C2-NiCr. The microstructural analysis of the Cr3C2-25(Ni20Cr) coatings in cross-section (Fig. 8) reveals that the ceramic particles in the DL coatings are smaller and distributed evenly, while those in the DT coatings are larger and arranged in the strands-like elongated forms. Small ceramic particles uniformly embedded in a coating are considered beneficial because they can act as a barrier against crack propagation. They also play important roles of peening and roughing of the layers during deposition, modifying the coating growth and, as a consequence, the coating properties [29], Metallic particles, in contrast, harden and strengthen under the impact of the subsequent layers, which may cause higher hardness and rigidity of the coatings [23], Ni20Cr particles are plastically deformed and bonded both with the substrate and each other as a result of mechanical interlocking and adiabatic shear instability. The significant loss of the ceramic phase in relation to the metallic phase is explained by Fernandez and Jodoin [9], who argue that, when ceramic particles crack upon impact, some are lost as a result of weak bonding between them and are blown away by the gas stream. They also report an increase in the volume fraction of the NiCr phase. The coating thickness is strictly related to the deposition efficiency of the process [31], Sova et al. [15] state that, for coatings made from A1203-A1, SiC-Al, A1203-Cu and SiC-Cu powders, the deposition efficiency strongly depends on the size of powder grains and is higher for smaller particles. This study confirmed this relationship, i.e., that the thickness of the Cr3C2-25(Ni20Cr) coatings decreased with increasing powder grain size. The DL coatings, which were made from the smallest particles, were the thickest. Large ceramic particles phase are brittle, which hinders the bonding of the sprayed material with the previously deposited layer [9], The analysis shows that regardless of the size of the powder grains, an increase in the standoff distance causes a decrease in the thickness of the coatings (Fig. 4) and a slight decrease in the volume fraction of the Cr3C2 ceramic phase, the exception being the D05 coating sprayed at a distance of 50 mm (Fig. 11). The difference in thickness between the coatings obtained at the spraying distance of 10 mm and those produced at 50 mm was about 100 |.un. no matter how large the powder particles were (DL, DO or DT). In cold spraying, longer standoff distances result in thinner coatings because of the bouncing of particles off the surface during the formation of subsequent coating layers, which may be due to an increase in the particle momentum, being a result of a longer standoff distance. Smaller particles achieve higher velocities than larger ones, and since powders are a mixture of particles of different sizes, some of the powder is deposited, while the remainder bounces off [16], The greater size of grains in the DT3 coatings resulted in a surface roughness of at least 18 % higher than for the other coatings. Higher roughness can affect coatings both positively and negatively. While the subsequent layers deposited on the substrate bond easily, the significant fragmentation of ceramic particles increases the undesirable porosity [9]. The high value of the roughness parameter Sa observed for the DT3 coating was due to the higher volume fraction of ceramic particles. The microstructural examination showed that porosity occurred mainly where the particles of the ceramic phase crumbled (Fig. 10c). The preparation of the specimen cross-sections for the SEM examinations also increased the porosity of the composite materials. During that process, very small fragments of the ceramic particles broke off and fell out. All the coatings studied were characterised by low porosity (less than 3 % vol.). Sevillano et al. [23] observed a similar relationship for Ni-Al203 coatings, which had low porosity (about 3 % vol.) with pores found mainly in cracks in the ceramic phase. As for the mechanical properties, the DO coatings exhibited the highest hardness, ranging between 615 ± 66 and 635 ± 42 HV0.3, while the DT ones had the lowest, varying from 552 ± 100 to 575 ± 49 HV0.3 (Fig. 12). According to the literature [26], the higher volume fraction of the ceramic phase, the higher hardness of the Ni-WC composite coatings. No such observations were made in this study. The reason for that was the cracking of large Cr3C2 particles during the spraying process, which resulted in lower hardness of the coatings. In the Cr3C2-25(Ni20Cr) coatings produced from the DO and DL powders, lower volume fractions of the ceramic phase did not lead to lower hardness because the microstructures of the coatings were compact (Fig. 8). They showed a small difference in hardness of about 4.5 % (Fig. 12). The coatings made from the finest powder showed higher mass loss during the abrasive wear tests (Fig. 14). 502 Trelka, A. - Zörawsto, W. - Gora/, A. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 494-504 This was attributed to the lower volume fraction and smaller size of the ceramic phase particles (Fig. 11). From the results concerning the coating thickness, hardness, abrasion resistance, porosity and content of ceramic phase, it can be concluded that the Cr3C2-25(Ni20Cr) coatings produced from the original powder (DO) sprayed at a distance of 30 mrn have the best properties. They are compact and reveal the high hardness (635 ± 42 HV0.3), a high volume fraction of the ceramic phase (31.8 ± 1.8 % vol.), low porosity (1.7 ± 0.4 % vol.), and considerable thickness (676 ±21 (un), which confirms good relative deposition efficiency. The coatings made from smaller (9.5 |im to 40 |.un) powder grains at different standoff distances (DL1, DL3, DL5) had higher thickness, but they showed lower hardness. The DT3 coating, in contrast, had a high content of the ceramic phase but low hardness and highest porosity; it was also the thinnest of all the deposits. The experimental data indicate that the desired mechanical properties of Cr3C2-25(Ni20Cr) coatings can be achieved economically by selecting the optimal parameters, i.e., powder grain size ranging from 9.5 (un to 55.3 |im and a standoff distance of 30 mm. 4 CONCLUSIONS The purpose of this study was to determine how the size of the Cr3C2-25(Ni20Cr) powder particles and the nozzle-substrate standoff distance (10 mm, 30 mm or 50 mm) affect the micro structure and mechanical properties of the Cr3C2-25(Ni20Cr) cermet coatings cold sprayed on the A1 7075 substrate. The following conclusions have been drawn from the study results: 1. All the Cr3C2-25(Ni20Cr) coatings had a compact microstructure; however, those that were produced from the powder with the largest particles had higher porosity. 2. The cold-sprayed coatings obtained from the powder with particles greater than 40 (un (DT) had the highest volume fraction of the Cr3C2 ceramic phase. The volume fraction of the Cr3C2 ceramic phase was the lowest for the coatings produced from powder with particles smaller than 40 (un (DL). The standoff distance had little effect on the content of the ceramic phase in the coatings. 3. The thickness of the cold-sprayed Cr3C2-25(Ni20Cr) coatings decreased with increasing powder grain size and increasing standoff distance. 4. The coatings made from the powder with the original grain size (DO: 9.5 (un to 55.3 (un) had the highest hardness. The hardness of the DT coatings made of the powder with the grain size 40 (un to 55.3 (un, in contrast, was the lowest. 5. The cold-sprayed coatings obtained from the powders with larger particles (DO and DT) had higher resistance to abrasive wear when tested with loose abrasive particles than the coatings obtained from the finest powder. 6. The Cr3C2-25(Ni20Cr) coatings produced from the DO powder at a standoff distance of 30 mm had the best properties, i.e., high thickness, high hardness, high resistance to abrasive wear, a relatively high volume fraction of the Cr3C2 phase and low porosity. 5 ACKNOWLEDGEMENTS This work was supported by the National Science Centre, Poland (Project No 2017/25/B/ST8/ 02228). 6 REFERENCES [1] Yi, H, Liu, X, Che, J, Liang, G. (2020). Thermochemical compatibility between La2(Cel-xZrx)207 and 4 mol% Y203 stabilized zirconia after high temperature heat treatment. Ceramics International, vol. 46, no. 4, p. 4142-4147, D0l:10.1016/j.ceramint.2019.10.130 [2] Goral, A, Zorawski, W, Lityriska-Dobrzyriska, L. (2014). Study of the microstructure of plasma sprayed coatings obtained from AI203-13Ti02 nanostructured and conventional powders. Materials Characterization, vol. 96, p. 234-240, D0l:10.1016/j.matchar.2014.08.016 [3] Pawtowski, L. (2008). TheScienceandEngineeringof Thermal Spray Coatings. 2nd ed. John Wiley & Sons, Chichister. [4] Goral, A., Lityriska-Dobrzyriska, L., Zorawski, W., Berent, K., Wojewoda-Budka, J. (2013). Microstructure of AI203-13Ti02 coatings deposited from nanoparticles by plasma spraying Archives of Metallurgy and Materials, vol. 58, no. 2, p. 335339, D0l:10.2478/vl0172-012-0194-l [5] Kusiriski, J., Kac, S, Kowalski, K, Dosta, S, Georgiou, E.P, Garcia-Forgas, J., Matteazzi, P. (2018). Microstructure and properties of TiC/Ti coatings deposited by the supersonic cold gas spray technique. Archives of Metallurgy and Materials, vol. 63, no 2, p 867-873, D0l:10.24425/122416 [6] Wolfe, D.E., Eden, T.J., Potter, J.K., Jaroh, A.P (2006). Investigation and characterization of Cr3C2-based wear-resistant coatings applied by the cold spray process. Journal of Thermal Spray Technology, vol. 15, p. 400-412, D0l:10.1361/105996306X124400 [7] Matthews, S, James, B, Hyland, M. (2009). The role of microstructure in the mechanism of high velocity erosion of Cr3C2-NiCr thermal spray coatings: Part 1 - As-sprayed coatings. Surface and Coatings Technology, vol. 203, no. 8, p. 10861093, D0l:10.1016/j.surfcoat.2008.10.005 [8] Roy, M, Pauschitz, A, Bernardi, J, Koch, T, Franek, F. (2006). Microstructure and mechanical properties of Microstructure and Property Modification of Cold Sprayed Coatings Using Different Grain Sizes of Cr3C2-25(Ni20Cr) Composite Powder 503 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 494-504 HVOF sprayed nanocrystalline Cr3C2-25(Ni20Cr) coating. Journal of Thermal Spray Technology, Vol. 15, p. 372-381, DOI:10.1361/105996306X124374. [9] Fernandez, R., Jodoin, B. (2017). Effect of particle morphology on cold spray deposition of chromium carbide-nickel chromium cermet powders. Journal of Thermal Spray Technology, vol. 26, p. 1356-1380, DOI:10.1007/s11666-017-0580-3. [10] Zavareh, M.A., Sarhan, A.A.D.M., Razak, B.B., Basirun, W.J. (2015). The tribological and electrochemical behaviour of HVOF-sprayed Cr3C2-NiCr ceramic coating on carbon steel. Ceramics International, vol. 41, no. 4, p. 5387-5396, D0I:10.1016/j.ceramint.2014.12.102. [11] Zhang, L., Hou, J.B. (2013). Study of microstructure and phase of plasma sprayed Cr3C2-NiCr coating before and after the sparking plasma sintering. Physics Procedía, vol. 50, p. 293296, D0I:10.1016/j.phpro.2013.11.047. [12] Wolfe, D., Eden, T. (2007). Cold spray particle deposition for improved wear resistance. Champagne, V.K. (ed). The Cold Spray Materials Deposition Process. Fundamentals and Applications. Woodhead Publishing, Cambridge, p. 264-301, DOI:10.1533/9781845693787.3.264. [13] Wang, Y., Normand, B., Mary, N., Yu, M., Liao, H. (2017). Effects of ceramic particle size on microstructure and the corrosion behavior of cold sprayed SiCp/Al 5056 composite coatings. Surface and Coatings Technology, vol. 315, p. 314325, D0I:10.1016/j.surfcoat.2017.02.047. [14] Silva, F.S., Cinca, N., Dosta, S., Cano, I.G., Guilemany, J.M., Benedetti, A.V. (2017). Cold gas spray coatings: basic principles, corrosion protection and applications. Eclética Química Journal, vol. 42, no. 1, p. 9-32, D0I:10.26850/1678-4618eqj.v42.1.2017.p09-32. [15] Sova, A., Papyrin, A., Smurov, I. (2009). Influence of ceramic powder size on process of cermet coating formation by cold spray. Journal of Thermal Spray Technology, vol. 18, no. 4, p. 663-641, D0I:10.1007/s11666-009-9359-5. [16] Champagne, V.K. (2007). Introduction, Champagne, V.K. (ed.) The Cold Spray Materials Deposition Process. Fundamentals and Applications. Woodhead Publishing, Cambridge, p. 1-7, D0I:10.1533/9781845693787.1. [17] Schmidt, T., Gartner, F., Stoltenhoff, T., Kreye, H., Assadi, H. (2005). High velocity impact phenomena and coating quality in cold spraying. Proceedings of the International Thermal Spray Conference, p. 232-238. [18] Klassen, T., Kliemann, J.O., Onizawa, K., Donner, K., Gutzmann, H., Binder, K., Schmidt, T., Gartner, F., Kreye, H. (2009). Cold spraying - new developments and application potential. 8th Kolloquium HVOF Spraying, 9-16. [19] Grujicic, M. (2007). Particle/substrate interaction in the cold-spray bonding process. Champagne, V.K. (ed). The Cold Spray Materials Deposition Process. Fundamentals and Applications. Woodhead Publishing Limited, Cambridge, p. 148-177, D0I:10.1533/9781845693787.2.148. [20] ASM International Handbook Committee (2004). Thermal Spray Technology ASM Handbook. ASM International, Materials Park. [21] Schmidt, T., Gärtner, F., Assadi, H., Kreye, H. (2006). Development of a generalized parameter window for cold spray deposition. Acta Materialia, vol. 54, no. 3, p. 729-742, D0l:10.1016/j.actamat.2005.10.005. [22] Góral, A., Zórawski, W., Czaja, P., Lityñska-Dobrzyñska, L., Makrenek, M., Kowalski, S. (2019). Effect of powder morphology on the microstructure and properties of cold sprayed Ni coatings. International Journal of Materials Research, vol. 110, no. 1, p. 49-59, D0I:10.3139/146.111698. [23] Sevillano, F., Poza, P., Munez, C.J., Vezzu, S., Rech, S., Trentin, A. (2012). Cold-sprayed Ni-Al203 coatings for applications in power generation industry. Journal of Thermal Spray Technology, vol. 22, no. 5, p. 772-782, D0I:10.1007/s11666-013-9890-2. [24] Celotto, S., Pattison, J., Ho, J.S., Johnson, A.N., O'Neill, W. (2007). The economics of the cold spray process. Champagne, V.K. (ed.). The Cold Spray Materials Deposition Process. Fundamentals and Applications. Woodhead Publishing, Cambridge, p. 72-101, D0I:10.1533/9781845693787.1.72. [25] Kosarev, V.F., Klinkov, S.V., Papyrin, A.N. (2007). Supersonic jet/substrate interaction in the cold spray process. Champagne, V.K. (ed.). The Cold Spray Materials Deposition Process. Fundamentals and Applications. Woodhead Publishing, Cambridge, p. 178-216, D0I:10.1533/978184569 3787.2.178. [26] Alidokht, S.A., Vo, P., Yue, S., Chromik, R.R. (2017). Cold spray deposition of Ni and WC-reinforced Ni matrix composite coatings. Journal of Thermal Spray Technology, vol. 26, p. 1908-1921, D0I:10.1007/s11666-017-0636-4. [27] Lee, Y.T.R., Ahrafizedeh, H., Fisher, G., McDonald, A. (2017). Effect of type of reinforcing particles on the deposition efficiency and wear resistance of low-pressure cold-sprayed metal matrix composite coatings. Surface and Coatings Technology, vol. 324, p. 190-200, D0I:10.1016/j. surfcoat.2017.05.057. [28] Raoelison, R.N., Verdy, Ch., Liao, H. (2017). Cold gas dynamic spray additive manufacturing today: deposit possibilities, technological solutions and viable applications. Materials & Design, vol. 133, p. 266-287, D0I:10.1016/j. matdes.2017.07.067. [29] Irissou, E., Legoux, J.-G., Arsenault, B., Moreau, C. (2007). Investigation of Al-Al203 cold spray coating formation and properties. Journal of Thermal Spray Technology, vol. 16, no. 5-6, p. 661-668, D0I:10.1007/s11666-007-9086-8. [30] EN ISO 6507-1:2018. Metallic Materials - Vickers Hardness Test - Part 1: Test Method. International Organization for Standardization, Geneva. [31] Helfritch, D.J. (2007). Electromagnetic interference shielding by cold spray particle deposition. Champagne, V.K. (ed.). The Cold Spray Materials Deposition Process. Fundamentals and Applications. Woodhead Publishing, Cambridge, p. 315-326, D0I:10.1533/9781845693787.3.315. 504 Trelka, A. - Zórawski, W. - Góral, A. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 505-512 © 2020 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2020.6665 Original Scientific Paper Received for review: 2020-03-11 Received revised form: 2020-06-18 Accepted for publication: 2020-07-14 Bristle Motion, Forces, and Related Vertical Translation for a Novel Electric Toothbrush Design Frank Goldschmidtboeing1*- Uwe Pelz1 - Karen Claire-Zimmet2 - Michael Wolf2 - Ralf Goerlach2 - Peter Woias1 1 University of Freiburg, Department of Microsystems Engineering, Germany 2 Procter & Gamble Service GmbH, Germany This paper presents a combination of theoretical and experimental techniques applied to characterize the bristle motion, forces, and related vertical translation for a novel electric toothbrush design with a linear drive system. Results of the theoretical description, based on a single filament, were successfully compared with laboratory-based investigations: force measurements and high-speed video analysis, and tracking the toothbrush motion. This work describes the vertical translation induced in the toothbrush head, of up to 250 pm, when the toothbrush bristles are applied against a contact surface at brushing loads of approximately 1N to 2.5 N. Using these techniques, including Fast-Fourier transform analysis, it is shown that the vertical motion of the head is composed of the driving frequency and its harmonics. Keywords: single filament theory, electric toothbrush, bristle motion, bristle forces, inclined toothbrush filaments Highlights • Single filament based theory for tooth brushing considiring bristle driven vertical translation of the brush head. • Theoretical prediction and experimental validation of bristle peak forces at turning points of the oscillating-rotating motion. • High-speed videography visualisation of brush head motion. • Frequency analysis of the brush head motion for different drive types. 0 INTRODUCTION It may be surprising to learn that toothbrush designs have evolved over centuries. What may seem, at first thought, to be quite simple devices, are a key component of mechanical plaque removal and, therefore, oral hygiene in developed regions globally. Furthermore, far from being simple, many improvements have been made in how toothbrushes aid in plaque removal, including the introduction of electric toothbrushes in the mid-20th century. In general, electric toothbrushes include a motor and power source in the handle, along with a gear system that drives the mechanical motion of the toothbrush head and the bristles. Electric toothbrushes vary in design, including not only the type of motion delivered, but also toothbrush head size and shape, and even in the size, type, and arrangement of bristles included within the head. A wide variety of electric toothbrush models are currently found on the market, with by and large two main types of brush motion seen. One type of motion can be described as utilizing oscillating-rotating motion, in which the brush head (typically round or slightly oval) moves back and forth, at a given angle, around the central point of the brush head bristle surface. A second type of motion, often called "sonic" by toothbrush manufactures, typically involves a rectangular or oblong shaped head, which moves back and forth by rotating at a given angle relative to the long axis of the handle. The cleaning efficacy and benefits to oral health by electric toothbrushes [1] to [5], and in particular rechargeable toothbrushes with oscillating-rotating motion [6] and [7], have been extensively demonstrated via clinical measurement techniques Additional developments in electric toothbrushes include features to aid user performance, including for instance toothbrush timers, which assist the user in brushing for an appropriate duration. More recent developments include pressure sensors which may respond, including signalling to the user, when brushing is performed at a given, higher pressure threshold - thereby reducing time spent brushing with excessive force [8]. Most recently, electric toothbrushes have also been developed that incorporate position detection and interactivity with smartphone apps, helping users brush not only with appropriate duration and pressure but also ensuring they brush all areas in the mouth evenly [9] and [10]. Regarding toothbrushes with oscillating-rotating (O/R) motion, the mechanical drive systems used to move the toothbrush head have also evolved. In the 1990s, Oral-B launched a novel drive system which introduced a driven three-dimensional (3D) pulsation in addition to O/R motion, which has been demonstrated to improve cleaning efficacy via both laboratory [11] and clinical [12] testing methodologies. *Corr. Author's Address: University of Freiburg, Department of Microsystems Engineering, Laboratory for Design of Microsystems, Germany, fgoldsch@imtek.de 505 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 505-512 While cleaning efficacy is typically examined by clinical studies that are published in dentistry journals, there is almost no literature on the mechanical foundations or physics of toothbrush filament motion. This lack of publicly available, peer-reviewed research is likely due to the complicated interaction of the different filaments with each other and with the complex tooth surface. Our approach is to model single filaments and to determine their displacement and interaction force with the tooth. While the interaction between different filaments is neglected, high-speed videos show that toothbrush filaments typically move synchronously and, therefore, a basic model may provide insight, to aid in understanding the behaviour of the complex system. This approach was also used in [13], Regarding toothbrush head design, a study on the impact of bristle positioning in an O/R head, and specifically the angling of bristles, has been presented with a semi-analytical model and discussion on impact to brushing shear force and penetration depth [8], In that work, a filament with fixed length L was treated as a cantilever, and the model assumed a fixed distance between the brush and a planar tooth surface d. (see Fig. 1). Beyond the developments described above, further work has been performed to improve electric toothbrush designs, with a twofold goal: increased performance and improved brushing experience. To accomplish this goal, a completely new motor and drive train was designed, incorporating a linear motor drive system, which enables a higher brush frequency and operates with a constant oscillation angle even under load. In developing the new motor drive system, extensive technical, consumer, and clinical testing [14] to [17] were performed. Elements under investigation included overall brush motion, vibration and related experience perception, all of importance, particularly since the 3D pulsation drive introduced previously is not included within the new drive train. This work explains theoretical and experimental methods applied to describe and characterize the vertical translation, and related bristle motion and forces, for this novel electric toothbrush design. 1 METHODS The main novelty in the theoretical description is the instruction of a bristle-driven brush head motion into the semi-analytic theory developed in [13], In that publication, the distance between the brush head and tooth surface d was assumed to be constant. However, high-speed videography shows that d oscillates by typically about a few 100 |im. We have changed the boundary conditions of our model to account for this observation. The handle is still kept at a fixed position, but the brush head can be moved due to the finite stiffness of the handle. We furthermore do not label the different load cases by their distance d but by their average force that is applied by the user to the tooth. This average normal force FN0 compresses the filament such that the brush-head-to-tooth surface distance takes the value d0. During the filament movement, the normal force FN and the distance d might change. These changes can be related to each other by assuming a linear elastic behaviour of the handle. FN=FN0+K-(d-d0). (1) Here K is the scaled spring constant of the handle, i.e., the ratio of the normal force of a single filament to the brush head displacement divided by the number of filaments. This approximation implies that all filaments of a particular type act synchronously on the single brush head. This assumption is a simplification that renders the complex problem accessible for a single filament-based theory. As a typical brush head consists of filaments of different dimensions and inclinations, the different filaments act with different forces on the brush head, which is not covered in Eq. (1). Nevertheless, one type of filament typically dominates over the other types, which makes this simplification reasonable. Second, though inertial forces do not play an essential role for a single filament, the brush itself lias a considerable mass, so that inertia plays a role for the brush head dynamics. This will lead to a delayed vertical translation of the brush head compared to this simplified theory. Furthermore, the geometric relation between projected length /.,.. tip displacement wtip and brush-head-to-tooth surface distance d holds as in [8] (Fig. 1): d = cosaLx - wtjp ■ sina. (2) 506 Goldschmidtboeing, F. - Pelz, U. - Claire-Zimmet, K. - Wolf, M. - Goerlach, R. - Woias, P. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 505-512 Fig. 1. Sketch of the theoretically examined filament shapes We distinguish between the local coordinate system of the filament, where x is the coordinate along the filament and z is the coordinate perpendicular to it and the local coordinate system of the tooth surface with shear force 4 and normal force fn. Both coordinate systems are linked by the angle a (Fig. 1). The forces fx and f: are expressed by the shear force /•'-, and the average normal force fm with the help of the Eqs. (1) and (2). fx = ft sina + + (F„0 + K ■ (cos aLx - wtip ■ sina - dQ)) • cos a, (3) 4 = ft - cos a + + + K -(cosaLx -w/ip -sina-£/„))-sina, (4) Now the filament's bending line can be expressed and /•'-, by solving the Euler- w in terms of Lx Bernoulli-equation, including normal force. tip v( x,w,. ,4,4) =---^-- tan (À-Lx)-À-Lx | sin (A • x) - A ■ x + tan (A • Lx ) • (l - cos (A • x)) j, (5) , 4-sina A~ =-- ei {fno + K ■ (cosaLx - wtip ■ sina -dQ ))• cosa ei • (6) To obtain a bending line that only depends on one parameter, two equations have to be found to express the mutual dependence of /.,.. wtip and /•'-,. With these equations, all parameters (i.e., forces and the bending line) of the filament can be expressed in dependence on tip displacement wtip. The first equation states that the bending force at the tip /•'_- calculated from the bending line Eq. (5) equals that from Eq. (3). The second equation is given by the fact that the filament length does not change. Here ei is the flexural rigidity of the bristle. ei ■ w. A' ■ = ft - cosa - "p A-4-tan(A-4) (¿V0 + K ■ (cosa4 ~ • sina - dQ )) • sina, (7) L- \ l1 + ^ldxaL.Lt+i^l = 0. (8) J V dx 2 dx o As mentioned, solving Eqs. (7) and (8) for given w^ yields all information about the filament's forces and displacements. The above calculations are all done under the assumption that the average distance d0 is known. We assume that the average distance d0 is the same as the static distance at the same normal force. In static conditions, the tangential force ft should be near zero, which means that the static case is given by a displacement wtip that fulfils Eq. (9). 4 =o. (9) From this solution, the equilibrium brush-head-to-tooth surface distance d0 and the average normal force per bristle fn0 can be calculated. When modelling the motion of the filaments, it turns out that the tangential stiction force fstjc has to be considered, as it delivers the reason that the filament's tip position becomes stuck at each turning point of the rotational motion. The filament can only move along the tooth surface as long as it overcomes the maximum stiction force fstjc. a common model for the stiction force is given by Eq. (10). \f.„ \ = h-fn. (10) Here pi is the static friction coefficient. The stiction force is negative, if the filament moves in a forward direction, i.e., from left to right in Fig. 1, and consequently positive for backward movement. Solving Eq. (10) with m;i,. = ft results in the maximum normal force per bristle 4?,max- If a bristle moves from a turning point in a forward direction, the tip becomes stuck until the resulting shear force 4 overcomes the maximum stiction force in a negative direction, according to Eq. (10). Consequently, the positive maximum stiction force has to be overcome if coining from a turning point and moving in a backward direction. When the maximum stiction force is overcome, dynamic friction occurs. Its value Bristle Motion, Forces, and Related Vertical Translation for a Novel Electric Toothbrush Design 507 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 505-512 depends not only on the normal force but also on the velocity. Furthermore, a typical brush head consists of tufts with different inclinations, directions, and even different filament lengths. The full quantitative modelling of this complicated friction pattern is extremely challenging; nevertheless, a simplified model that describes the main features of the brush head is useful to understand and optimize the brushing performance. As an example, we model the new toothbrush head (Fig. 2) used with the linear drive system, which has approximately nF = 2900 bristles angled in the forward direction, nB = 600 bristles angled in the backward direction and n0 = 200 non-inclined bristles. To further simplify the model, we assume a constant length of L = 8.8 mm, a radius of r = 76 |im and an inclination angle of a = 15° for each filament, though a few filaments are slightly longer. The filaments are made of polyamide (PA 6.12), which has a Youngs' modulus of E = 2.7 GPa. A friction coefficient of |i = 0.3 was assumed. Most of the time, all the bristles will exert a total normal force /■'-, - rT1, n on the tooth. At the respective turning points, forces of /*'-,-, and Fm will occur. Therefore, we expect two force peaks with the height of Fm and Fm with a phase shift of 180°. F„ = (* FN,I={"B+"0)-FN0+"F-FN, FN,2={"F+"0)-FN0+"B-FN (11) (12) (13) 2 EXPERIMENTAL The experimental methods to examine the bristle driven head motion and the related forces are twofold: First, the motion of the brush head was visualized and tracked by high-speed video imaging; second, force measurements were conducted. To achieve repeatable conditions, all measurements were done under dry conditions, i.e., without water or toothpaste. Toothpaste and water will have an impact on the friction coefficient and the Youngs' modulus of the filaments. Therefore, we expect that the forces while cleaning will be lower than those measured in our laboratory set-up. In preparation for high-speed videography, laser-structured marking dots (o 450 |im) and scale bars (1930 |im) on tape were carefully attached to relevant moving and stationary parts of the brush heads, enabling scaled movement tracking within the video frames (Fig. 3). Fig. 3. Brush head with tracking markers in a) side- and b) front view The brush handle was secured in two pairs of custom, 3D printed clamps, ensuring the possibility for vibration only in the brush head. The clamps were attached to an aluminium stand, which was mounted to a movable stage. By using the stage movement capabilities and pressing the brush head on a flat, polished acrylic block on top of a scale (Sartorius, ED822), the force applied to the brush head could be adjusted (Fig. 4). Fig. 4. Measurement set-up for high-speed videography 508 Goldschmidtboeing, F. - Pelz, U. - Claire-Zimmet, K. - Wolf, M. - Goerlach, R. - Woias, P. direction inclined in backward direction Fig. 2. Filament pattern of the novel brush head Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 505-512 Brush head driving signal. As mentioned in the theory part, this is attributed to the inertia of the brush head. Fig. 5. a) Measurement set-up for transient force measurements; b) using a flat or c) profiled surface High-speed video of the brush head in front and side view was recorded on a Chronos 1.4 highspeed camera from Kron Technologies Inc. (Burnaby, Canada) at 4000 frames per second. All recordings have been analysed using the free-to-use software Tracker (www.physlets.org/tracker/, version 5.0.6) to extract the relative movement of all micro markers in each frame. A similar set-up was used for the force measurement (Fig. 5). The brush is pressed to a flat aluminium plate (Fig. 5b) and alternatively on an aluminium plate with a tooth-like topography (Fig. 5c). The aluminium plate is fixed on a force sensor (Kistler Model 9217A) with response threshold of 1 mN. The senor signal is amplified with a Kistler Model 5018 charge amplifier and recorded with a sample time of 50 |is. The brush head is free to move according to its elasticity as only the other end of the handle is fixed. All experiments were recorded for a total duration of 2 s. The averaged normal force Fm is determined by averaging the transient force signal. 3 RESULTS & DISCUSSION The vertical translation of the brush head and the rotation angle versus time for the average normal force between 1 N and 2.5 N are shown in Fig. 6. The periodicity of the up-and-down motion (black curve) is the same as that of the linear drive (red curve), but the force peaks are delayed compared to the Fig. 6. Vertical translation (black) and rotation angle (red) versus time Furthermore, the total stroke versus average normal force is plotted in Fig. 7. The difference between the maximum and minimum vertical translation was determined for each period and averaged for the whole signal. Theory predicts an increase of stroke with increasing normal force. The experiment, however, shows a maximum value of about 250 |im and a decreasing stroke for forces higher than 2 N. This is explained by a sidewise bulging of the filaments, which is not modelled by the theory. 1 normal f Fig. 7. Total vertical stroke vs. average normal force Fm A more direct comparison of the experiment to theory can be made from the force data. Fig. 8 shows a time series of force data for an average normal force of 2 N. The black curve shows Bristle Motion, Forces, and Related Vertical Translation for a Novel Electric Toothbrush Design 509 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 505-512 the unfiltered force data and the red marks show the automatically detected force peaks with a circular shape for the force Fm, and a diamond shape for the force Fm. The dashed lines show the averaged values. Fig. 8. Transient force data for the novel brush head at an average normal force of 2 N These data show much more details compared to the displacement data, as short spikes of force do not directly transform to displacement due to the inertia of the brush head. This transient curve shows two force maxima per period which are separated by exactly half a period and show two different values. These values are interpreted as force maxima /*'-,-, (Eq. (12)) for the turning point to forward motion and Fm (Eq. (13)) for the turning point to backward motion. The averaged experimental values for these force peaks and the theoretical predictions are compared in Fig. 9, where the red marks indicate the force maximum /*'-,-, while the black marks indicate the force maximum Fm. Each point was measured three times to assess the repeatability. The lines indicated the theoretical values (Eqs. (12) and (13)). The results scatter significantly, as it was relatively difficult to adjust the same conditions in every experiment. Nevertheless, both curves match well, especially for high forces. The experimental maximum forces, however, show relatively high deviations at the lowest normal force of 1 N. The reason for this deviation is probably a vibration of the sensor plate itself that was observed at low normal force values. Fig. 9. Comparison of experimental and theoretical magnitudes of the force peaks at the turning points To further examine the vibrations of the brush head under realistic conditions, force measurements were conducted on the aluminium plate with the toothlike topography. Fig. 10 shows a fast Fourier transform (FFT) of the force data for an average force of 2 N with a linear O/R dive. Only the drive frequency of about 147 Hz and its harmonics are present. Fig. 11 shows the FFT of a comparable experiment with a brush with additional 3D pulsation (see [8] for details on the brush design). In this case, both drives, the O/R drive with a frequency of 93 Hz and the 3D pulsation drive with a frequency of 417 Hz interact nonlinearly, which leads to complicated frequency pattern with many different frequency components in the sensible frequency range of a few hundred Hz. Not only those two mismatched frequencies but also the sums and difference of those and their multiples are present. 200 300 Frequency [Hz] Fig. 10. FFT analysis of the frequency spectrum of the normal force with the novel brush design and the linear drive 510 Goldschmidtboeing, F. - Pelz, U. - Claire-Zimmet, K. - Wolf, M. - Goerlach, R. - Woias, P. Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 505-512 100 200 300 400 Frequency [Hz] Fig. 11. FFT analysis of the frequency spectrum of the normal force with the O/R drive with 3D pulsation 4 CONCLUSIONS A single filament model was utilized to describe the complex bristle motion of an electric toothbrush, including the effect of bristle driven vertical translation of the brush head. The main predictions of the simple model, such as the value of the force peaks at the turning points, the periodicity and the order of magnitude of the vertical translation were confirmed by the experiment. Thus, the simple model helps the designers of toothbrushes to optimize their design. Other details, such as the reduction of vertical translation amplitude at high average normal forces or the exact transient behaviour of the brush head were experimentally examined while the model was not capable of predicting these effects quantitively, because effects like sidewise bulging and the inertia of the brush were not included into the single-filament-based model. Therefore, a mixed theoretical and experimental approach is necessary to design and understand the complex movement patterns of driven brush heads. The experiments have shown that an electric brush with linear O/R-drive produces force peaks at the turning points that can be adjusted by the ratio of filaments inclined into different directions. A vertical translation of up to 250 ^m amplitude is generated by the novel design. This amplitude is influenced by properties, such as the stiffness, orientation, and quantity of filaments, as well as the handle stiffness, and can therefore be modulated. The current value has been perceived as pleasing by test consumers. The force signal of the brush with linear drive consists only of the drive frequency and its harmonics. Though toothbrushes are in daily use and have a history of centuries, there are still open questions that may be addressed by theoretical and experimental methods. 5 ACKNOWLEDGEMENTS This work was funded by Procter & Gamble Service GmbH, Germany. 6 REFERENCES [1] Rosema, N., Slot, D.E., van Palenstein Helderman, W.H., Wiggelinkhuizen, L., Van der Weijden, G.A. (2016). The efficacy of powered toothbrushes follong a brushing exercise: a systematic review. International Journal of Dental Hygiene, vol 14, no. 1, p. 29-41, D0l:10.1111/idh.12115. [2] Heanue, M., Deacon, S.A., Deery, C., Robinson, P.G., Walmsley, A.D., Worthington, H.V., Shaw, W.C., Shaw, B.C. (2003). Manual versus powered toothbrushing for oral health. Cochrane Database of Systematic Reviews, no. 1, art. no. CD002281, D0I:10.1002/14651858.cd002281. [3] Yaacob, M., Worthington, H.V., Deacon, S.A., Deery, C., Walmsley, A.D., Robinson, P.G., Glenny, A.M. (2014). Powered versus manual toothbrush for oral health. Cochrane Database of Systematic Reviews, no. 6, art. no. CD002281, D0I:10.1002/14651858.CD002281.pub3. [4] Van der Weijden, F.A., Slot, D.E. (2015). Efficacy of homecare regimens for mechanical plaque removal in managing gingivitis a meta review. Journal of Clinical Periodontology, vol. 42, no. S16, p. S77-S91, D0I:10.1111/jcpe.12359. [5] Pitchika, V., Pink, C., Völzke, H., Welk, A., Kocher, T., Holtfreter, B. (2019). Long-term impact of powered toothbrush on oral health: 11-year cohort study. Journal of Clinical Periodontology, vol. 46, no. 7, p. 713-722, D0I:10.1111/jcpe.13126. [6] Grender, J., Adam, R., Zou, Y. (2020). The effects of oscillating-rotating electric toothbrushes on plaque and gingival health: A meta-analysis. American Journal of Dentistry, vol 33, no. 1, p. 3-11. [7] Clark-Perry, D, Levin, L. (2020) Systemic review and metaanalysis of randomized controlled studies comparing oscillating rotating and and other powered toothbrushes. The Journal of the American Dental Association, vol. 151, no. 4, p. 265-275, D0I:10.1016/j.adaj.2019.12.012. [8] Janusz, K., Nelson, B., Bartizek, R.D., Walters, P.A., Biesbrock, A. (2008). Impact of a novel power toothbrush with SmartGuide technology on brushing pressure and thoroughness. Journal of Contemporary Dental Practice, vol. 9, no. 7, p. 1-8. [9] Erbe, C., Klees, V., Ferrari-Peron, P., Cahuana-Vasquez, R.A., Timm, H., Grender, J., Cunningham, P., Adam, R., Farrell, S., Wehrbein, H. (2018). A comparative assessment of plaque removal and toothbrushing compliance between a manual and an interactive power toothbrush among adolescents: a single-center, single-blind randomized controlled trial. BMC Oral Health, vol. 18, no. 1, p. 130, D0I:10.1186/s12903-018-0588-1. [10] Erbe, C., Klees, V., Braunbeck, F., Ferrari-Peron, P., Ccahuana-Vasquez, R.A., Timm, H., Grender, J., Cunningham, P., Adam, R., Wehrbein, H. (2019). Comparative assessment of plaque Bristle Motion, Forces, and Related Vertical Translation for a Novel Electric Toothbrush Design 511 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 505-512 removal and motivation between a manual toothbrush and an interactive power toothbrush in adolescents with fixed orthodontic appliances: A single-center, examiner-blind randomized controlled trial. American Journal of Orthodontics and Dentofacial Orthopedics, vol. 155, no. 4, p. 462-472, DOI:10.1016/j.ajodo.2018.12.013. [11] Driesen, G.M., Warren, P.R., Hilfinger, P. (1998). Cleaning efficacy of a new electric toothbrush. American Journal of Dentistry, vol. 11, special issue, p. S7-S11. [12] Ernst, C.P., Nauth, C., Willershausen, B., Warren, P.R. (1998). Clinical plaque removing efficacy of a new power toothbrush. American Journal of Dentistry, vol. 11, special issue, p. S13-S16. [13] Goldschmidtboeing, F., Doll, A., Neiss, S., Doll A., Stoerkel, U., Woias, P. (2014). Comparison of vertical and inclined toothbrush filaments: Impact on shear force and penetration depth. Strojniški vestnik - Journal of Mechanical Engineering, vol. 60, no. 7-8, p. 449-461, DOI:10.5545/sv-jme.2013.1513. [14] Adam, R. (2020). Introducing the Oral-B iO electric toothbrush: next generation oscillating-rotating technology. International Dental Journal, vol. 70, no. S1, p. S1-S6, DOI:10.1111/ idj.12569. [15] Grender, J., Goyal, C.R., Qaqish, J., Adam, R. (2020). An 8-week randomized controlled trial comparing the effect of a novel oscillating-rotating toothbrush versus a manual toothbrush on plaque and gingivitis. International Dental Journal, vol. 70, No. S1, p. S7-S15, DOI:10.1111/idj.12571. [16] Adam, R, Goyal, C.R., Qaqish, J., Grender, J. (2020). Evaluation of an oscillating-rotating toothbrush with microvibrations versus a sonic toothbrush for the reduction of plaque and gingivitis: results from a randomized controlled trial. International Dental Journal, vol. 70, no. S1, p. S16-S21, DOI:10.1111/idj.12569. [17] Adam, R., Erbe, J., Grender, J. (2020). Randomized controlled trial assessing plaque removal efficacy of a novel oscillating-rotating electric toothbrush with micro-vibrations. International Dental Journal, vol. 70, no S1, p. S22-S27, DOI:10.1111/ idj.12568. 512 Goldschmidtboeing, F. - Pelz, U. - Claire-Zimmet, K. - Wolf, M. - Goerlach, R. - Woias, P. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 513-522 © 2020 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2020.6667 Original Scientific Paper Received for review: 2020-03-13 Received revised form: 2020-06-07 Accepted for publication: 2020-07-27 Research on the Lubrication Characteristics of Harmonic Gear Transmission Meshing Areas Xiangyang Xu1 * - Xupeng Fan1 - Peitang Wei2 - Baojun Yang3 1 Chongqing Jiaotong University, School of Mechatronics and Vehicle Engineering, China 2 Chongqing University, College of Mechanical Engineering, China 3 Chongqing Huashu Robot Co., China To analyse the lubrication characteristics of harmonic gears and lay the foundation for the study of its gear tooth failure performance and dynamic characteristics, based on the tooth contact geometry of harmonic gear, the integrated curvature radius, tooth load, and entrainment velocity at the meshing point of the gear teeth in the harmonic gear transmission are analysed. A finite-length line contact elastohydrodynamic lubrication (EHL) model for harmonic gears is established. The numerical calculation method is used to solve the oil film thickness and pressure distribution in the lubricating area, and the effects of rotational speed and temperature on the contacting load ratio and film thickness ratio of the meshing area are studied, as well as the change of oil film stiffness under different working conditions. The results show that along the meshing direction, the pressure is small at the end and reaches a peak at the centre, and the film thickness is the largest in the entrance area and is evenly distributed in the centre contact area. As the speed increases, the gear tooth contact load ratio decreases, the oil film thickness ratio increases, the stiffness of the oil film decreases significantly and the lubrication effect is improved; but the temperature has the opposite effect. Proper increase of rotation speed and decrease of oil temperature can effectively improve the lubrication characteristics of the system. Keywords: harmonic gear transmission; elastohydrodynamic lubrication; meshing area; lubrication characteristics Highlights • The integrated curvature radius, the entrainment velocity and the load at the meshing point of harmonic gear are analysed and extracted for the harmonic gear lubrication. • A finite line contact elastohydrodynamic lubrication model for harmonic gears is established. • A complex iterative method is used to solve the lubrication equation to obtain high precision. • The lubrication characteristics and the change of oil film stiffness under different working conditions are explored. 0 INTRODUCTION Harmonic gear transmission is a new type of transmission mode, which is different from traditional gear transmission. It mainly transfers power through the deformation of flexible gears, with small volume, high precision, large transmission ratio and other excellent performance features. Therefore, it is widely used in the field of high-precision machinery, such as industrial robots, aerospace equipment, ect. [1] and [2], However, factors such as nonlinear stiffness, hysteresis and friction nonlinear in the harmonic gear drive restrict its improvement of the accuracy and efficiency. The interaction between lubrication oil and gears affects the power loss of gears and transmission accuracy [3] and [4], Thus, it is necessary to study the lubrication performance of the gear system to improve its efficiency and accuracy. Therefore, the analysis of the geometry and motion parameters of harmonic gear teeth and the research of harmonic gear lubrication control can promote the analysis of the friction and failure mechanism of the harmonic gear transmission and improve the nonlinear behaviour of the system. Concli and Goria [5] and [6] used an improved CFD simulation method to predict the power loss caused by the interaction between gears and lubricant, simulated the wind resistance, agitation and cavitation with a comprehensive, time-saving numerical simulation method, and compared the experimental data with the numerical simulation results in terms of power loss and lubricant distribution with good agreement. In the study of elastohydrodynamic lubrication, various numerical solutions were mainly aimed at infinite line contact and point contact problems; however, in engineering practice, most lubrication problems are finite line contact problems. In the static contact elasticity theory [7], the contact pressure distribution of the finite line contact pair with straight generatrix profile is different from that of the classic Hertz infinite line contact theory, and there is an "edge effect" at the end part. Furthermore, the contact width of the gear is limited, so the infinite line contact model cannot adequately simulate the actual situation of the end part. Researchers have carried out some studies on the problem of finite-line wire contact *Corr. Author's Address: Chongqing Jiaotong University, ChongQIng, 400074, China, xxlangyang@hotmail.com 513 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 513-522 elastohydrodynamic lubrication. Wymer and Cameron [8] studied the elastohydrodynamic lubrication of tapered roller utilizing an optical interference method, and obvious edge effect was observed at the end. Park and Kim [9] established the elastohydrodynamic lubrication model of finite length line contact and found that in the middle of the roller, the pressure distribution is almost the same as that of the infinite line contact. However, near the end part, the results of EHL are quite different from that of the infinite solution, with which the pressure reaches its peak. Zhu et al. [10] established a finite length EHL model for helical gears and found that the second pressure peak occurred at the end part. Wei et al. [11] established a mixed lubrication model for cycloidal pinwheels and discussed the effects of contact load and radius of curvature on lubrication conditions. Liu et al. [12] developed a line contact model of thermoelastic hydrodynamic lubrication under the condition of an ultrathin film and studied the influence of speed and dynamic load on the lubrication performance of a spur gear pair. Ouyang et al. [13] established an integrated friction dynamic lubrication model of spur gear pairs considering dynamic load and elastohydrodynamic lubrication and analysed the relationship between the dynamic model and the elastohydrodynamic lubrication model through an iterative algorithm. Xiao et al. [14] discussed the distribution of central pressure and film thickness in the lubrication area before and after considering the thermal effect and analysed the influence of heat on the film damping, and the influence of contact force, speed and number of teeth on the film damping. Sikorski and Pawlowski [15] added graphite and molybdenum disulphide into the lubricant, and through the experimental verification of the two kinds of ball bearings, measured the power of the internal friction of the ball bearings under the action of both at low temperature with the change of temperature, and found that the two additives can effectively improve the lubrication performance. Scaraggi et al. [16] analysed the lubrication of pin-pulley interface in CVT of gear chain industrial chain and found that the lubrication mode at the interface was determined by a series of hydrodynamic and hybrid lubrication stages. Ciulli et al. [17] presented a model for predicting scuffing considering various lubrication states during the operation of the machine, and the effectiveness of the model is verified through experiments. Li et al. [18] used the elastohydrodynamic lubrication (EHL) method to calculate oil film thickness and pressure distribution in predicting roller bearings' skidding for obtaining more accurate oil film resistance. Li et al. [19] carried out a life test on the space harmonic reducer, finding that the flexspline and bearings would wear under the condition of separate grease lubrication, and analysed the effects of different factors such as load and temperature on lubrication. Li et al. [20] calculated the stiffness of the elastohydrodynamic oil film in the point contact area, established an oil film stiffness calculation model, and verified the accuracy of the model through experiments. However, at present, few scholars have quantitatively analysed the influence of working conditions on the lubrication characteristics and oil film stiffness of gear meshing areas based on tooth contact geometry. Therefore, the entrainment velocity, the integrated curvature radius and the load distribution in the meshing area of harmonic gears are analysed in the present paper. Then, the elastohydrodynamic lubrication model of harmonic gear transmission based on the finite line contact is established. And the numerical method is used to solve the lubrication area. The lubrication characteristics and the oil film stiffness in meshing area are quantitatively analysed with the change of rotational speed and oil temperature. 1 METHODS The elastic deformation equation, lubricant viscosity equation, and density equation are functions of pressure, which are solved by the compound iterative method after being combined with the Reynolds equation. The specific steps are as follows: First, the basic parameters, initial pressure and initial geometric film thickness are input, and the initial elastic deformation is obtained by integral solution, and the initial film thickness is obtained by substituting the obtained results into the film thickness formula. Second, the initial film thickness is substituted into Reynolds equation, and all the nodal pressures can be obtained. Then, the elastic deformation and film thickness are recalculated by the obtained pressure, and the iterative process is repeated using the above calculation results as the input of the next iteration until the convergence solution with accuracy better than 10-4 is obtained. 2 EXPERIMENTAL 2.1 Principle of Harmonic Gear Transmission As shown in Fig. 1, the harmonic gear transmission system is mainly composed of a wave generator, a flexspline and a circular spline. Before assembly, 514 Xu, X - Fan, X. - Wei, P. - Yang, B. Strojnlskl vestnik - Journal of Mechanical Engineering 66(2020)9, 513-522 the original profile of the flexspline is circular, and the circular pitch of the flexspline is the same as that of the circular spline, but the number of teeth of the flexspline is slightly less than that of the circular spline, while the maximum diameter of the wave generator is slightly larger than that of the inner circle of the flexspline. When the wave generator is installed in the flexspline, the deformation of the flexspline is forced, and the initial deformation force is generated. Thus, the power is transmitted through the meshing in and out of deformation of the flexspline as a flexible component. The lubrication characteristics of flexspline and circular spline in the conjugate meshing area the main research target in the next step. Fig. 1. Schematic diagram of harmonic gear drive The basic parameters of harmonic gear transmission are shown in Table 1. Table 1. fias/c parameters of harmonic gear transmission Parameters Value Number of flexspline 200 Number of circular spline z 202 Modulus [mm] 0.5 Tooth width [mm] 12 Transmission ratio 100 Pressure angle 20° Rated speed of wave generator [rpm] 3000 Output torque [N m] 400 Circular spline material density [kg/m3 7.84x103 Elastic modulus of circular spline jT2 [Pa] 2.06x1011 Poisson's ratio of circular spline material v2 0.3 Material density [kg/m3 7.83x103 Elastic modulus of flexspline _£\ [Pa] 2.01x1011 Poisson's ratio of flexspline material Vj 0.3 Surface roughness of flexspline [)um] 0.34 Surface roughness of circular spline [/jm] 0.35 2.2 Establishment of Elastohydrodynamic Lubrication Model 2.2.1 Geometric and Kinematic Analysis For the convenience of discussion of lubrication characteristics between flexspline and circular spline in the conjugate meshing area, as shown in Fig. 2, in the harmonic gear transmission, the meshing of the flexible and the circular spline be approximated as two finite-length cylindrical rollers contacting each other, where x is the direction of motion and Rki, Rk2 are the curvature radii of meshing point of the flexspline and the circular spline, respectively. Fig. 2 Finite-length line contact model Fig. 3 shows thctransmission relationship during meshing. R and G are the tooth curves of the flexspline and the circular spline, k is the meshing point, Obl and Ob2 are the corresponding base circle centres of the flexspline and the circular spline, respectively, rM and rb2 are the base circle radii of the meshing point k on the flexspline and the circular spline, and T is the common normal line of the two base circles, and the meshing point k is on the common normal line T. Research on the Lubrication Characteristics of Harmonic Gear Transmission Meshing Areas 515 neutral of flexspline 0{Oh2) Fig. 3. Schematic diagram of the transmission relationship of the meshing point Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 513-522 According to the above transmission relationship diagram, it can be concluded that when the gear teeth of the harmonic reducer mesh, the curvature radii of the meshing points of the flexspline and the circular spline are shown as follows: R, R, (1) where Rt is the distance from the centre of the flexspline base circle to the meshing point, and pk is the polar radius of the meshing point k. Then the comprehensive curvature radius R at the meshing point required for the lubrication calculation is obtained by the following Eq. (2): R = Rn -RH2 Rk2 Rkl (2) Fig. 4 is a schematic diagram of the velocity at the meshing point k during the movement of the teeth of the flexspline, where m represents the angular velocity; C represents the intersection of the line between the coordinate origin and the meshing point and the neutral line of the flexspline. Vek, Vrk, Vtk are the velocity of the involved motion, the radial and circumferential elastic deformation, respectively. Then, the required meshing point entraimnent velocity h in the lubrication equation is the component of the resultant velocity along the tangential direction of the tooth profile. In the meshing direction, the variation of comprehensive curvature radius R and the entraimnent speed h at each meshing position are shown in Fig. 5. Meshing jii>int position Fig. 5. The variation of comprehensive curvature radiusR and entrainment speed u Different from ordinary gear load distribution, the load distribution between the teeth of harmonic gear drive is complex in the theoretical calculation, especially when considering the effect of the lubricating oil film. Therefore, on the basis of summarizing experimental data results for when there is no lubrication, it is found that the distribution law is similar to the cosine distribution shown in the Fig. 6. Fig. 4. Schematic diagram of meshing point speed Fig. 6. Simplified curve of load distribution Given that the angular velocity of the wave generator is o)H and the angular velocity of the flexspline is '>>,. The relative motion velocity of the k-point is the combined velocity of the implication velocity (Vek) and elastic deformation velocity (Vrk, I',/, ). which can be shown by the Eq. (3): vK=vek+vrk+v,k (3) Therefore, the tooth load F can be expressed by: n{(x, v). According Fig. 7, the film thickness can be obtained as: h = h0+ — + n(x,y) + 8(x, v), IK (9) where the comprehensive elastic deformation ¿>(x, v) can be expressed as: ""'■»-M P(sj) xE JJn ^(x-s)2 +{y-t) =dsdi. (10) where e is the comprehensive elastic modulus and Q is the calculation domain. The equivalent elastic modulus is calculated by the following Eq. (11): (ID 2 e i-vr , i r. en F 2 where /•.',. Vj, v2 are the elastic moduli and Poisson's ratios of the flexspline and the circular spline materials, respectively. The viscosity and density equations are expressed as: '? = '? oexP (lnij0+9.67)(l + 5.1xlO"9p); -9..W 71-138 -1 ■(12) P = Po \+ 0-6X10-V 1-0.7x10-3(7;-r„), (13) 1 + 1.7x10/? where 7, is the lubricant temperature, t0 is the ambient temperature, ;/0 is environmental viscosity, p0 is the enviromnental density, z =-—^---, 5. lx 10 (lni70 + 9.67) a is the viscosity pressure s = -138) hi770 +9.67 " coefficient and /; is the viscosity temperature coefficient, respectively. The basic initial parameters of the lubricant are shown in Table 2. Table 2. Lubricant parameters Parameters Value Ambient temperature 7}, [K] 313 Ambient viscosity ;/n [Pa-s] 0.04 Ambient densityp(l [kg/m3] 870 Viscosity pressure coefficient a [m2/N] 2.2x10-8 Viscosity temperature coefficient p [K-1] 0.0476 The equilibrium equation of the contact area pressure and external load is: Fig. 7. Schematic diagram of film thickness Research on the Lubrication Characteristics of Harmonic Gear Transmission Meshing Areas 517 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 513-522 F = jj p (x, y )dxdy. (14) 2.2.3 Oil Film Stiffness To calculate the oil film stiffness, this paper uses the average film thickness method, as shown in Eq. (15): k -AL k ~ Ah ' (15) where Af is the load increase of the oil film, and Ah is the average oil film thickness increase from the entrance area to the exit area. Ah=I hi ( f ) ^(«2 -ni) I hi (f + Af ) («2 - ni) ' (16) where nj and n2 are the nodes where the pressure in the inlet area is generated and the pressure in the outlet area disappears. The number of nodes in x and y directions is 129 x 129. 3 RESULTS According to the elastohydrodynamic lubrication control equation, the pressure distribution and oil film thickness in the meshing area of harmonic gear are obtained by numerical calculation method, and the influence of rotation speed and oil temperature conditions on the lubrication performance parameters, such as contact load ratio, film thickness ratio and oil film stiffness of harmonic gear, is analysed. 3.1 Pressure and Oil Film Distribution First, the numerical method is used to calculate the numerical solution of the elastohydrodynamic lubrication of the harmonic gear transmission. The oil film thickness and pressure distribution between the contact interfaces are shown in Fig. 8. It can be seen that the pressure at both ends in the axial >>-direction is significantly higher than that in the middle, resulting in an end leakage effect. While in the x-direction, the pressure is greater as it approaches a) b) Fig. 8. Numerical solution of pressure and film thickness in the meshing area; a) pressure distribution, and b) oil film thickness distribution 518 Xu, X - Fan, X. - Wei, P. - Yang, B. Strojnlskl vestnik - Journal of Mechanical Engineering 66(2020)9, 513-522 the centre of the meshing point and the pressure peaks at x = 0, indicating that the pressure along the meshing line is the largest at the centre of meshing. The oil film thickness reaches its peak in the entrance area in the x-dircction. and the distribution is flat in the contact area. The oil film at both ends in the v-direction is smaller than the middle part, producing an end effect, which is because the main difference between the finite line contact and the infinite line contact is the end leakage effect, which makes the pressure at both ends higher than the middle part, while the film thickness trend is the opposite. The relevant research results have been introduced in another paper of the author [10], and similar phenomena have been found in other scientific literature [8] and [9]. 3.2 Influence of working condition on lubrication characteristics of contact area To study the actual meshing situation of the gear tooth meshing area, considering the surface roughness, two points a (conjugate tooth root meshing point) and b (conjugate the meshing point of the tooth top circle) are selected, and the analysis of contact load ratio w (roughness peak bearing load as a percentage of total load) and oil film thickness ratio a (ratio of average oil film thickness to surface roughness) are performed under different working conditions. The analysis results shown in Figs. 9 and 10 were obtained. Fig. 9 shows the change curve of the contact load ratio and film thickness ratio of the points a and b with the rotation speed. Fig. 9. Effect of rotation speed on lubrication As can be seen from Fig. 9, when the rotation speed is less than 100 r/min, the load ratio w at each point is larger (more than 80 %). When the speeds of points a and b are 300 r/min and 2500 r/min. respectively, the film thickness ratio a is close to 1, and at this time, if the rotation speed continues to decrease, a will be less than 1, and the lubrication effect will be inadequate. When the rotation speed increases, the contact load ratio w decrease, the film thickness ratio a becomes larger, and the lubrication condition in the meshing area is improved. However, when the rotation speed is too large, the effect will gradually become smaller. At low speed, the film thickness is relatively small, and the rough peak contact between the contact surfaces increases. When the speed is too low, the film thickness ratio is less than 1, and the lubrication will turn to boundary lubrication, which will accelerate gear wear. Therefore, the harmonic reducer should appropriately increase the rotation speed to prevent poor lubrication. To study the effect of temperature on the lubrication characteristics, without changing the speed and other factors, by changing the lubricant temperature, the analysis results are shown in Fig. 10. Fig. 10. Effect of oil temperature on lubrication It can be seen from Fig. 8 that as the temperature increases, the contact load ratio gradually increases, and the film thickness ratio gradually decreases, which shows that the increase in temperature reduces the thickness of the lubricant film and increases the roughness peak load. When the temperature at point a is greater than 60 °C and the temperature at point b is greater than 30 °C, the film thickness ratio will be less than 1, the oil film thickness will decrease, the lubrication will change from mixed lubrication to boundary lubrication, the rough peak contact between the contact surfaces will increase, and the lubrication effect will deteriorate. Therefore, avoiding long-term high-temperature operation is conducive to improving the lubrication of harmonic gears and increasing the gears' lifespan. Research on the Lubrication Characteristics of Harmonic Gear Transmission Meshing Areas 519 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 513-522 3.3 Effect of Working Conditions on the Change of Oil Film Stiffness To further analyse the influence of working conditions on the oil film stiffness, the changes in oil film stiffness in the meshing area of the gear teeth along the direction of the meshing line under different rotation speed and oil temperature were calculated, as shown in Fig. 11. a) b) Fig. 11. Effect of working conditions on oil film stiffness; a) effect of speed on oil film stiffness, and b) effect of temperature on oil film stiffness It can be seen that the single-tooth meshing area of the oil film in the entrance area is relatively thin due to the large film thickness. However, the stiffness of the multi-tooth meshing region is greater due to the larger pressure peak, and sudden variation occurs at the alternate points of single-tooth and multi-teeth. As the speed increases, the stiffness of the oil film decreases significantly, and the sudden variation in stiffness gradually decreases. While the effect of temperature on stiffness is just the opposite: as the temperature increases, the oil film stiffness increases significantly, and the sudden variation in stiffness also becomes apparent, which shows that properly increasing the rotation speed and lowering the oil temperature can effectively reduce the oil film stiffness, thereby improving the overall stiffness of the gear teeth and reducing the system vibration. 4 CONCLUSION For harmonic gear transmission, the parameters such as the integrated curvature radius, the entrainment velocity and the load at the meshing point of harmonic gear are analysed and extracted, and the mathematical analysis model of EHL of harmonic gear transmission based on the finite-length line contact is established. Using the numerical calculation method to solve the model, the pressure and film thickness distribution in the lubrication area are obtained. The characteristic lubrication parameters, such as the contact load ratio and film thickness ratio of the meshing points of the harmonic gear transmission under different rotational speeds and oil temperature conditions, were analysed, as well as the changes of the oil film stiffness. The results show that the pressure along the meshing direction reaches a peak at the centre meshing point, the film thickness is larger in the entrance area, and it is evenly distributed in the centre meshing area. Both the rotational speed and oil temperature have significant effect on the lubrication characteristics. A proper increase of the rotational speed can increase the film thickness ratio of meshing points and reduce the contact load ratio, which is beneficial to improving the lubrication characteristics of the meshing area and reducing the oil film stiffness. The influence of oil temperature on it is just the opposite. Increasing the temperature will reduce the film thickness ratio of the meshing point and reduce the contact load ratio, which will make the lubrication performance worse. Therefore, appropriately increasing the rotational speed and reducing the lubricant temperature can effectively improve the lubrication characteristics of the harmonic gear transmission system. The research results can provide a reference for the subsequent tooth wear failure mechanism of harmonic gear transmission and the development of nonlinear dynamic analysis. 5 ACKNOWLEDGEMENTS This research is supported by the Foundation of the National Key Research and Development Plan of China (No. 2018YFB2001300), the Foundation of 520 Xu, X - Fan, X. - Wei, P. - Yang, B. Strojnlskl vestnik - Journal of Mechanical Engineering 66(2020)9, 513-522 National Natural Science of China (No. 51975078), the Natural Science Foundation Project of Chongqing (No. cstc2018jcyjAX0087), the Key Research and Development Project of Chongqing (cstc2017rgzn-zdyfX0038), and the Technology Innovation and Application Development Project of Chongqing (No.cstc2018jszx-cyzdXO 159 and No.cstc2019jscx-mbdxX0049). 6 NOMENCLATURES c the intersection of the line between the coordinate origin and the meshing point and the neutral line of the flexspline e equivalent modulus of elasticity, [Pa] /•.',. e2 the elastic moduli of the flexspline and the circular spline materials, respectively, [Pa] f load between the flexspline and the circular spline, [N] /•',, the maximum load, [N] A/ the load increase of the oil film, [Pa] h the oil film thickness, [m] h„_ the central film thickness, [m] Ah the average oil film thickness increase from the entrance area to the exit area, [m] k the meshing point A0ll the oil film stiffness, [N/m] N the parameter of the number of meshing teeth «,. ;?2 the node numbers where the pressure in the inlet area is generated and the pressure in the outlet area disappears ()h]. Ob2 the corresponding base circle centres of the flexspline and the circular spline, respectively p the oil film pressure, [Pa] rbl, rb2 the base circle radii of the meshing point k on the flexspline and the circular spline, respectively, [Pa] R the comprehensive curvature radius at the meshing point, [m] R\ the distance from the centre of the flexspline base circle to the meshing point, [m] /<7,i- ri,2 the curvature radii of meshing point of the flexspline and the circular spline, respectively, [m] r, g the tooth curves of the flexspline and the circular spline, respectively t the common normal line of the two base circles t0 the ambient temperature, [K] 71 the lubricant temperature, [K] tm Torque, [N m] ii the entraimnent velocity, [m/s] Ithe relative motion velocity of the A-point, [m/s] Vek the velocity of the involved motion, [m/s] I,'./, the velocity of the radial elastic deformation, [m/s] I',/, the velocity of the circumferential elastic deformation [m/s] w the contact load ratio vj, v2 the Poisson's ratios of the flexspline and the circular spline materials, respectively o) the angular velocity, [rad/s] o)H the angular velocity of the wave generator, [rad/s] o)R the angular velocity of the flexspline, [rad/s] x the direction of motion v the direction perpendicular to x z the number of the circular spline teeth a the viscosity pressure coefficient /; the viscosity temperature coefficient k oil film thickness ratio p the lubricant density, [kg/m3] p0 the environmental density, [kg/m3] pk the polar radius of the meshing point k, [m] // the lubricant viscosity, [Ns/m2] //o enviromnental viscosity, [Ns/m2] u (x.y) the surface roughness, | |im | ¿>(x, v) the comprehensive elastic deformation £2 the calculation domain ([> the position of meshing point cpm the position of the maximum load 9>0 the size of the actual meshing area 7 REFERENCES [1] Tao, M„ Chen, Y„ Chen, D„ Wu, J„ Chen, F„ Li, B„ Mei, J. (2018). Research present status and outlook of harmonic reducer testing technology. Journal of Mechanical Transmission, vol. 42, no. 7, p. 175-180, D0l:10.16578/j. issn.1004.2539.2018.07.034 (in Chinese) [2] Hu, Q„ Liu, Z„ Cai, L„ Yang, C„ Zhang, T„ Wang, G. (2019). Research on prediction method of transmission accuracy of harmonic drive. International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, vol. 10, p. V010T11A006, D0l:10.1115/ DETC2019-97214 [3] Preumont P., Szewczyk, R. (2018). Key factors influencing the accuracy of harmonic gears for space applications. Conference on Automation, p. 483-489, 001:10.1007/978-3-319-77179-3_45 [4] Concli, F. (2016). Windage, churning and pocketing power losses of gears: differentmodeling approaches for different goals. Forschung im Ingenieurwesen, vol. 80, p. 85-99, DOI:10.1007/sl0010-016-0206-9 Research on the Lubrication Characteristics of Harmonic Gear Transmission Meshing Areas 521 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 513-522 [5] Concli, F., Gorla, C. (2017). CFD simulation of power losses and lubricant flows in gearboxes. American Gear Manufacturers Association Fall Technical Meeting, p. 2-14. [6] Concli, F, Gorla, C. (2016). Numerical modeling of the power losses in geared transmissions: Windage, churning and cavitation simulations with a new integrated approach that drastically reduces the computational effort. Tribology International, vol. 103. p. 58-68, D0l:10.1016/j. triboint.2016.06.046. [7] Jonson, K.L. (1985). Contact Mechanics. London: Cambridge University Press, p. 129-152. [8] Wymer, D.G., Cameron, A. (1974). Elastohydrodynamic lubrication of a line contact. Proceedings of the Institution of Mechanical Engineers, vol. 188, no. 1, p. 221-238, D0I:10.1243/PIME_PR0C_1974_188_024_02. [9] Park, T.J., Kim, K.W. (1998). Elastohydrodynamic lubrication of a finite line contact. Wear, vol. 223, no. 1-2, p. 102-109, D0I:10.1016/S0043-1648(98)00317-2. [10] Zhu, C.C., Liu, M.Y., Liu, H.J., Xu, X.Y., Liu, L.B. (2013). A thermal finite line contact EHL model of a helical gear pair. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 227, no. 4, p. 299-309. D0I:10.1177/1350650112462324. [11] Wei, B., Wang, J., Zhou, G., Yang, R., Zhou, H., He, T. (2016). Mixed lubrication analysis of modified cycloidal gear used in the RV reducer. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 230, no. 2, p. 121-134, D0I:10.1177/1350650115593301. [12] Liu, H., Mao, K., Zhu, C., Chen, S., Xu, X., Liu, M. (2013). Spur gear lubrication analysis with dynamic loads. Tribology Transactions, vol. 56, no. 1, p. 41-48, D0I:10.1080/1040200 4.2012.725805. [13] Ouyang, T., Chen, N., Huang, J., Huang, H. (2016). Analysis of lubricating performance for spur gear pairs applying tribo- dynamic model. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 230, no. 10, p. 1244-1257, D0l:10.1177/1350650116631713. [14] Xiao, Z., Li, Z., Shi, X., Zhou, C. (2018). Oil film damping analysis in non-Newtonian transient thermal elastohydrodynamic lubrication for gear transmission. Journal of Applied Mechanics, vol. 85, no. 3, D0I:10.1115/1.4038697. [15] Sikorski, J., Pawlowski, W. (2020). Internal Friction of Ball Bearings at Very Low Temperatures. Strojniški vestnik -Journal of Mechanical Engineering, vol. 66, no. 4, p. 235-242, D0I:10.5545/sv-jme.2019.6398. [16] Scaraggi, M., De Novellis, L., Carbone, G. (2010). EHL-Squeeze in Highly Loaded Contacts: The Case of Chain CVT Transmissions. Strojniški vestnik - Journal of Mechanical Engineering, vol. 56, no. 4, p. 253-260. [17] Ciulli, E., Bartilotta, I., Polacco, A., Manconi, S., Vela, D., Saverio, F., Paleotti, G. (2010). A model for scuffing prediction. Strojniški vestnik - Journal of Mechanical Engineering, vol. 56, no. 4, p. 231-238. [18] Li, J., Chen, W., Zhang, L., Wang, T. (2016). An improved quasi-dynamic analytical method to predict skidding in roller bearings under conditions of extremely light loads and whirling. Strojniški vestnik - Journal of Mechanical Engineering, vol. 62, no. 2, p. 86-94, D0I:10.5545/sv-jme.2015.2848. [19] Li, J., Wang, J., Zhou, G., Zheng, J., Hao, H. (2013). Study on failure mechanism of spatial lubrication harmonic reducer. Tribology, vol. 33, no. 1, p. 44-48, D0I:10.16078/j. tribology.01.011. (in Chinese) [20] Li, S. S., Shao, J., Yu, F., Ling, J., Chen, B. (2014). Research on modeling method for oil-film stiffness calculation of point contact based on EHL. Advanced Materials Research, vol. 945-949, p. 802-810, D0I:10.4028/www.scientific.net/ amr.945-949.802. 522 Xu, X - Fan, X. - Wei, P. - Yang, B. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 523-533 © 2020 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2020.6678 Original Scientific Paper Received for review: 2020-03-20 Received revised form: 2020-06-10 Accepted for publication: 2020-07-27 Dji amic Instability^ a W nd Turbine Blade Due to Large Deflections: An Ep erimental Validation Andres Lopez-Lopez1 - Jose Billerman Robles-Ocampo2 * - Perla Yazmin Sevilla-Camacho2 -Orlando Lastres-Danguillecourt1 - Jesús Muniz3 - Bianca Yadira Perez-Sariñana2 - Sergio de la Cruz2 1 University of Sciences and Arts of Chiapas, Mexico 2 Polytechnic University of Chiapas, Energy and Sustainability Academic Group, Mexico 3 National Autonomous University of Mexico Renewable Energy Institute, Mexico Wind turbine blades are designed to be thin and flexible elements. Because unstable dynamic behaviour can affect the life of the rotor, it is crucial to understand the instability of non-linear behaviour caused by large deflections. The present study undertakes both a stability analysis of the non-linear response and an experimental validation of a simplified model for a wind turbine blade based on a cantilever beam. The model is formulated taking into account large geometric deflections and assuming a Galerkin approach. The model is validated experimentally in a wind tunnel with aluminium beams of differing geometry. Analysis of the dynamic response using phase planes reveals that the degree of instability is related to the amplitude of the excitation and the stiffness characteristics. Keywords: non-linear model, large deflection, wind turbine blade, phase planes, instability Highlights • Slenderness and wind speed are parameters that have a significant influence on the stability of the dynamic behaviour of large wind turbine blades. • The non-linear dynamic model enables the dynamic behaviour of wind turbine blades to be obtained, taking into account the dynamic and gravitational loads. • The phase planes enable the level of stability of the dynamic behaviour ofHAWT blades to be detected qualitatively. • The nonlinearity of the model is linked to the large deflections present in the blade. 0 INTRODUCTION The strategy generally applied to obtain greater efficiency and increase the levels of power generated by a wind turbine involves the use of large diameter rotors and the manufacture of thinner, larger, lighter, and more flexible blades, resulting in a rotor diameter of 168 m [1], However, new blade designs cause greater sensitivity to dynamic excitations, which reduce wind turbine life span and efficiency, and increased vibrations and large deflections; these cause a high level of rotor instability [2], Furthermore, the blades are usually subject to random and complex mechanical stresses [3], In dynamic wind turbine blade models, elastic deformation is often considered in four directions, including flapwise (out-plane) and edgewise (in-plane or lead/lag) transverse deflections as the two main directions of vibration, as well as axial (longitudinal) and torsional (feather) deformations. As the blade is sufficiently stiff in both a torsional and longitudinal direction, the deformation occurring in these directions can usually be ignored. Moreover, as the cross-section of the blade is designed to make it stiffer in an edgewise direction than in a flapwise direction [4], researchers have focused more on studying the flapwise vibration of the blades [4] to [6], Recent studies conducted on dynamic blade models in a flapwise direction include that of Jokar et al. [7], who used the Hamiltonian principle to derive a non-linear partial differential equation for a wind turbine blade. Although they did linearize and simplify the non-linear model to find natural frequencies, their model did not include the effect of large deformations. However, for modern large-scale slender and flexible blades, consideration must also be given to the deformation caused by external loads which, in many cases, results in a non-linear geometric effect on the blade. Both this effect and the aeroelastic effects of extreme wind conditions are becoming more and more common with the increasing flexibility of wind turbine blades. Using the multi-body method, Xu et al. [8] built a blade analysis model that describes the geometric nonlinearity and complex geometry of the blades, for which they obtained accurate results. They found that the accuracy of their results increased as the number of rigid bodies increases, although they also observed that the calculation time increases exponentially. *Corr. Author's Address: Cuerpo Académico de Energía ySustentabilidad, Universidad Politécnica de Chiapas, Carretera Tuxtla Gutiérrez. - Portillo Zaragoza Km 21+500, Col. Las Brisas: Suchiapa, Chiapas. CP29150, México, jrobles@upchiapas.edu.mx 523 Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 523-533 Cruz et al. [9] presented the experimental validation of a simplified one-dimensional model of a cantilever beam formulated considering large geometric deformations and assuming a Galerkin approach. They aimed to ascertain the non-linear behaviour caused by the large deformations of thin structures, such as wind turbine blades. One-dimensional blade models have a low computational cost and have been used to efficiently model nonlinear vibration analysis for a wind turbine blade, for both in-plane [10] and out-plane motion [11], These models consider the geometric nonlinearity of blades, namely that large deformations cause a significant change in force, with a non-linear relationship between force and displacement. Gerstmayr and Irschik [12] represented large deformations using the elastic line approach. They used a cubic polynomial to represent the displacement; a similar approach is also applied in the present work. Maktouf et al. [13] undertook the modal and nonlinear dynamic analysis of a rotating wind turbine blade in order to study the effects of rotation speed on natural frequencies and displacements, finding that the natural frequencies and displacements of the blade increase when the angular rotation speed is increased. While the sensitivity of the blade's vibration characteristics to geometric nonlinearity was significant, the study did not consider aerodynamic effects. Liu [14] studied the effect of the aerodynamic load and its interaction with the structure under uniformly distributed wind flow conditions, determining an expression for the wind force per unit based on both average and fluctuating wind speed. Additionally, researchers have used both blade element momentum (BEM) theory and actuator disk theory in the development of aeroelastic models [15] to [17], Moreover, Akbiyik et al. [18] evaluated an air foil in a wind tunnel with a test section of 570 mm x 570 mm x 1000 mm, measuring and monitoring the response signal using an oscilloscope. While phase planes have been applied as a qualitative technique for analysing the dynamic stability of thin structures under considerable excitation amplitudes [19] and [20], phase planes should better pictorially represent the chaotic motion. The present study seeks to evaluate the unstable response of wind turbine blades with large deflections, assuming maximum flapwise deflections in the presence of extreme wind loads and the blade tip at the highest part of the rotor. The model is based on a one-dimensional cantilever beam Lagrangian formulation that considers the following: uniform wind distribution, an effective blade surface normal to wind incidence; non-linear curvature of the beam and internal viscous damping. Also, an experimental set based on a wind tunnel is proposed for model validation. The present research is of relevance for improving the design of wind turbine blades that present excessive vibrations and substantial deviations, as it provides an efficient tool for analysing the nonlinear response via phase planes. In addition, it offers a visual interpretation (qualitative behaviour) of the evolution of the instability under increasing wind speed. The main advantage of the method proposed is that, although it considers a similar shape to the blade and, solely deformation in a flapwise direction, the non-linear dynamic behaviour is described efficiently and with a low computational cost. Furthermore, the method determines the degree of blade stability in different conditions, ranging from low wind speeds to extreme loads. 1 DYNAMIC MODEL 1.1 Description of the Dynamic System The blades undergo variable tension-compression stresses during one rotation of a wind turbine rotor. These stresses are caused by the combined loads of the gravity and dynamic pressure of the fluid on the structure. Moreover, they depend on the angular position of the blade in the plane of rotation, namely that the angular position of the blade changes during its rotation in the plane of rotation, which also causes the stresses to change from tension to compression during each rotation and vice versa. However, when the blade is positioned vertically, the tip reaches the highest part of the wind turbine. In this instantaneous angular position, the bending moments that act on the blade are at a maximum because their interactions with wind speeds are more intense and provoke large blade deflections. Within this time frame, the dynamic behaviour of the rotor is strongly influenced by the large blade deflections. In contrast, the blade, which is joined to the wind turbine hub, resembles a cantilever beam. Under this analogy, the mathematical modelling for the blade is considered a thin blade rigidly embedded at one of its ends, which causes the other free end to bend laterally in the direction of the wind, as seen in Fig. 1. As the coordinate system (x. i j is perpendicular to the plane of rotation (x.z). the out-plane deflection is defined by £(v). the slope of deflection by 0(y). and the normal force of the wind on the rotation plane as f. 524 Lopez-Lopez, A. - RoMes-Ocampo, J.B. - Sevilla-Camacho, P.Y. - Lastres-Danguillecourt, 0,- Muniz, J. - Perez-Sariñana, B.Y. - de la Cruz, S. Strojnlskl vestnik - Journal of Mechanical Engineering 66(2020)9, 523-533 1.2 Large Blade Deflections In the dynamic system proposed, the wind force is assumed to be perpendicular to the plane of rotation, while the displacements out-plane are determined to be parallel to the wind direction, and the rotational displacements and torsional deformations are considered null [7], Thus, the movement of the blade in the coordinate plane (x, y) is defined as one-dimensional and restricts the movement of the deflection £(v). For the continuous vibratory system, the deflection of the blade at any point along its length is detennined by discretizing it into small one-dimensional beam-type elements by means of finite element method (FEM). This discretization enables the deflection £(v) to be defined as a function of the nodal displacements of each one-dimensional element, expressed in Eq. (1). \',(>') is called function form, while is nodal displacement and £(v) is the number of discretization nodes [21], n (1) The deflection of the blade is defined in Eq. (2) as a cubic polynomial dependent on the vertical position, where the coefficients g, relate the deflection '=/.. the Eqs. (19) and (20) are obtained. These equations describe the distribution of fluid over the blade. r=f, pg= 0. 1.6 Structural Damping Forces (19) (20) The structural damping force of the blades is caused by the dissipation of energy that occurs during deflection. This energy is not derived from a potential function and depends on the deflection speed. The dissipative energy of the blade ed is modelled using Eq. (21), where c is the structural damping coefficient [10]. (21) ED=\cèl The Rayleigh function models the internal dissipation forces of bodies with viscous damping. 526 Lopez-Lopez, A. - RoMes-Ocampo, J.B. - Sevilla-Camacho, P.Y. - Lastres-Danguillecourt, 0,- Muniz, J. - Perez-Sariñana, B.Y. - de la Cruz, S. Strojnlskl vestnik - Journal of Mechanical Engineering 66(2020)9, 523-533 D = -?*>> dgi (22) These structural damping forces are determined for +—tue +-/> +—-de +-A,(27) mL * mL' mL * mL' mL * where: Uf = EI (-1O4405L5 +540£d4L4 -3645fe'Ê + HOL1 +16200sw}, msw) = {y u/r(y) ak X> .(6) Step 4. The rough boundary interval assessment matrix is constructed as ©'= . Then, ©' is split into two matrices as follows. (7) (8) where cp' is the rough lower boundary matrix and rf is the rough upper boundary matrix. According to the gravity centre principle of triangular fuzzy number,

'm \. The synthesized attribute importance vector is as follows. co = | o>l, co2 (12) where co, = X(!>',' P ■ Then the attribute value matrix /;=i / with importance information is obtained as 538 Li, Y - Li, L Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 534-543 Next, we propose a modified TOPSIS by replacing Euclidean distance with relational vector distance based on set pair analysis theory. In traditional TOPSIS [16] and [17], the Euclidean distance from assessment object to the ideal solution is adopted for closeness calculation. However, the objects on the perpendicular bisector of two ideal solutions have the same closeness and cannot be distinguished. The principles of using proposed modified TOPSIS to realize the PSS scheme assessment are as follows. (1) An assessment object (i.e., a PSS scheme) and two ideal solutions (i.e., positive and negative ideal solutions) constitute two set pairs, respectively. (2) The set pair is then decomposed into multiple element pairs. (3) For each element pair, the sameness, contrariety and difference relationships are analysed to obtain the relationship between assessment object and ideal solution. Then the relational vector distance from the assessment object to the ideal solution is obtained. (4) According to the arranging rule of TOPSIS, the closeness of assessment object is calculated based on its relational vector distances. Finally, all assessment objects are arranged according to their closeness values. The detailed process of PSS scheme assessment based on modified TOPSIS is as follows. Step 1. Set pairs formed by assessment object and two ideal solutions are constructed. Based on y = [ v,, ]/ , the positive ideal solution is obtained as )' = [>■•, .>\.....ymJ and the negative ideal solution is obtained as y~ = [x\yl .....v~], where yl = max {>,, ..v.......y..\ and y- = mill¡.r. .. v. ......K = [vu,yv,.,„] represents PSS scheme i. Based on set pair analysis theory [18] to [20], )) and )' form a set pair, which is expressed by <)). y >. while )) and )' form a set pair which is expressed by <)].)' >. Through comparing the values, S. element pairs have tiny difference (sameness relationship), C,+ element pairs have huge difference (contrariety relationship) and I)j element pairs have not very obvious difference (difference relationship). Here, .S'/ +() l)i =m. Therefore, the relational degree of set pair <)). Y >. which characterizes the uncertain quantitative relationship between )) and )' . is expressed as follows: P,+ =s;A'+c;A"+d;Am. (13) where A', A", A'" indicate the sameness, contrariety and difference relationships respectively, and si = SI / m , cl = CI / m and d+t = £>,+ / m indicate the sameness, contrariety and difference coefficients, respectively. The relational vector of set pair <)). y > is as follows: n! =\s! ,c! ,dn (14) Similarly, the relational degree of set pair <)).y > is expressed as follows: p; = sTA' + c; A" + d; A'". (15) The relational vector of set pair <)).)' > is as follows: Vi = [sj,c;,d;]. (16) Step 2. Set pair is divided into several element pairs. Set pair <)). Y > consists of m element pairs • .>•: •.>• ••• .>; •.>•; •.....• •• For element pair < >', ,. y, >, its relational degree can be expressed as follows: Pu = + C!A" + d!tA"\ (17) where A', A", A'" indicate the sameness, contrariety and difference relationships respectively, and ,v( (, cl, and i/ ; indicate the sameness, contrariety and difference coefficients of element pair <>',,.>', > respectively. If y„ = y, . then ,v( ( = 1. cn =0 and i/ ; = 0. If v„ = v, . then ,v(( =0. cn = 1 and 4 = 0. If yl'„)/(}', -}', ). Therefore, the relational degree of < v,.,,yl >is: m m m IX, IX, Pi -A' + —-A" + —-A'". m m m Similarly, for set pair <))' > fonned by )', and )' . it consists of m element pairs < V,,i-.vr > - < yU2,v~2 ■.....• „..(•„, > . The relational degree of element pair is expressed as follows: />,. -c:r\" ■ d:rV. (18) where su, cu and du indicate the sameness, contrariety and difference coefficients of element pair respectively. If yu = v( . then, then su = 0, cu = 1 and dn = 0. If yu = yl, then su = 1, cu = 0 and dl, =0. If yl is: P, = HSU IX 1=!-A' + _i=i-A" + —-A'". m m m Step 3. Set pair formed by assessment object and itself is constructed. Yj and Yf can form a set pair )..). . Because a set is the same as itself, the relational degree of set pair is p] = 1 • A' + 0• A" + 0• A'". Therefore, its relational vector is u] = [1,0,0]. Step 4. The relational vector distance from the assessment object to the ideal solution is calculated. Because the relational vector of set pairs <)). Y > and <}■,}■> are ^ =[s%c%dj] and u,' =|UU)|. respectively, the relational vector distance from assessment object )) to ideal solution )' is calculated as follows. r; = V«-i)2+(C,+)2+«)2. (19) Similarly, the relational vector distance from assessment object )) to ideal solution )' is calculated as follows. r; =yl(s;-if+(c;)1 + (d;)1. (20) Step 5. The closeness of assessment object to positive ideal solution is calculated and adopted as the arranging basis. By replacing the Euclidean distance with the relational vector distance, TOPSIS is modified. For )]. its closeness to the positive ideal solution is calculated as follows. X, =- r. +r, (21) According to the arranging rule of TOPSIS, all assessment objects are arranged according to their closeness values. The PSS scheme with the biggest closeness value is the best one, and the PSS scheme selection optimization is achieved. 2 CASE STUDY In order to enhance the competitiveness of the enterprise, to meet the growing personalized needs of consumers, and to improve the efficiency of product sales, an air purifier manufacturing enterprise began to plan and build an air purification PSS (AirP-PSS). With the AirP-PSS, the enterprise sells the promised level of air purification effectiveness instead of selling air purifiers. In the design phase, several AirP-PSS schemes with different characteristics are provided by PSS design engineers. To achieve AirP-PSS scheme selection optimization, the enterprise will carry out AirP-PSS scheme assessment and determine the best one. There are five alternate AirP-PSS schemes: intelligent series (1\). energy saving series (,P2), portable series (P3), sterilization series (P4) and humidification series (P5). Based on the membership function of trapezoid fuzzy number, natural number N can be converted to trapezoid fuzzy number N as follows. (1,1,-,2), N = l 2 N-l TV + 1 (N-l-,-,TV +1), 1 < TV < 9. N = 17 (8, —,9,9), 2 N = 9 According to the arithmetic rules of the trapezoid fuzzy number, the commonly used nine-level scale assessment comments and values are fuzzified to get the corresponding trapezoid fuzzy number, as shown in Table 1. Table 1. Comments and fuzzifieded values of nine-level scale assessment Comment Fuzzed value Extremely superior 9/1 =(4,17/3,9,9) Strongly superior 8/2 = (7/3,3,17/3,9) Obviously superior 7/3 = (3/2,13/7,3,4) Slightly superior 6/4 = (1,11/9,13/7,7/3) Equal 5/5 =(1,1,1,1) Slightly inferior 4/6 = (3/7,7/13,9/11,1) Obviously inferior 3/7 =(1/4,1/3,7/13,2/3) Strongly inferior 2/8 = (1/9,3/17,1/3,3/7) Extremely inferior 1/9 = (1/9, 1/9,3/17, 1/4) Based on the digital twin-driven PSS scheme selection optimization framework shown in Fig. 1, the AirP-PSS scheme assessment is done after the experts master the digital twin information comprehensively. The detailed process is as follows. Three experts participate in the AirP-PSS scheme assessment. For attribute .!,. the fuzzy reciprocal assessment matrices given by three experts are i?u, E2'1 and E3-1 as shown in Tables 2, 3, and 4. 540 Li, Y - Li, L Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 534-543 Table 2. Fuzzy reciprocal assessment El,x Pi Pi P3 P4 P5 Pi 5/5 5/5 6/4 5/5 6/4 P2 5/5 5/5 6/4 5/5 6/4 4/6 4/6 5/5 4/6 5/5 P4 5/5 5/5 6/4 5/5 6/4 Ps 4/6 4/6 5/5 4/6 5/5 Table 3. Fuzzy reciprocal assessment E2,1 Pi Pi PJ P 4 P5 Px 5/5 7/3 5/5 6/4 7/3 P2 3/7 5/5 3/7 4/6 5/5 P3 5/5 7/3 5/5 6/4 7/3 P4 4/6 6/4 4/6 5/5 6/4 P s 3/7 5/5 3/7 4/6 5/5 Table 4. Fuzzy reciprocal assessment E^'1 Pi Pi P3 P4 P5 Px 5/5 6/4 6/4 7/3 6/4 P2 4/6 5/5 6/4 6/4 5/5 P3 4/6 4/6 5/5 6/4 4/6 P4 3/7 4/6 4/6 5/5 4/6 Ps 4/6 5/5 6/4 6/4 5/5 El-\ E2'1 and it3,1 are all qualified by consistency inspection. The group assessment matrix is constructed as ^=KvL,inwhich su According to Eqs. (5) and (6), the rough boundary interval of canbe obtained as: RBI{e\2) =[(1.0556, 1.1570, 1.4603, 1.7037),(1.3056, 1.5855,2.4603, 3.2037)]. After calculating the rough boundary intervals of other elements in El = \el. 1 , the rough L 'J J5x5 boundary interval assessment matrix is constructed as 0' = \RBI(e)t) 1 . Then, ©' is split into rough L -J J5x5 . r . -| lower boundary matrix

'J5x5 Table 5. Rough lower boundary matrix lz/1. Table 8. Attribute value matrix with Importance Information Y=[y» L a! a 2 a 3 a4 a5 py 0.1030 0.0494 0.0714 0.0870 0.0932 p?. 0.0697 0.1190 0.1077 0.0401 0.1057 pi 0.0771 0.1429 0.0796 0.0248 0.0777 p4 0.0677 0.0687 0.0445 0.2017 0.0601 p5 0.0603 0.1552 0.0294 0.1236 0.1003 Next, the proposed modified TOPSIS is adopted for AirP-PSS scheme assessment to achieve scheme optimization. According to )' = [>■',,]_ _ ■ the positive ideal solution is obtained as y+ = [0J030, 0.1552, 0.1077, 0.2017, 0.1057] and the negative ideal solution is obtained as y~ = [0.0603, 0.0494, 0.0294, 0.0248, 0.0601], Set pair <}'1,}'+> consists of five element pairs < vu,yl >,< y1X •.....• v. ...r. >. For element pair < >|3. >', >, the relational degree of < • >3 > is obtained as = 0-A'+ 0-A"+ 0.5364-A". The relational degrees of other element pairs in can be obtained similarly. Then, the relational vector of set pair <}'1,}'+> is obtained as =[0.2000,0.2000,0.3229]. In the same manner, the relational vector of set pair <}'b )' > is obtained as = [0.2000,0.2000,0.2771]. The relational vector of the set pair <)',. )\> is jui = [1,0,0], Based on the relational vectors of <}'1,}'+> , <}'!,}'"> and <}'!,}'!>, we obtain that the relational vector distances from )', to )' and from )', to y~ are r, =0.8856 and r, =0.8699. Therefore, the closeness of Y^ to the positive ideal solution is calculated as Xx =0.4955. Similarly, the closeness values of other assessment objects to the positive ideal solution are calculated as A2 = 0.6463, ¿3 = 0.4371, A4=0.5356 and = 0.4060. According to the arranging rule of TOPSIS, five AirP-PSS schemes are arranged as /J2>7J4>7J, >P3>P5. p2 is the best AirP-PSS scheme, and the scheme selection optimization is achieved. 3 CONCLUSION Because the influencing factors of PSS scheme selection optimization have dynamic multidimensional characteristics and a complex coupling relationship in the optimization process, static assessment and selection optimization are with greater restrictions and unreasonable. Aiming at this problem, a digital twin driven framework to enhance PSS scheme selection optimization is presented, and its feasibility is verified by a case of air purification PSS scheme selection optimization. This paper explores the application of digital twin-driven frameworks in the optimization of PSS scheme selection. Although the work enhances PSS scheme selection optimization, some limitations call for future research. For example, the assessment data modelling of digital twin-enabled PSS scheme selection could be considered to provide a detailed information-analysing methodology. More data and cases should be used to improve the approach's feasibility and practicality. 4 ACKNOWLEDGEMENTS This work is supported by the Curriculum Reform Project of Teacher Education in Henan Province (B054), Ningxia Natural Science Foundation (2020AAC03202) and the Third Batch of Ningxia Youth Talents Supporting Program (TGJC2018048). The authors also would thank Hongxia Sun, Bingbing Lei, Chuntao Zhang and Chunlei Mao for their contributions to this work. 5 REFERENCES [1] Chiu, M.-C, Chu, C.-Y, Kuo, T.-C. (2019). Product service system transition method: building firm's core competence of enterprise. International Journal of Production Research, vol. 57, no 20, p. 6452-6472, D0l:10.1080/00207543.2019.156 6670 [2] Kim, S, Christiaans, H, Baek, J.-S. (2019). Smart Homes as product-service systems: Two focal areas or developing competitive smart home appliances. Service Science, vol. 11, no 4, p. 292-310, D0l:10.1287/serv.2019.0248 [3] Liu, Z, Ming X. (2019). A framework with revised rough-DEMATEL to capture and evaluate requirements for smart industrial product-service system of systems. International Journal of Production Research, vol. 57, no. 22, p. 7104-7122, D01:10.1080/00207543.2019.1577566 542 Li, Y - Li, L Strojniski vestnik - Journal of Mechanical Engineering 66(2020)9, 534-543 [4] Meier, H., Roy, R., Seliger, G. (2010). Industrial Product-Service systems-IPS2. CIRPAnnals, vol. 59, no. 2, p. 607-627, DOI:10.1016/j.cirp.2010.05.004. [5] Sun, H., Wang, Z., Zhang, Y., Chang, Z., Mo, R., Liu, Y. (2012). Evaluation method of product-service performance. International Journal of Computer Integrated Manufacturing, vol. 25, no. 2, p. 150-157, DOI:10.1080/095119 2X.2011.627946. [6] Yoon, B., Kim, S., Rhee, J. (2012). An evaluation method for designing a new product-service system. Expert Systems with Applications, vol. 39, no. 3, p. 3100-3108, DOI:10.1016/j. eswa.2011.08.173. [7] An, X., Niu, C., Xue D., Mu, G., Lin, X. (2016). Optimal selection method for product service system scheme based on variable granule weight and group decision-making. Computer Integrated Manufacturing Systems, vol. 22, no. 1, p. 155-165, DOI: 10.13196/j.cims.2016.01.015. (in Chinese) [8] Lan, S., Zhang, H., Zhong, R.Y., Huang, G.Q. (2016). A customer satisfaction evaluation model for logistics services using fuzzy analytic hierarchy process. Industrial Management & Data Systems, vol. 116, no. 5, p. 1024-1042, DOI:10.1108/IMDS-09-2015-0389. [9] Chou, C.-J., Chen, C.-W., Conley, C. (2015). An approach to assessing sustainable product-service systems. Journal of Cleaner Production, vol. 86, p. 277-284, DOI:10.1016/j. jclepro.2014.08.059. [10] Liu, Q., Zhang, H., Leng J., Chen X. (2019). Digital twin-driven rapid individualised designing of automated flow-shop manufacturing system. International Journal of Production Research, vol. 57, no. 12, p. 3903-3919, DOI:10.1080/0020 7543.2018.1471243. [11] Zhang, K., Qu, T., Zhou, D., Jiang, H., Lin, Y., Li, P., Guo, H., Liu, Y., Li, C., Huang G.Q. (2020). Digital twin-based opti-state control method for a synchronized production operation system. Robotics and Computer-Integrated Manufacturing, vol. 63, p. 1-15, DOI:10.1016/j.rcim.2019.101892. [12] Leng, J., Zhang, H., Yan, D., Liu, Q., Chen, X., Zhang, D. (2019). Digital twin-driven manufacturing cyber-physical system for parallel controlling of smart workshop. Journal of Ambient Intelligence and Humanized Computing. vol. 10, p. 11551166, DOI:10.1007/s12652-018-0881-5. [13] LI, L.-H., Hang, J.-C., Gao, Y, Mu, C.-Y. (2017). Using an Integrated group decision method based on SVM, TFN-RS-AHP, and TOPSIS-CD for cloud service supplier selection. Mathematical Problems in Engineering, vol. 2017, p. 1-15, DOI:10.1155/2017/3143502. [14] Li, L.-H., Hang, J.-C., Sun, H.-X, Wang, L. (2017). A conjunctive multiple-criteria decision-making approach for cloud service supplier selection of manufacturing enterprise. Advances in Mechanical Engineering, vol. 9, no. 3, pp. 1-14, DOI:10.1177/1687814016686264. [15] Geng, X., Xu S., Li Y., You, X. (2015). Product service system concept evaluating approach considering combinatorial effect between product and service. Computer Integrated Manufacturing Systems, vol. 21, no. 10, p. 2798-2806, DOI:10.13196/j.cims.2015.10.029. (in Chinese) [16] Yuvaraj, T., Suresh, P. (2019). Analysis of EDM process parameters on inconel 718 using the grey-Taguchi and Topsis methods. Strojniški vestnik - Journal of Mechanical Engineering, vol. 65, no. 10, p. 557-564, DOI:10.5545/sv-jme.2019.6194. [17] Alhabo, M., Zhang, L. (2018). Multi-Criteria Handover Using Modified Weighted TOPSIS Methods for Heterogeneous Networks. IEEE Access, vol. 6, p. 40547-40558, DOI:10.1109/ ACCESS.2018.2846045. [18] Li, L.-H., Mao, C.-L., Sun, H.-X., Yuan, Y.-P., Lei, B.-B. (2020). Digital twin driven green performance evaluation methodology of intelligent manufacturing: Hybrid model based on fuzzy rough-sets AHP, multistage weight synthesis, and PROMETHEE II. complexity, vol 2020, p. 1-24, DOI:10.1155/2020/3853925. [19] Yang, M., Zhang, L., Cui, Y., Zhou, Y., Chen, Y., Yan, G. (2020). Investigating the wind power smoothing effect using set pair analysis. IEEE Transactions on Sustainable Energy, vol. 11, no. 3, p. 1161-1172, DOI:10.1109/TSTE.2019.2920255. [20] Garg, H., Kumar, K. (2020). A novel possibility measure to interval-valued intuitionistic fuzzy set using connection number of set pair analysis and its applications. Neural Computing & Applications, vol. 32, no. 8, p. 3337-3348, DOI:10.1007/s00521-019-04291-w. Enhancing the Optimization of the Selection of a Product Service System Scheme: A Digital Twin-Driven Framework 543 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9,544-553 © 2020 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2020.6776 Original Scientific Paper Received for review: 2020-05-15 Received revised form: 2020-08-19 Accepted for publication: 2020-08-20 Numerical Anaty is with Experimental Validation of Single-Phase Fluid Flow in a Dimple Pattern Heat Ek hanger Channel Urban Mocnik12 - Bogdan Blagojevic3 - Simon Muliic14 * JUniversity of Novo mesto, Faculty of Mechanical Engineering, Slovenia 2Danfoss Trata d.o.o., Slovenia 3Plinovodi d.o.o., Slovenia 4SIMUTEH s.p., Slovenia A plate heat exchanger with a dimple pattern heat plate has a large number of dimples. The shape of dimples defines the characteristics of the plate heat exchanger. Although such heat exchangers have become increasingly popular due to their beneficial characteristics, knowledge of the flow characteristics in such kind of channel is poor. A good knowledge of the flow conditions inside of such channel is crucial for the successful and efficient development of new products. In this paper single-phase water flow in dimple pattern plate heat exchanger was investigated with application of computational fluid dynamics and laboratory experiments. Numerical analysis was performed with two turbulence models, Realizable k-e with enhanced wall treatment function and k-co SSI The first predicts a slightly smaller pressure drop and the second slightly larger compared to the results of laboratory measurements. Our research found that despite the relatively low velocity of the fluid, turbulent flow occurs in the channel due to its shape. We a/so found that there are two different flow regimes in the micro plate heat exchanger channel. The first regime is the regime that dominates the heat transfer, and the second is the regime where a recirculation zone appears behind the brazing point, which reduces the surface for heat transfer. The size of the second regime does not change significantly with the velocity of the fluid in the volume considered. Keywords: heat exchanger, dimple pattern, pressure drop, computational fluid dynamics, turbulence models Highlights • We carried out measurements of the pressure drop in channel formed by two dimple pattern heat plates of heat exchanger by a novel measuring method which allows direct comparison. • We presented a comprehensive analysis of the flow state in a plate heat exchanger with a dimple pattern heat plate structure. • We addressed the differences in CFD results between two major turbulence models. • We proposed the computational mesh structure to achieve a fair ratio between accuracy and computational resources. 0 INTRODUCTION A heat exchanger is a device that allows the transfer of heat between two or more liquids of different temperatures. Heat exchangers are useful in many systems in the industry and they are also widely represented in district heating systems. A plate heat exchanger is a type of heat exchanger in which hot and cold fluids are separated by a thin wall. These walls have a characteristic corrugated shape and are called heat plates. There can be several hundred thin walls in one plate heat exchanger where each of the other heat plate faces the opposite direction, forming a channel in which there is water flow. This creates a number of channels, half of which are connected to the primary and half to the secondary plate. The heat plates are made of stainless steel. Bonding material, most often copper, leaks into the edges of the product due to the capillary effect and forms a seal in a vacuum furnace. At the same time, copper also flows into the junctions of the plates, where it forms the joints referred as brazing points. The design of the heat plates determines the characteristics of the heat exchanger. Worldwide, the most popular type of plate is so called fishbone or chevron shape. These plates are shaped as shown in Fig. la. The angle and depth of the corrugations determine the characteristics of the heat exchanger. Products with such a plate are simply referred to as brazed plate heat exchangers (BPHE). This type of heat exchangers, with some improvements, have been in use since the 1970s [1], a) b} Fig. 1. a) The fish bone pattern heat plate, and b) dimple pattern heat plate Recently however, a new heat plate design has emerged that significantly improves the characteristics 544 *Corr. Author's Address: University of Novo mesto. Faculty of Mechanical Engineering, Na Loko 2, 8000 Novo mesto, Slovenia, simon.muhic@fs-unm.si Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 544-553 of heat exchangers. New structure of the plate does no longer have straight corrugations, but a large number of dimples, which is why it is called a dimple pattern. The shape of dimples defines the characteristics of this type of heat exchanger. The dimple pattern heat plate is shown in Fig. lb. While there is widely available literature about the fluid flow through a duct in the field of BPHE, there is no available literature of dimple pattern type heat exchangers. The flow of liquid at a micro-level is very difficult to observe in the heat exchanger due to small channels and the diverse structure [2], Some have succeeded to some extent, for example Focke and Knibbe [3], But with the spread of computational fluid dynamics (CFD) researchers finally got a useful tool to investigate fluid flow conditions inside plate heat exchangers. In the beginning, they first studied flows in two dimensions, just like Metwally and Manglik [4], They investigated the distribution of flow and temperature, friction factor and the Nusselt number in the sinusoidal fraction of the heat exchanger channel. Kanaris et al. [5] has used more geometry and observed similar phenomena. With the increasing capacity of computers, scientists have been able to start studying the entire BPHE channel. Han et al. [6] analysed a few consecutive channels and observed the heat transfer and flow distribution. Their numerical results deviated from the laboratory at some points up to 35 %, but some uncertainty was attributed also to the laboratory experiment. Gherasim et al. [7] used a computer-modelled model of two BPHE plates in their calculations, but excluded the geometry of the distribution channel and instead predicted a uniform distribution of fluid flow velocity at the inlet and outlet in the channel between the plates. They calculated the case with different turbulence models at different flow velocities. The closest results compared to experimental results were achieved with Realizable k-e turbulence model with Non-equilibrium wall functions and an average dimensionless wall distance v+ between 2.8 (for Re = 400) and 11.25 (for Re = 3000). In general, their numerical results differ from the experimental results by 10 % to 16 %. Gullapalli and Sunden [8] performed simulations of fluid flow throughout the channel of a plate heat exchanger with several different corrugation angles, and added an inlet and outlet connection to create similar inlet or outlet conditions as in the laboratory experiment. They used the LRR-IP (Reynolds Stress) turbulence model. They performed several calculations with different boundary conditions on the wall of the plate when calculating heat transfer (e.g. constant heat flux, constant wall temperature, etc.). The results of the heat transfer calculations were 20 % to 30 % and the pressure differences were 10 % to 35 % below the values from the laboratory experiment. Nonetheless, they found CFD to be a useful tool in the relative comparison of the various shapes of the heat plates, for detection of possible poorly constructed geometric details of corrugation, and the determination of velocity profiles inside of the channels. Tiwari et al. [9] analysed a case using CFD tools, with two channels in which a fluid with a nanoparticle (nanofluid) flowed counter-current. They also used a - turbulence model in the study. They also carried out a laboratory experiment with the same configuration of the plate heat exchanger elements as they used in the numerical calculation, with the latter results differing from the experimental results by a maximum of 3.75 %. Lee and Lee [10] performed non-stationary CFD analyses using the large eddy simulation (LES) approach on a small part of a fish bone heat plate heat exchanger. They found oscillations in the flow at a turbulent flow regime, which increase the heat transfer and also increase the pressure drop. Sarraf et al. [11] carried out laboratory and CFD analyses of two-channel sample. They compared numerical results with the results obtained with a thermal camera, where they recorded the exterior of the plate pack during laboratory experiments. The temperature distribution difference between CFD results and experimental results was within 10 %. They obtained the results with a k-e turbulence model. They also tried to analyse the channel flow with a k-m turbulence model, but were not successful due to convergence problems. They distinguished two types of fluid flow, „helical" and „cross" flow, that coexists in the channels formed by the fish bone patterned heat plates. Helical flow appears in direction from the inlet to outlet passing grooves while cross flow appears along the groove. Which flow prevails depends on the mass flow and the angle between the grooves on the heat plate. Piper et al. [12] presented a numerical study of a pillow plate heat exchanger using a plate with a periodic pillow-like pattern. They found that fluid flow in such a channel is characterized by two distinct areas - the meander core, which is crucial for heat transfer, and the recirculation zone that occurs behind the welded joint and reduces the efficiency of the heat exchanger. They suggested several possible solutions to improve the efficiency of such heat exchanger, which could increase efficiency by up to 37 %. In numerical analysis, they used a domain derived from a transformation simulation. The domain was a small part of the entire channel which, however, repeats periodically throughout the volume. They simplified Numerical Analysis with Experimental Validation of Single-Phase Fluid Flow in a Dimple Pattern Heat Exchanger Channel 545 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 544-553 the calculation with all possible symmetries. They used the k - s turbulence model and constant heat flux as a boundary condition on the domain wall. In a study presented by Yogesh et al. [13] investigating the influence of pipe geometry in a tubular heat exchanger with cooling fins, authors used k - m turbulence model instead of the commonly used k - s turbulence model to obtain results that were more consistent with their reference. For the k - m turbulence model, the results were 5.2 % and 6.6 % lower than the reference ones, unlike the k - s turbulence model results where the difference was 26.9 % and 44.8 %. Al zahrani et al. [14] performed a numerical analysis of a plate heat exchanger with five plates, two channels on the warm side and two on the cold side. They also used a k - s turbulence model. They compared the numerical results with the experimental ones obtained by Muley and Manglik [15]. The deviations were <10 %. As mentioned above, despite the fact that numerous investigations were conducted on chevron plates heat exchanger, no literature on dimple pattern heat exchanger can be found. In our study, we have analysed the fluid flow in a channel formed by dimple pattern heat plate numerically. The obtained results were validated with laboratory results. These were obtained by newly developed measuring method, which allows direct comparison. 1 NUMERICAL ANALYSIS OF FLUID FLOW 1.1 Channel Geometry Modeling The three- dimensional (3D) model of the heat plate was obtained by 3D scanning of the real plate. One-sixteenth of the cell was cut from the scanned 3D model. With the microscope measuring system, we also determined the appropriate shape of the brazing point and transferred it to the 3D model. The whole-cell model was created by mirroring over the corresponding edges. Symmetric cell was then multiplied by mirroring to obtain the desired numerical domain. Fig. 2. Cell model (red) and part of the channel; green arrow indicates flow direction, domain inlet is darker, cell inlet/outlet is dark green The geometry presented in Fig. 3 was made from the base cell, which served as a basis for modelling narrower or wider domains. The number of cells in basic geometry is 28. When referred on e.g. 14*2 cells means 14 successive cells in direction and 2 parallel cells in x direction. During the study, it was shown that due to the uniform velocity distribution at the inlet surfaces (darker in Fig. 2), the fluid flow in the first few cells was not yet developed to such an extent that it could be considered as representative. Similar things happen in the last few cells. For the aforementioned reasons, rather than for longitudinal periodic boundary conditions, we decided to extend the domain to the maximum number of consecutive cells in which we still obtained useful results using various turbulence model, that is 14 consecutive cells. For longer domains, this was not possible. There is an extension at the end of the geometric model to prevent backflow at the outlet surface. The analysis of the impact of width of the numerical domain on the simulation result was performed on domains that also have a length of 14 cells, the widths of which are presented in Fig. 3. Domains with multiple parallel cells (2, 3 and 6) have the same amount of inlet and outlet surfaces, each outlet has an extension. When we refer to water velocity in this article, we refer to the average velocity of water through the largest cross-section of the cell, which is equal to the area of one inlet or outlet surface of the domain. mY_i--1--.--1 WHHW I 1 2 3 10 11 12 13 14 Fig. 3. The basic geometry of the domain with numbering of a) narrower or wider domain, and b) successive cell 1.2 CFD Simulation The CFD simulation was performed using the commercial Ansys Fluent 19.2 CFD software, which operates on the finite volume method. The turbulent flow was modelled using Reynolds-averaged Navier-Stokes equations (RANS). Reynolds stresses were calculated using Realizable k-s turbulence model with Enhanced wall treatment (EWT) wall function (hereinafter k-s) and k-m SST turbulence models 546 Močnik, U. - Blagojevič, B. - Muhič, S. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 544-553 (hereinafter k-m) [16], The coupling between velocity and pressure was performed using a coupled scheme. During calculations, we monitored the change in significant quantities such as: residuals, area average velocity through the area between the 12th and 13th cell (Fig. 3b), area average pressure between the 2nd and 3rd cell (Fig. 3b), the maximal v+ on the wall and velocities in left and right points, (Fig. 3b, green points between 10th and 11th cell). The basic relations for the motion of incompressible fluid are the continuity and momentum equations [16]: ^ + V\pu) = Sm, (1) at d - — - — (pi!) + V • (puu) = -Vp + V • (f) + pg + F. (2) The dimensionless distance between the mesh node closest to the wall and the wall is presented by the equation [16]: y ~ • where "'■(p)'- <4) The domain was divided into seven elements. Planes which divided the elements were positions on which the results were monitored. This is shown in Fig. 3b. The lines in the figure represent the planes on which the area average pressure was monitored. The difference in pressures on one element, on two cells, was then calculated by subtracting the pressure at the beginning of the element from the pressure at the end of the element. Thus, we obtained a pressure difference (pressure drop) on all individual domain elements. These results of the pressure difference were compared with those obtained through laboratory experiments. 1.3 Mesh Generation The mesh was created in Ansys Fluent (with Fluent meshing) software enviromnent. During research, the flow conditions inside the described geometry were proved to be challenging. Knowing this, we implemented a polyhedron mesh. Tests have shown that calculations with polyhedral meshes need up to four times fewer elements, half the computational memory and one-tenth to one-fifth of the computational time compared to tetrahedral meshes to achieve the same accuracy [17], In addition to the aforementioned properties, during the research, it was also found that a convergent solution cannot be easily achieved by a tetrahedron mesh. Taking this into account we conducted a mesh density independence study. We produced ten different mesh densities, which are presented in Table 1 (number of finite volumes is valid for basic domain presented on Fig. 3b). In Table 1 there is also data for maximal >• for each mesh density at different velocities. Table 1. Mesh properties for basic domain Mesh No. of finite Maximal v*" at different velocities volumes 0.2 m/s 0.45 m/s 0.7 m/s M1 107126 1.3 2.3 3.4 M2 143829 1.4 2.5 3.4 M3 183973 1.3 2.4 3.3 M4 270753 1.0 1.9 2.6 M5 533612 0.5 1.0 1.3 M6 743587 0.6 0.9 1.2 M7 1158886 0.5 0.8 1.1 M8 1464905 0.6 0.9 1.2 M9 2408576 0.5 0.8 1.1 M10 5372272 0.5 0.7 1.0 Average v^ for the mesh M4 using the k-e based turbulence model is 1.01 at the velocity of 0.7 m/s. Using k—co based turbulence model, y* is on the same level. The ratio of maximal to average v+ is approximately 2.6 and is maintained in all cases. 1.4 Boundary Conditions In practice, the water velocity in the channel is between 0.2 m/s and 0.7 m/s. For the calculations, where we checked the adequacy of the mesh and the computational model and thus the consistency of the calculated results with measurements, we used velocities between 0.2 m/s and 0.7 m/s for the inlet boundary condition. The considered flow was singlephase, incompressible, stationary, three-dimensional, and turbulent. Working fluid is water with constant physical properties (water at temperature 15.2 °C). Natural convection and radiation are neglected. The symmetry boundary conditions were determined for the 14x1/4, 14xl/2H, 14xl/2V and 14x1 domains. The 14x2, 14x3, and 14x6 domains do not have symmetry boundary conditions applied (if 14x1 cells means domain with 14 successive cells with extension as presented in Fig. 3b, 14x3 cells would mean that there are three main domains bonded in parallel). The surfaces on which the symmetry boundary conditions Numerical Analysis with Experimental Validation of Single-Phase Fluid Flow in a Dimple Pattern Heat Exchanger Channel 547 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 544-553 were determined are indicated in green in Fig. 3 a under marks for 1/4, 1/2H, 1/2V and 1 where results were obtained with, as well as without symmetry. 2 EXPERIMENTAL ANALYSIS OF FLUID FLOW The test bench shown in Fig. 4 consists of a frequency-controlled pump that drives the water in a closed circuit, an expansion vessel that maintains a uniform water pressure, a flowmeter, a water tank, a cooling heat exchanger to maintain the temperature of the water during the measurement and a large vessels (LV) assembly with sample and gauge pressure and pressure drop sensor. The LV assembly consists of two large stainless-steel vessels. These are built in the same way, with exception that one has an opening for installing the temperature sensor. In use, however, they are facing each other so that the test sample can be clamped to them, as shown in Figs. 4 and 5. At the inlet and/or outlet, respectively, the flow rectifier is installed with the purpose to direct the flow of water towards the walls of the vessels, thereby mixing the water inside and preventing a portion of the flow from having direct access to the inlet surfaces of the test sample, which could cause uneven flow to the inlet surfaces of the test sample. At the same time, water with a homogeneous temperature is obtained throughout the whole volume. The temperature sensor is positioned in the centre of the vessel and measures the temperature at which water enters the sample. The pressure is measured on the wall of the vessel as shown in Fig. 4. The pressure drop sensor is mounted on the edge of the vessels where no vortices are expected, in order to provide a stable pressure measurement. Because the fluid flow is extremely slow in the vessels, the pressure field is also homogeneous. This way, the same inlet conditions are ensured on all inlet surfaces of the sample. The measurements were performed at discrete points at 0.2 m/s, 0.3 m/s, 0.4 m/s, 0.45 m/s, 0.5 m/s, 0.6 m/s, and 0.7 m/s average water velocity through the inlet surfaces. At each point, a steady state was established and the data of all measured values was simultaneously recorded for 30 s at a frequency of 5 Hz. The result is the average of all measured quantities. We also carried out measurements of the pressure drop at a different water temperatures, which were 15.2 °C, 16 °C, 18 °C, 20 °C and 25 °C with the maximum standard deviation of temperature 0.19 °C. Measurements of pressure drop at 15.2 °C were performed on three samples of 19 cells (inlet surfaces) of width and 31, 24 and 16 consecutive cells long, respectively. 1- Test sample, 2 - Big chamber, 3 - Deflector, 4 - Pressure drop sensor, 5 - Gauge pressure sensor, 6 - Flowmeter, 7 - Temperature sensor, 8 - Pump, 9 - Cooling heat exchanger, 10 - Water tank, 11 - Expansion vessel Fig. 4. Test bench scheme The measurement uncertainty of the results for flow is 1 %, the pressure drop is 1 % and the temperature is 0.3 °C. Fig. 5. Photo of the large vessel assembly with the test sample during the measurement 3 RESULTS AND DISCUSSION As described, we conducted a mesh density independence study for 10 different meshes with 3 boundary conditions for two turbulence models (k- s and k- ®). Fig. 6 shows the dependence of the simulation results on the mesh density for the area average pressure drop between 9th and 10th cell marked on Fig. 3b. The mesh number 4 (M4) with 270753 finite volumes is used as a reference in this graph. The entire XB12L heat plate has about 1500 dimple cells, which means that it is necessary to economically choose the density of the mesh with which to calculate solutions that will cover conditions in the whole channel in the future, or even in several channels of the considered heat exchanger. In our case, we decided to perform the calculations with the density of M4 shown in Fig. 6, which means about 10000 finite volumes per dimple cell. This is a manageable number of mesh elements, and the difference from the mesh independent solution is a 548 Močnik, U. - Blagojevič, B. - Muhič, S. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 544-553 reasonable 3 % to 4 % for both turbulence models and for all velocities in the channel. The only difference is the k - a> solution at 0.7 m/s, which no longer follows other solutions at mesh densities >6. Fig. 6. Dependence of numerical results on mesh density Fig. 7 shows the dependence of the domain width on the results of pressure drop (Fig. 3). As a reference, in this graph the k- s turbulence model uses a 14x3 domain and the k- a> turbulence model a 14*1 domain. In Fig. 3, the green line indicates the surfaces where symmetric boundary conditions were used. 1/4 1/2H 1/2V Mxlsvmm 14» 1 14*2 14x3 14*6 Fig. 7. Dependence of numerical result on domain width Numerical results analysis showed that the use of k - s turbulence model requires the use of a domain with three parallel cells in our case 14*3 domain (Fig. 3a, no. 3), and for k - a> turbulence model, one cell type is sufficient (Fig. 3a, no. 1). From the graph in Fig. 7, it can be seen that the predicted results for the domains 1/4 and 1/2H (that is for those two cells in which the cells are cut horizontally) are showing pressure drop which is significantly too high. Slightly better results are obtained for vertically sectioned cells, the 1/2V domain, where the results are still high for k - s and closer to the measurement results for k- a> turbulence model. It is interesting to note that the results for the "14x1 symm" domain, where symmetry was used as the wall boundary condition, are worse than for the same domain where we did not use symmetry boundary condition. The domain geometry in this part was defined as a wall on which, like everywhere else on the wall, a boundary layer mesh was built. With the k- e, the result is slightly different with the expansion of the domain, but then stabilizes. The 14*3 domain was chosen as the most appropriate. However, with the k- a> turbulence model, the result does not change with the expansion of the domain. Therefore, the 14*1 domain was chosen as the most suitable one. Fig. 8 shows the dependence of the pressure drop result on the distance (number of cells) from the inlet surface (Fig. 3). As a reference, the result on cells 9 and 10 is used in this graph. -9.0 Fig. 8. Dependence of the result of the pressure drop on the domain length As mentioned above, the pressure difference along the domain varies slightly due to the inlet and outlet boundary conditions. It stabilizes satisfactorily on cells between 7 and 10 as shown in Fig. 8. Therefore, the pressure difference on cells 9 and 10 was chosen for validation. Fig. 9. Dependences of measurement result on water temperature in laboratory measurement Pressure drop is also a function of water temperature, so a laboratory analysis was performed to evaluate the dependence of pressure drop on temperature. Measurements required for this analysis were conducted with the measuring method as described in chapter 2. Fig. 9 shows the dependence of the result of pressure drop on the water temperature of the laboratory measurement. The water temperature at 15.2 ° C is used as a reference in this graph. It can Numerical Analysis with Experimental Validation of Single-Phase Fluid Flow in a Dimple Pattern Heat Exchanger Channel 549 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 544-553 Fig. 10. a) Measurement and CFD results with associated power-law function, and b) relative difference between laboratory measurements and different numerical solutions of k- e and k- m turbulence models be seen that the dependence in this domain is expected to be linear and also that the deviation increases nonlinearly as water flow decreases. The uncertainty of the measured pressure drop due to temperature deviation during the measurement is significantly lower than 0.5 %. Black colored curve in Fig. 10 on the left chart shows the results of laboratory measurements on a sample of 31*19 cells. As CFD results, also laboratory measurements are presented for two consecutive cells calculated by dividing the measured pressure drop result with 31 and multiplied by 2. There are also results from CFD calculations obtained by using k- s (blue) and k- m (orange) turbulence model as presented above. All curves have associated power-law functions. In Fig. 10b chart shows the relative difference between numerical calculations and laboratory results, where results on abscissa contain the results of said laboratory measurements. Fig. 11 shows the relative difference between the calculated values (from the regression equations) and the laboratory measurements or CFD calculations. Velocity fm/sj Fig. 11. Relative difference between calculated values (from regression equations) and laboratory measurements or CFD calculations 550 The results of laboratory measurements on samples with lengths of 24 and 16 cells are not presented in this paper. However, the results of pressure drop for a sample with 24 cells are on average about 1 % higher and for a sample with 16 cells about 2 % higher than for the longest sample. The difference increases as the flow decreases. We selected the results for the longest sample as the most appropriate one. Due to the deliberate design of mounting the sample into the test bench, it can be assumed that the pressure conditions at each individual inlet (or outlet) surface are the same. Therefore, the water velocities are the same in all parts of the sample. The influence of the wall on the pressure drop is negligible in relation to the relatively large width of the test sample. Also, the influence of possible temperature inhomogeneity of water in the sample is neglected. However, the influence of the sample length is not negligible. We could assume that the fluid flow in a large vessel is completely laminar. It is laminar also just before it enters the sample. Even when it enters the cells, the flow does not immediately turn completely to fully developed turbulent flow, but it does require some cells. Although the exact length of this inlet region remains unknown, we can say that it is the same for all samples. It can also be said that the pressure drop in this region is larger than with the fully developed turbulent flow inside the sample. Thus, the influence of the inlet region on the pressure drop per cell in the short sample is greater than in the longer one. In this way, we can explain why the pressure drop per cell increases with the shortening of the samples. Močnik, U. - Blagojevič, B. - Muhič, S. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 544-553 The graph in Fig. 10a shows relatively nice curves, which also match the corresponding power-law functions, indicating that the laboratory measurements and numerical calculations were performed correctly. How well the power-law functions fit with the measured data or CFD calculations is shown in the graph in Fig. 11. We can conclude that all values are within the interval ±1.6 %. From the graph in Fig. 10a, we can see that the k- s on average predicts an 8.4 % lower pressure drop compared to the measured values. This difference is approximately the same in the entire area considered, which can be clearly seen from the graphs in Fig. 10b and in Fig. 11. The results of numerical calculations with the k- m turbulence model predict an average difference of 9.7 % in pressure drop compared to the results of laboratory measurements. This difference is not constant, but varies with speed, as can be seen from the graph in Fig. 10b. The agreement with the power curve is in a similar ratio to the measurements but with the opposite shape, which can be seen in the graph in Fig. 11. Fig. 12 shows the volume of domain in which the velocity of water is positive in the >>-direction (that is in the direction from inlet to outlet). The volume in the Fig. 12a is prediction of the calculation with k- s and the volume in the Fig. 12b is prediction of the calculation with the use of k- m turbulence model. Zone where the velocity in the >>-direction is positive was named "live zone", zone where the velocity in the >>-direction is negative was named "dead zone". The shape of "live zone" can be seen on Fig. 12. Fig. 12. The volume of the domain in which the velocity of water in they-dlrectlon Is positive: a) k- s, and b) k- m From the Fig. 10b we can see that the use of the k- m turbulence model predicts on average 19.8 % greater pressure drop than the prediction from model where k - s was used. This may also be because the k- m turbulence model predicts a significantly different "dead zone" behind the brazing point as could be seen in Fig. 12. The "dead zone" prediction of k - s turbulence model is not only different in shape, it is also significantly smaller when compared with prediction of k- a> turbulence model. The "dead zone" is 22.3 % of whole volume for the model when k - s turbulence model was used and 31.6 % of the whole volume when the k- a> turbulence model was used. The size of the "dead zone" practically does not change with velocity in the analysed velocity range. In our case, use of the k- s turbulence model (in the case of steady-state simulation) proved to be a robust and reliable tool for predicting the pressure drop in the analysed geometry under consideration. We had no major convergence problems, no matter what mesh density we used for calculations. This is also evidenced by the study in which we determined the dependence of the result from the computational mesh density, where the results are smoothly limited to the true value (Fig. 6). Even along with the domain, the result does not change significantly (Fig. 8), and as a result, the domain could be significantly shorter, which in turn means a shorter computation time or a higher mesh density and thus a better result. The results at different velocities behave similarly to the measurement results (Figs. 10 and 11). It yielded useful results comparable to those of laboratory measurements. With the use of the k- a> turbulence model we did not obtain solutions for long samples, because they behave unexpectedly with increasing mesh density, the results along the domain are unusual, and the results behave differently at different velocities than the measurement results (Figs. 10 and 11). It seems that with use of the k- a> turbulence model we have obtained fairly unreliable results. However, it must be understood that the flow (at given velocities in the channel) is turbulent. This can be assured with certainty, since the relation between the pressure drop and the flow rate is not linear, which would mean laminar flow. With the increase of mesh density, the solution approaches the true value due to the higher resolution in solving the equations, but there are more and more problems with convergence, which indicates a strongly turbulent flow, which in the densest mesh becomes too intense and the time independent solution finally fails. It can be concluded that steady-state simulations with k- a> turbulence model are not suitable for predicting the pressure drop in analysed heat exchanger. However, this means that in this case the k- s turbulence model is a safe choice to predict Numerical Analysis with Experimental Validation of Single-Phase Fluid Flow in a Dimple Pattern Heat Exchanger Channel 551 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 544-553 pressure drop in channels formed by dimple pattern heat plate. 4 CONCLUSIONS In this article, we have presented a comprehensive analysis of the flow state in a plate heat exchanger with a dimple pattern heat plate structure. The research was conducted in two interrelated fields: laboratory and numerical. We carried out several laboratory measurements of the channel that makes up the XB12L. The methodology of work was developed and the result was a completely new methodology for conducting laboratory measurements with the so-called large vessels. This measurement method has yielded the results we can trust. A similar method of performing measurements on a heat exchanger channel cannot be found in the literature. We have found that the fluid flow in the channel under consideration is turbulent at all velocities considered. In the numerical part of the research, we have determined the optimal mesh density, which lias a very favourable relationship between the accuracy of the results and the consumption of computer resources. We have determined the optimal numerical domain size for use of different turbulence models. We have found that in steady-state simulations use of k-e realizable turbulence model with Enhanced wall treatment is more appropriate than use of the k-m SST turbulence model. Future research will cover the implementation of time-dependent simulations and introducing energy equations for heat transfer in numerical work, and by completing laboratory test bench in a way that allows the measurement of heat transfer. 5 NOMENCLATURES f External body forces, [N] g Gravitational acceleration, [m/s2] p Pressure, [Pa] u Velocity, [m/s] s„, Mass added to the continuous phase, [kg] uT Velocity at a distance Av, [m/s] Av Distance from the wall, [m] v+ Dimensionless distance from the wall, [-] u Dynamic viscosity, [kg/(ms)] p Density, [kg/m3] t0} Wall shear stress, [Pa] t Stress tensor [-] 6 REFERENCES [1] Fernandes, C.S, Dias, R.P, Maia, J.M. (2008). New plates for different types of plate heat exchangers. Recent Patents on Mechanical Engineering, vol. 1, no. 3, p. 198-205, D0l:10.217 4/2212797610801030198 [2] Kandlikar, S.G, Grande, W.J. (2003). Evolution of MicroChannel flow passages - Thermohydraulic performance and fabrication technology. Heat Transfer Engineering, vol. 24, no. 1, p. 3-17, 2003, D0l:10.1080/01457630304040 [3] Focke, W.W, Knibbe, P.G. (1986). Flow visualization in parallelplate ducts with corrugated walls. Journal of Fluid Mechanics, vol 165, p. 73 77, D0l:10.1017/S0022112086003002 [4] Metwally, H.M, Manglik, R.M. (2004). Enhanced heat transfer due to curvature-induced lateral vortices in laminar flows in sinusoidal corrugated-plate channels. International Journal of Heat and Mass Transfer, vol. 47, no. 10-11, p. 2283-2292, 2004, D0l:10.1016/j.ijheatmasstransfer.2003.11.019 [5] Kanaris, A.G, Mouza, A.A., Paras, S.V. (2005). Flow and heat transfer in narrow channels with corrugated walls a CFD code application. Chemical Engineering Research and Design, vol. 83, no 5, p. 460-468, D0l:10.1205/cherd.04162 [6] Han, X.H., Cui, L.Q., Chen, S.J., Chen, G.M., Wang Q. (2010). A numerical and experimental study of chevron, corrugated-plate heat exchangers. International Communications in Heat and Mass Transfer, vol. 37, no. 8, p. 1008-1014, D0l:10.1016/j. ¡cheat masstransfer.2010.06.026 [7] Gherasim, I, Galanis, N, Nguyen, C.T. (2011). Heat transfer and fluid flow in a plate heat exchanger. Part II: Assessment of laminar and two-equation turbulent models. International Journal of Thermal Sciences, vol. 50, no. 8, p. 1499-1511, D01:10.1016/j.ijthermalsci.2011.03.017 [8] Gullapalli, V.S, Sunden, B. (2014). CFD simulation of heat transfer and pressure drop in compact brazed plate heat exchangers. Heat Transfer Engineering, vol. 35, no. 4, p. 358366, D0l:10.1080/01457632.2013.828557 [9] Tiwari, A.K, Ghosh, P, Sarkar, J, Dahiya, H, Parekh, J. (2014). Numerical investigation of heat transfer and fluid flow in plate heat exchanger using nanofluids. International Journal of Thermal Sciences, vol. 85, p. 93-103, D0l:10.1016/j. ijthermalsci.2014.06.015 [10] Lee, J, Lee, K.-S. (2014). Flow characteristics and thermal performance in chevron type plate heat exchangers. International Journal of Heat and Mass Transfer, vol. 78, p. 699-706, D0l:10.1016/j.ijheatmasstransfer.2014.07.033 [11] Sarraf, K, Launay, S, Tadrist, L. (2015). Complex 3D-flow analysis and corrugation angle effect in plate heat exchangers. International Journal of Thermal Sciences, vol. 94, p. 126138, D01:10.1016/j.ijthermalsci.2015.03.002 [12] Piper, M, Zibart, A„ Tran, J.M., Kenig E.Y. (2016). Numerical investigation of turbulent forced convection heat transfer in pillow plates. International Journal of Heat and Mass Transfer, vol. 94, p. 516-527, D0l:10.1016/j. ijheatmasstransfer.2015.11.014 [13] Yogesh, S.S., Selvaraj, A.S, Ravi, D.K, Rajagopal, T.K.R. (2018). Heat transfer and pressure drop characteristics of inclined elliptical fin tube heat exchanger of varying ellipticity ratio using CFD code. International Journal of 552 Močnik, U. - Blagojevič, B. - Muhič, S. Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, 544-553 Heat and Mass Transfer, vol. 119, p. 26-39, DOI:10.1016/j. ijheatmasstransfer.2017.11.094. [14] Al zahrani, S., Islam, M.S., Saha, S.C. (2019). A thermo-hydraulic characteristics investigation in corrugated plate heat exchanger. Energy Procedia, vol. 160, p. 597-605, DOI:10.1016/j.egypro.2019.02.211. [15] Muley, A., Manglik, R.M. (1999). Experimental study of turbulent flow heat transfer and pressure drop in a plate heat exchanger with chevron plates. Journal of Heat Transfer, vol. 121, no. 1, p. 110-117, DOI:10.1115/1.2825923. [16] ANSYS Inc. (2018). ANSYS Fluent Theory Guide, Release 19. ANSYS Inc., Canonsburg. [17] Peric, M., Ferguson, S. (2005). The advantage of polyhedral meshes, CD-adapco. Melville. Numerical Analysis with Experimental Validation of Single-Phase Fluid Flow in a Dimple Pattern Heat Exchanger Channel 553 Strojniški vestnik- Journal of Mechanical Engineering 66(2020)9 Vsebina Vsebina Strojniški vestnik - Journal of Mechanical Engineering letnik 66, (2020), številka 9 Ljubljana, september 2020 ISSN 0039-2480 Izhaja mesečno Razširjeni povzetki (extended abstracts) Andraž Lipolt, Brane Širok, Marko Hočevar, Lovrenc Novak: Konvektivno sušenje plasti odpadnega blata s tokom skozi plast SI 61 Anna Trelka, Wojciech Zórawski, Anna Góral: Modificiranje mikrostrukture in lastnosti hladno napršenih prevlek z različnimi velikostmi zrn kompozitnega prahu Cr3C2-25(Ni20Cr) SI 62 Frank Goldschmidtboeing, Uwe Pelz, Karen Claire-Zimmet, Michael Wolf, Ralf Goerlach, Peter Woias: Gibanje ščetin, sile in vertikalni pomik pri novi zasnovi električne zobne ščetke SI 63 Xiangyang Xu, Xupeng Fan, Peitang Wei, Baojun Yang: Raziskava mazanja v področju ubiranja pri harmoničnih zobniških prenosnikih SI 64 Andres Lopez-Lopez, Jose Billerman Robles-Ocampo, Perla Yazmin Sevilla-Camacho, Orlando Lastres-Danguillecourt, Jesús Muniz, Bianca Yadira Perez-Sariflana, Sergio de la Cruz: Dinamična nestabilnost lopatic vetrnih turbin zaradi velikih odklonov: eksperimentalna validacija SI 65 Yan Li, Lianhui Li: Izboljšava optimizacije izbire programa produktno-servisnega sistema z okviijem na osnovi digitalnega dvojčka SI 66 Urban Močnik, Bogdan Blagojevič, Simon Muliič: Numerična analiza z eksperimentalno validacijo enofaznega toka tekočine v kanalu ploščnega prenosnika toplote z jamičasto strukturo SI 67 Strojniški vestnik-Journal of Mechanlcal Englneering 66(2020)9, SI 61 © 2020 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2020-04-12 Prejeto popravljeno: 2020-07-21 Odobreno za objavo: 2020-08-12 Konvektivno sušenje plasti odpadnega blata s tokom skozi plast Andraž Lipolt1 - Brane Širok2 - Marko Hočevar2 - Lovrenc Novak2 * 1 Petrol d.d., Slovenija 2 Univerza v Ljubljani, Fakulteta za strojništvo, Slovenija Odpadno blato s čistilnih naprav vsebuje mnoge kontaminante, a tudi mnoge koristne elemente (predvsem fosfor), zato se metodam obdelave blata posveča vedno več pozornosti. Tradicionalno odlaganje na deponijah se nadomešča z naprednejšimi metodami obdelave, ki vključujejo sežig blata ter izločitev nutrientov. Zaradi velike vsebnosti vode je blato pred sežigom potrebno posušiti, za kar se najpogosteje uporablja metoda konvektivnega sušenja. V prispevku je obravnavano konvektivno sušenje blata z lastnostmi, ki so značilne za tračne sušilnike. Medtem, ko večina objavljenih raziskav na področju sušenja blata obravnava majhne vzorce z maso nekaj gramov, so bili v predstavljeni študiji izvedeni eksperimenti na laboratorijski sušilni napravi z vzorci mase več kilogramov, ki so bili oblikovani v 4 cm debele plasti velikosti 0,5 m x 0,5 m. Izbrana velikost vzorca je omogočala večjo izraznost vplivov na sušenje, ki izhajajo iz časovne spremenljivosti porozne plasti zaradi krčenja in pokanja. V nasprotju z večino objavljenih raziskav je bil tok sušilnega zraka usmeijen skozi plast. V takih razmerah poroznost plasti tako na lokalni kot integralni skali znatno vpliva na tok zraka in s tem na hitrost sušenja blata. Sušilni procesi so potekali pri temperaturi zraka 65 °C in 80 °C in pri hitrostih zraka 0,53 m/s in 0,83 m/s. Ustrezna relativna vlažnost sušilnega zraka je bila zagotovljena s pomočjo konstantnega nadomeščanja relativno vlažnega zraka iz sušilne komore s svežim zrakom. Merilni sistem je omogočal stalen nadzor parametrov sušilnega zraka nad in pod plastjo, opazovanje površine plasti z IR kamero ter tehtanje plasti blata. Tehtana sila je poleg mase blata vsebovala tudi silo upora zaradi prehoda zraka skozi plast. Meritev dejanske mase blata je bila izvedena z občasnim ustavljanjem toka zraka skozi plast in s tem izločitvijo sile upora. Analiza rezultatov je bila osredotočena na ovrednotenje kinetike sušenja blata in v razvoj modelov za opis sušenja blata na integralnem nivoju. Zmanjševanje sile upora v odvisnosti od časa je bilo modelirano s preprosto funkcijo, kije dala odlično ujemanje modela z izbranimi meritvami. Sušilna kinetika plasti blata je bila določena s prileganjem običajnih sušilnih modelov meijenim podatkom. Na osnovi izbranih statističnih kriterijev (vrednosti RMSE, R2, z2) ni bilo mogoče definirati enotnega modela sušenja, primernega za vse obratovalne pogoje (različne temperature in hitrosti sušilnega zraka), zato sta bila kot statistično najprimernejša določena dva modela in sicer model Wang Singh in modificirani model Nadhari. Čas sušenja na kilogram blata je bil modeliran kot funkcija temperature sušilnega zraka, hitrosti sušilnega zraka in začetne vsebnosti suhe snovi v blatu. Koeficient detenninacije (R2) za razviti model znaša 0,944. Znotraj intervala preučevanih sušilnih pogojev so bili dobljeni sušilni časi med 43 in 76 minut na kilogram blata. Ker so razmere v predstavljeni študiji zelo podobne razmeram v industrijskih tračnih sušilnikih so predstavljeni rezultati primerni za skaliranje in uporabo na dejanskih sušilnih napravah. Rezultati so prikazali pomembnost spreminjanja lastnosti plasti blata, predvsem poroznosti, ki se odraža v padcu tlaka skozi plast in posledično v sili upora. Z analizo trenutnih podatkov ni bilo mogoče določiti modelov, ki bi povezali poroznost plasti s parametri sušilnega zraka in sušilno kinetiko. Izvedeni bodo nadaljnji eksperimenti z nadgrajenim merilnim sistemom, s katerimi bo mogoče analizirati spremembo strukture plasti zaradi krčenja in pokanja kot funkcijo časa ter zasnovati odnose med temi pojavi in sušilnim procesom. Ključne besede: odpadno blato s čistilnih naprav, sušenje v tanki plasti, porozna plast, sušilna kinetika, sila upora, modeliranje *Naslov avtorja za dopisovanje: Univerza v Ljubljani, Fakulteta za strojništvo, Aškerčeva 6,1000 Ljubljana, Slovenija, lovro.novak@fs.unMj.si SI 61 Strojniški vestnik-Journal of Mechanlcal Englneering 66(2020)9, SI 62 © 2020 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2020-03-24 Prejeto popravljeno: 2020-05-19 Odobreno za objavo: 2020-07-15 Modificiranje mikrostrukture in lastnosti hladno napršenih prevlek z različnimi velikostmi zrn kompozitnega prahu Cr3C2-25(Ni20Cr) Anna Trelka1 * - Wojciech Žorawski2 - Anna Goral1 1 Poljska akademija znanosti, Inštitut za metalurgijo in materiale, Poljska 2 Tehniška univerza Kielce, Raziskovalno središče za lasersko obdelavo, Poljska Pričujoča študija prinaša nova spoznanja o hladno napršenih kompozitnih prevlekah, še posebej o vplivu velikosti zrn prahu in odmika na mikrostrukturo in lastnosti prevlek iz cenneta Cr3C2-25(Ni20Cr) na substratu iz zlitine Al 7075. Namen študije je opredelitev vpliva treh različnih velikosti zrn prahu (DL: 9,5-40 |im. DO: 9,5-55,3 |im in DT: 40-55,3 |im) in treh različnih vrednosti odmika (10 mm, 30 mm in 50 mm) na mikro strukturo (sestava faz, topografija površine, volumski delež in porazdelitev faze Cr3C2 v osnovi Ni20Cr ter poroznost) in mehanske lastnosti (trdota in obstojnost proti abrazivni obrabi) prevlek iz Cr3C2-25(Ni20Cr). Proces nabrizgavanja v hladnem je bil opravljen s sistemom Impact Innovations 5/8, nameščenim na robotski roki Fanuc M-20iA. Vsebnost keramične faze v prahovih je bila določena z rentgensko difrakcijo in s slikovno analizo. Analiza debeline prevlek je bila opravljena s programsko opremo ImageJ na posnetkih, ustvaijenih z optičnim mikroskopom. Mikrostruktura hladno napršenih prevlek Cr3C2-25(Ni20Cr) je bila določena z vrstično elektronsko mikroskopijo (FEI/Philips XL30). Analiza topografije površine je bila opravljena z mikroskopom Keyence VHX-7000. Trdota prevleke po prerezu je bila določena s sistemom za meijenje trdote po Vickersu z majhno silo (HV0.3) CSM Instruments. Preizkusi abrazivne obrabe so bili opravljeni z napravo ITEE T-07 z gumijastim kolesom in suhim abrazivom. Študije so pokazale, da imajo vse hladno napršene prevleke Cr3C2-25(Ni20Cr) kompaktno mikrostrukturo. Najbolj porozne so bile prevleke iz prahu z največjimi delci. Prevleke, izdelane iz prahov z delci velikosti nad 40 (nn (DT), so imele največji volumski delež keramične faze Cr3C2, najmanjši volumski delež pa je bil ugotovljen pri delcih, manjših od 40 |im (DL). Odmik je imel majhen vpliv na vsebnost keramične faze v prevlekah. Debelina hladno napršenih prevlek Cr3C2-25(Ni20Cr) se zmanjšuje z rastjo velikosti zrn prahu in s povečevanjem odmika. Največja trdota je bila ugotovljena pri prevlekah iz prahu z izhodiščno velikostjo zrn (DO: 9,5 (nn do 55,3 |im). Najnižja trdota je bila ugotovljena pri prevlekah DT, izdelanih iz prahu z zrni velikosti 40 (nn do 55,3 (nn. Hladno napršene prevleke, izdelane iz prahov z večjimi delci (DO in DT), so bile v primeijavi s prevlekami, izdelanimi iz najfinejših prahov, obstojnejše proti obrabi, povzročeni z delci abraziva. Najboljše lastnosti so imele prevleke Cr3C2-25(Ni20Cr), izdelane iz prahu DO pri odmiku 30 mm: veliko debelino, visoko trdoto, visoko obstojnost proti abrazivni obrabi, razmeroma velik volumski delež faze Cr3C2 in majhno poroznost. Raziskovalci se trenutno soočajo z izzivom razvoja novih hladno napršenih kompozitnih prevlek z lastnostmi, primernimi za avtomobilsko in letalsko industrijo. Ta študija premošča vrzeli v znanju s karakterizacijo vpliva odmika in velikosti zrn prahu na mikrostrukturo in mehanske lastnosti prevlek iz Cr3C2-25(Ni20Cr). Ključne besede: kompozitna prevleka Cr3C2-25(Ni20Cr), nabrizgavanje v hladnem, mikrostruktura, velikost zrn prahu, odmik, mehanske lastnosti SI 62 *Naslov avtorja za dopisovanje: Poljska akademija znanosti. Inštitut za metalurgijo In materiale, ul. Reymonta 25, 30-059 Krakow, Poljska, a.trelka@lmim.pl Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, SI 63 © 2020 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2020-03-11 Prejeto popravljeno: 2020-06-18 Odobreno za objavo: 2020-07-14 Gibanje ščetin, sile in vertikalni pomik pri novi zasnovi električne zobne ščetke Frank Goldschmidtboeing1*- Uwe Pelz1 - Karen Claire-Zimmet2 - Michael Wolf2 - Ralf Goerlach2 - Peter Woias1 1 Univerza v Freiburgu, Oddelek za mikrosisteme, Laboratorij za konstruiranje mikrosistemov, Nemčija 2 Procter & Gamble Service GmbH, Nemčija Članek preučuje dinamiko nove zobne ščetke z rotacijsko-oscilacijskim pogonom, s poudarkom na zunajravninskem gibanju glave ščetke. Gibanje ščetin, sile in vertikalni pomiki pri novi zasnovi električne zobne ščetke z linearnim pogonom so okarakterizirani s kombinacijo teoretičnih in eksperimentalnih tehnik. Razvita je bila teorija na osnovi posameznih vlaken in opravljena je bila eksperimentalna določitev sil in gibanj pri ščetkanju s pomočjo senzorja sile in visokohitrostne kamere. Eksperimenti so pokazali, da električna ščetka z linearnim rotacijsko-oscilacijskim pogonom ustvari najvišjo silo v točkah obračanja, ki jih je mogoče prilagajati z razmerji med različno usmerjenimi vlakni. Nova zasnova zagotavlja vertikalni pomik z amplitudo do 250 ^m. Signal sile pri ščetki z linearnim pogonom sestoji samo iz pogonske frekvence in iz njenih harmonikov. Pri modeliranju dinamike glave ščetke je bila zanemarjena njena vztrajnost. Vpliv vztrajnosti bo zato treba vključiti v novejše različice teorije posameznih vlaken. V prihodnjih raziskavah bodo potrebni tudi dodatni eksperimenti za primerjavo zasnov ščetk z različno razporeditvijo čopov oz. z različnimi koti ščetin. Učinkovitost čiščenja zob običajno preučujejo klinične študije, ki so objavljene v specializiranih zobozdravstvenih revijah, zato praktično ni virov o mehanskih osnovah ali o fiziki gibanja vlaken zobnih ščetk. Razlog za pomanjkanje javno objavljenih in recenziranih raziskav je verjetno v kompleksnih interakcijah ščetin med seboj in s kompleksno površino zoba. Naš pristop uporablja model posameznih vlaken in razširitev teorije posameznih vlaken za vključitev vertikalnega gibanja glave ščetke, kakor tudi določitev odmika in sile interakcij s površino zoba. Pri zobnih ščetkah kljub vsakodnevni uporabi in večstoletni zgodovini še obstajajo odprta vprašanja, ki se jih je mogoče lotiti teoretično in eksperimentalno. Globlje razumevanje dinamike vlaken bo osnova za prihodnje izboljšave zasnove zobnih ščetk. Ključne besede: teorija posameznih vlaken, električna zobna ščetka, gibanje ščetin, sile ščetin, poševna vlakna zobne ščetke *Naslov avtorja za dopisovanje: Univerza v Freiburgu, Oddelek za mikrosisteme, Laboratorij za konstruiranje mikrosistemov, Nemčija, fgoldsch@imtek.de SI 63 Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, SI 64 © 2020 Strojniški vestnik. Vse pravice pridržane. Raziskava mazanja v področju ubiranja pri harmoničnih zobniških prenosnikih Xiangyang Xu1 * - Xupeng Fan1 - Peitang Wei2 - Baojun Yang3 1 Univerza Chongqing, Šola za mehatroniko in avtomobilsko tehniko, Kitajska 2 Univerza Chongqing, Šola za strojništvo, Kitajska 3 Chongqing Huashu Robot Co., Ltd., Chongqing 400714, Kitajska Harmonični zobniški prenosniki so nova vrsta prenosnikov, ki se razlikujejo od klasičnih zobniških gonil. Moč prenašajo predvsem z deformacijo upogljivega zobnika, njihova prednost pa je v majhnih izmerah, visoki točnosti, visokem prestavnem razmeiju, stabilnosti prenosa in visokem izkoristku. Za harmonične zobniške prenosnike pa je značilna tudi nelinearna togost, histereza in nelinearno trenje. To so dejavniki, ki omejujejo njihovo točnost in izkoristek. V literaturi je na voljo le malo kvantitativnih analiz vpliva delovnih pogojev na mazalne lastnosti in togost oljnega filma v področjih ubiranja v odvisnosti od kontaktne geometrije zob. Zato obstaja potreba po preučitvi zmogljivosti mazanja sistemov harmoničnih zobniških prenosnikov za izboljšanje njihove učinkovitosti in točnosti. Pričujoči članek za analizo mazanja harmoničnih zobnikov in postavitev osnov za preučevanje odpovedi zob zobnikov ter njihovih dinamičnih lastnosti analizira krivinski radij, obremenitev zob in hitrosti v točki ubiranja zob v odvisnosti od kontaktne geometrije zob harmoničnih zobnikov. Postavljen je model elastohidrodinamičnega mazanja (EHL) z linijskim kontaktom omejene dolžine. Enačbe za elastično deformacijo, viskoznost maziva in gostoto so v članku postavljene kot funkcije tlaka ter kombinirane z Reynoldsovo enačbo in razrešene po kompozitni iterativni metodi. Rezultati po velikem številu iteracij konvergirajo k vrednosti, ki izpolnjuje zahtevano natančnost. Z numeričnimi računskimi metodami sta bili preučeni debelina oljnega filma in porazdelitev tlaka v področju mazanja, kakor tudi vpliv vrtilne hitrosti in temperature na razmeije kontaktnih obremenitev in razmeije debelin filma v področju ubiranja, poleg tega pa tudi sprememba debeline oljnega filma v različnih delovnih pogojih. Rezultati so pokazali, daje tlak v smeri ubiranja na začetku majhen in nato doseže največjo vrednost v središču, debelina filma pa je največja v predelu vstopa in enakomerno porazdeljena v sredinskem kontaktnem območju. Z naraščanjem hitrosti se znižuje razmeije kontaktne obremenitve zoba, povečuje se razmeije debelin oljnega filma, togost oljnega filma se občutno zniža in izboljša se učinek mazanja; medtem ko je učinek temperature prav nasproten. Z ustreznim povišanjem vrtilne hitrosti in znižanjem temperature olja je torej mogoče učinkovito izboljšati mazalne lastnosti sistema. Namen pričujoče študije je numerična analiza lastnosti mazanja v predelu ubiranja harmoničnega zobniškega prenosnika ter vpliva delovnih pogojev nanje. Rezultati študije predstavljajo teoretično osnovo in referenčno orodje za nadaljnje raziskave dinamičnih lastnosti, načinov odpovedi in zanesljivosti sistemov harmoničnih prenosnikov v praktičnih delovnih pogojih. Rezultati bodo tako uporabni tudi za prihodnje študije mehanizmov odpovedi harmoničnih zobniških prenosnikov zaradi obrabe in za nelinearno dinamično analizo. Ključne besede: harmonični zobniški prenosnik, elastohidrodinamično mazanje, področje ubiranja, delovni pogoji, numerična metoda, lastnosti mazanja Prejeto v recenzijo: 2020-03-11 Prejeto popravljeno: 2020-06-18 Odobreno za objavo: 2020-07-14 SI 64 *Naslov avtorja za dopisovanje: Univerza Chongqing, Šola za mehatroniko In avtomobilsko tehniko. Kitajska, xxiangyang@hotmall.com Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, SI 65 © 2020 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2020-03-11 Prejeto popravljeno: 2020-06-18 Odobreno za objavo: 2020-07-14 Dinamična nestabilnost lopatic vetrnih turbin zaradi velikih odklonov: eksperimentalna validacija Andrés Lopez-Lopez1 - José Billerman Robles-Ocampo2 * - Perla Yazmin Sevilla-Camacho2 -Orlando Lastres-Danguillecourt1 - Jesús Muniz3 - Bianca Yadira Perez-Sariñana2 - Sergio de la Cruz2 1 Znanstveno-humanistična univerza v Chiapasu, Mehika 2 Politehnika v Chiapasu, Akademska skupina za energetiko in trajnostni razvoj, Mehika 3 Mehiška avtonomna nacionalna univerza, Inštitut za obnovljivo energijo, Mehika Vetrne turbine s horizontalno osjo (HAWT) so konstruirane z velikimi rotoiji, ki jim omogočajo, da izkoristijo kar se da velik del kinetične energije vetra. Vsak rotorje sestavljen iz pesta in lopatic. Lopatice so dolge, vitke, lahke, močno upogljive in prenašajo velike vetrne obremenitve. Med eksploatacijo so izpostavljene velikim odklonom in visokoamplitudnim vibracijam. Prav kombinacija teh dveh pojavov lahko povzroči nestabilnost lopatic, ki se nato prenaša naprej na gred rotoija in na ostale komponente vetrne turbine, to pa lahko privede do obsežnih poškodb sistema. Zato je pomembno podrobnejše razumevanje nestabilnosti kot odziva lopatic vetrne turbine na vibracije v prisotnosti velikih odklonov. Opravljenih je bilo več raziskav za razjasnitev nelinearnega dinamičnega vedenja lopatic. Nelinearno dinamično vedenje se lahko analizira npr. s kombiniranjem 3D modelov s končnimi elementi in računalniške dinamike fluidov (CFD). Te analitične metode so sicer učinkovite, so pa tudi zelo drage in še vedno ne omogočajo jasnega opisa dinamične nestabilnosti rotoijev pri različnih hitrostih vetra oz. vpliva večjih odklonov na nestabilnost. V pričujočem delu je bil za odpravo teh težav uporabljen poenostavljen aeroelastični model. Model omogoča identifikacijo, analizo in napovedovanje vpliva velikih odklonov na nelinearne vibracije lopatic vetrnih turbin. Uporabljen je bil model konzolnega nosilca z nelinearnimi členi, ki opisujejo elastično krivuljo. Oblikovanje bil dinamični model z upoštevanjem velikega geometrijskega odklona in uporabljen je bil pristop po Galerkinu. Privzete so bile enakovredne elastične in vztrajnostne lastnosti, vpliv sile vetra paje bil formuliran po načelu diska aktuatoija, ki omogoča spreminjanje hitrosti vetra. Aeroelastični model je bil nato eksperimentalno validiran. Eksperimentalni sistem je vključeval vetrovnik, pospeškomer, sistem za zajem podatkov in aluminijaste lopatice. Eksperimenti so bili opravljeni z lastnimi in z vsiljenimi vibracijami. Vsiljene vibracije so bile generirane v vetrovniku, ki ustvaija hitrosti vetra od 10 m/s do 30 m/s. Zajeti signali pospeškov so bili integrirani v hitrosti in odmike, ki so predstavljeni v obliki faznih diagramov. Ti diagrami omogočajo opredelitev stabilne, nestabilne oz. kaotične dinamike lopatic. Numerični in eksperimentalni rezultati analize dinamike lopatic razkrivajo, da je stabilnost odvisna od amplitude vzbujanja in od togosti. To pomeni, da je vedenje lopatic pri nizkih vetrnih obremenitvah stabilno, odkloni so majhni in sistem je linearen. Z naraščanjem vetrne obremenitve se nato pojavi nestabilnost, pri čemer so najbolj prizadete lopatice z nizko togostjo. Rezultati preskusov z lastnimi vibracijami so pokazali stabilen odziv. Večji odkloni lopatic vnašajo harmonične motnje oz. nelinearnosti in fazna ravnina izkazuje motnje v kinetični in potencialni energiji vibracij. Glavni prispevek pričujočega dela je v razvoju modela za določanje nelinearnih parametrov vibracij lopatic vetrnih turbin in identifikacijo njihove nestabilnosti pod vplivom različnih hitrosti vetra. Za model sta značilni računska ekonomičnost in učinkovitost pri identifikaciji geometrijskih nelinearnosti. Model bo v prihodnjih raziskavah mogoče razširiti z upoštevanjem sprememb preseka, ortotropnosti kompozitnih materialov in značilnih sekcijskih ojačitev lopatic velikih komercialnih vetrnih turbin, s tem pa bo dopolnjena ta študija nestabilnosti in nelinearnih vibracij lopatic. Ključne besede: nestabilnosti, vetrovnik, nelinearne vibracije, nelinearen model, velik odklon, lopatica vetrne turbine, fazna ravnina *Naslov avtorja za dopisovanje: Znanstveno-humanistična univerza v Chiapasu Carretera Tuxtla Gutierrez. - Portillo Zaragoza Km 21+500, Col. Las Brisas: Suchiapa, Chiapas. CP29150, Mehika, jrobles@upchiapas.edu.mx SI 65 Strojniški vestnik-Journal of Mechanical Engineering 66(2020)9, SI 66 © 2020 Journal of Mechanical Engineering. All rights reserved. Prejeto v recenzijo: 2020-02-17 Prejeto popravljeno: 2020-07-13 Odobreno za objavo: 2020-07-27 Izboljšava optimizacije izbire programa produktno-servisnega sistema z okvirjem na osnovi digitalnega dvojčka Yan Li1 - Lianhui Li2 * 1 Xinxiang univerza, Šola za informatiko in inženiring, Kitajska 2 Severna Minzu univerza, Šola za strojništvo in elektrotehniko, Kitajska Produktno-storitveni sistemi (PSS) so bili razviti za proizvodna podjetja kot orodje za zagotavljanje individualiziranih izdelkov in storitev za kupce. Namen pričujoče raziskave je izboljšanje optimizacije izbire programa PSS v fazi načrtovanja PSS. Glede na dinamične lastnosti večdimenzij skih vplivnih dejavnikov in njihove medsebojne odvisnosti je podan predlog okviija za izboljšavo optimizacije izbire programa PSS na osnovi digitalnega dvojčka. Predlagani okvirje razdeljen v tri sloje: digitalni dvojček, informacijski sloj in pristopni sloj. Podane so logične odvisnosti med sloji in zasnovan je kvantitativni mehanizem za optimizacijo izbire programa PSS. Za določitev vrednosti atributov v programu PSS sta vključena mehko število in grobi mejni interval. Za oceno programa PSS je bila uporabljena prilagojena metoda TOPSIS, v kateri je bila evklidska razdalja zamenjana z relacijsko vektorsko razdaljo. Preučen je primer optimizacije izbire programa PSS za čiščenje zraka s predlaganim okviijem na osnovi digitalnega dvojčka. Dokazana je uspešnost mehanizma za optimizacijo izbire programa PSS, ki jo je torej mogoče izboljšati s predstavljenim okviijem. Zaradi določenih omejitev pri optimizaciji bodo v prihodnje potrebne še dodatne raziskave. Z modeliranjem ocenjevalnih podatkov za izbiro programa PSS na osnovi digitalnega dvojčka bi bilo mogoče ustvariti metodologijo za analiziranje podrobnih podatkov. Za izvedljivost in praktično uporabnost pristopa bo potrebnih več podatkov in primerov. Logične odvisnosti med sloji okviija na osnovi digitalnega dvojčka lahko bolj dosledno opišejo dejanski proces odločanja pri optimizaciji izbire programa PSS. V oblikovanem kvantitativnem pristopu za optimizacijo so vključene ocene več strokovnjakov, ki so izražene z mehkim številom in z grobim intervalom. Alternativni programi PSS so ocenjeni s sintezo pomembnosti atributov in prilagojeno metodo TOPSIS za optimizacijo izbire programa PSS. Ključne besede: produktno-storitveni sistem, izbira programa, digitalni dvojček, mehka ocena, grobi mejni interval, relacijska vektorska razdalja SI 66 *Naslov avtorja za dopisovanje: Severna Minzu univerza. Šola za strojništvo in elektrotehniko, VVenchang Street 204, Yinchuan 750021, Kitajska, lianhuili@foxmail.com Strojniški vestnik - Journal of Mechanical Engineering 66(2020)9, SI 67 © 2020 Strojniški vestnik. Vse pravice pridržane. Prejeto v recenzijo: 2020-03-11 Prejeto popravljeno: 2020-06-18 Odobreno za objavo: 2020-07-14 Numerična analiza z eksperimentalno validacijo enofaznega toka tekočine v kanalu ploščnega prenosnika toplote z jamičasto strukturo Urban Močnik12 - Bogdan Blagojevič3 - Simon Muliič14 * iUniverza v Novem mestu, Fakulteta za strojništvo, Slovenija 2Danfoss Trata d.o.o., Slovenija 3Plinovodi d.o.o., Slovenija 4SIMUTEH s.p., Slovenija Ploščni prenosnik toplote z jamičasto strukturo grelne plošče (ang. dimple pattern ali tudi micro plate heat exclianger) vsebuje, v nasprotju s prevladujočo obliko grelnih plošč z obliko ribje kosti (ang. chevron plate heat exchanger), veliko število jamic. Število jamic in njihova oblika neposredno definirajo karakteristike prenosnika toplote. Pri jamičasti strukturi lahko jasno določimo element iz katerih je sestavljen celoten volumen prenosnika toplote skozi katerega teče tok tekočine. Tak element imenujemo celica. Kljub temu, da se prenosniki toplote z jamičasto strukturo vse bolj uveljavljajo, pa je poznavanje tokovnih razmer v tem tipu prenosnika toplote slabo raziskano. Dobro poznavanje stanja toka tekočine v ploščnih prenosnikih toplote z jamičasto strukturo grelne plošče (MPHE) pa je ključnega pomena za uspešen in učinkovit razvoj novih produktov. Analitični modeli za popis karakteristik ploščnih prenosnikov toplote se razvijajo že vrsto desetletij, vendar nam ti praviloma dajejo integralne rezultate za celoten prenosnik toplote. Za poglobljeno analizo stanja je treba poznati značilnosti toka tekočine na mikro nivoju, v celici MPHE. Tok tekočine na mikro nivoju je zaradi majhnih dimenzij celic in razgibane strukture v MPHE relativno težko opazovati. Nekatere značilnosti MPHE sicer lahko povezujemo z analitičnimi rezultati, kakovostnejši popis tokovnih razmer pa omogoča aplikacija numerične dinamike tekočin (ang. computational fluid dynamics; CFD). Za doseganje zanesljivih rezultatov s pomočjo aplikacije CFD je potrebno primarno določiti ustrezno geometrijo, ki jo želimo obravnavati, uporabiti primerno računsko mrežo, ustrezen fizikalni in matematični model ter prave robne pogoje. Prav tako je pomembno, da rezultate validiramo s kakovostnimi rezultati laboratorijskih eksperimentov. V članku je predstavljena analiza pri kateri smo z numerično dinamiko tekočin in laboratorijskimi eksperimenti raziskali stanje enofaznega toka vode v MPHE. Izvedli smo več laboratorijskih meritev kanala, ki sestavlja MPHE. Zasnovali smo nov način izvedbe laboratorijskih meritev MPHE kanala s tako imenovanimi velikimi posodami. S to merilno metodo smo dobili rezultate, kijih lahko neposredno primeijamo z rezultati numeričnih preračunov. Z metodami vzvratnega inženiringa smo izdelali 3D model celice s katerim smo izdelali različne geometrije numeričnih modelov, uporabljene v analizi. Izkazalo seje namreč, daje rezultat pri različnih numeričnih modelih lahko zavajajoč v kolikor ni uporabljena numerična domena ustrezne velikosti. Numerična analiza je bila izvedena s komercialno programsko opremo za numerično dinamiko tekočin, z dvema modeloma turbulence, in sicer Realizable k-e modelom z uporabo izboljšanih stenskih funkcij (ang. Enlianced wall treatment; EWT) in k-m SST modelom. Prvi model turbulence napove nekoliko manjšo tlačno razliko, drugi model turbulence pa nekoliko večjo v primeijavi z rezultati laboratorijskih meritev. Rezultati opravljene analize sicer izkazujejo dobro ujemanje rezultatov numeričnih preračunov z rezultati laboratorijskih meritev. Izvedena je bila tudi analiza vpliva gostote mreže na numerični rezultat. Izbrana je bila računska mreža, ki predstavlja dobro razmeije med točnostjo numeričnih rezultatov in uporabljeno računsko močjo. Ugotovili smo, da se v kanalu, kljub relativno majhni hitrosti tekočine pojavi turbulentni tok, kije posledica njegove oblike. Prav tako smo ugotovili, da v kanalu MPHE nastaneta dva različna tokovna režima. Prvi je živi režim, ki prevladuje pri prenosu toplote, drugi pa mrtvi, kjer se za lotno točko pojavi recirkulacijsko območje, ki zmanjšuje površino za prenos toplote. Velikost mrtvega režima se s hitrostjo tekočine v obravnavanem območju bistveno ne spreminja. Ugotovili smo tudi, daje pri stacionarnih simulacijah bolj primerna uporaba Realizable k-e modela turbulence z uporabo EWT kot pa uporaba k-m SST turbulentnega modela. Analiza, prikazana v članku, omogoča razumevanje tokovnih razmer v kanalu MPHE, ki so ključnega pomena pri načrtovanju novih in učinkovitejših geometrij prenosnikov toplote. Prihodnje raziskave bodo zajele izvajanje časovno odvisnih simulacij in uvedbo energijskih enačb za popis prenosa toplote pri numeričnem delu ter nadgradnjo laboratorijske testne proge na način, ki bo omogočal tudi meijenje prenosa toplote v MPHE. Ključne besede: prenosnik toplote, jamičasta struktura, tlačni padec, CFD, model turbulence *Naslov avtorja za dopisovanje: Univerza v Novem mestu. Fakulteta za strojništvo. Na Loko 2,8000 Novo mesto, Slovenija, simon.muhlc@fs-unm.si SI 67 Guide for Authors All manuscripts must be in English. Pages should be numbered sequentially. The manuscript should be composed in accordance with the Article Template given above. The maximum length of contributions is 12 pages (approx. 5000 words). Longer contributions will only be accepted if authors provide justification in a cover letter. For full instructions see the Information for Authors section on the journal's website: http://en.sv-jme.eu . SUBMISSION: Submission to SV-JME is made with the implicit understanding that neither the manuscript nor the essence of its content has been published previously either in whole or in part and that it is not being considered for publication elsewhere. All the listed authors should have agreed on the content and the corresponding (submitting) author is responsible for having ensured that this agreement has been reached. The acceptance of an article is based entirely on its scientific merit, as judged by peer review. Scientific articles comprising simulations only will not be accepted for publication; simulations must be accompanied by experimental results carried out to confirm or deny the accuracy of the simulation. Every manuscript submitted to the SV-JME undergoes a peer-review process. The authors are kindly invited to submit the paper through our web site: http://ojs.sv-jme.eu. The Author is able to track the submission through the editorial process - as well as participate in the copyediting and proofreading of submissions accepted for publication - by logging in, and using the username and password provided. SUBMISSION CONTENT: The typical submission material consists of: - A manuscript (A PDF file, with title, all authors with affiliations, abstract, keywords, highlights, inserted figures and tables and references), - Supplementary files: • a manuscript in a WORD file format • a cover letter (please see instructions for composing the cover letter) • a ZIP file containing figures in high resolution in one of the graphical formats (please see instructions for preparing the figure files) • possible appendicies (optional), cover materials, video materials, etc. Incomplete or improperly prepared submissions will be rejected with explanatory comments provided. In this case we will kindly ask the authors to carefully read the Information for Authors and to resubmit their manuscripts taking into consideration our comments. COVER LETTER INSTRUCTIONS: Please add a cover letter stating the following information about the submitted paper: 1. Paper title, list of authors and their affiliations. One corresponding author should be provided. 2. Type of paper: original scientific paper (1.01), review scientific paper (1.02) or short scientific paper (1.03). 3. A declaration that neither the manuscript nor the essence of its content has been published in whole or in part previously and that it is not being considered for publication elsewhere. 4. State the value of the paper or its practical, theoretical and scientific implications. What is new in the paper with respect to the state-of-the-art in the published papers? Do not repeat the content of your abstract for this purpose. 5. We kindly ask you to suggest at least two reviewers for your paper and give us their names, their full affiliation and contact information, and their scientific research interest. The suggested reviewers should have at least two relevant references (with an impact factor) to the scientific field concerned; they should not be from the same country as the authors and should have no close connection with the authors. FORMAT OF THE MANUSCRIPT: The manuscript should be composed in accordance with the Article Template. The manuscript should be written in the following format: - A Title that adequately describes the content of the manuscript. - A list of Authors and their affiliations. - An Abstract that should not exceed 250 words. The Abstract should state the principal objectives and the scope of the investigation, as well as the methodology employed. It should summarize the results and state the principal conclusions. - 4 to 6 significant key words should follow the abstract to aid indexing. - 4 to 6 highlights; a short collection of bullet points that convey the core findings and provide readers with a quick textual overview of the article. These four to six bullet points should describe the essence of the research (e.g. results or conclusions) and highlight what is distinctive about it. - An Introduction that should provide a review of recent literature and sufficient background information to allow the results of the article to be understood and evaluated. - A Methods section detailing the theoretical or experimental methods used. - An Experimental section that should provide details of the experimental set-up and the methods used to obtain the results. - A Results section that should clearly and concisely present the data, using figures and tables where appropriate. - A Discussion section that should describe the relationships and generalizations shown by the results and discuss the significance of the results, making comparisons with previously published work. (It may be appropriate to combine the Results and Discussion sections into a single section to improve clarity.) - A Conclusions section that should present one or more conclusions drawn from the results and subsequent discussion and should not duplicate the Abstract. - Acknowledgement (optional) of collaboration or preparation assistance may be included. Please note the source of funding for the research. - Nomenclature (optional). Papers with many symbols should have a nomenclature that defines all symbols with units, inserted above the references. If one is used, it must contain all the symbols used in the manuscript and the definitions should not be repeated in the text. In all cases, identify the symbols used if they are not widely recognized in the profession. Define acronyms in the text, not in the nomenclature. - References must be cited consecutively in the text using square brackets [1] and collected together in a reference list at the end of the manuscript. - Appendix(-icies) if any. SPECIAL NOTES Units: The SI system of units for nomenclature, symbols and abbreviations should be followed closely. Symbols for physical quantities in the text should be written in italics (e.g. v, T, n, etc.). Symbols for units that consist of letters should be in plain text (e.g. ms-1, K, min, mm, etc.). Please also see: http://physics.nist.gov/cuu/pdf/sp811.pdf . Abbreviations should be spelt out in full on first appearance followed by the abbreviation in parentheses, e.g. variable time geometry (VTG). The meaning of symbols and units belonging to symbols should be explained in each case or cited in a nomenclature section at the end of the manuscript before the References. Figures (figures, graphs, illustrations digital images, photographs) must be cited in consecutive numerical order in the text and referred to in both the text and the captions as Fig. 1, Fig. 2, etc. Figures should be prepared without borders and on white grounding and should be sent separately in their original formats. If a figure is composed of several parts, please mark each part with a), b), c), etc. and provide an explanation for each part in Figure caption. The caption should be self-explanatory. Letters and numbers should be readable (Arial or Times New Roman, min 6 pt with equal sizes and fonts in all figures). Graphics (submitted as supplementary files) may be exported in resolution good enough for printing (min. 300 dpi) in any common format, e.g. TIFF, BMP or JPG, PDF and should be named Fig1.jpg, Fig2.tif, etc. However, graphs and line drawings should be prepared as vector images, e.g. CDR, AI. Multi-curve graphs should have individual curves marked with a symbol or otherwise provide distinguishing differences using, for example, different thicknesses or dashing. Tables should carry separate titles and must be numbered in consecutive numerical order in the text and referred to in both the text and the captions as Table 1, Table 2, etc. In addition to the physical quantities, such as t (in italics), the units [s] (normal text) should be added in square brackets. Tables should not duplicate data found elsewhere in the manuscript. Tables should be prepared using a table editor and not inserted as a graphic. REFERENCES: A reference list must be included using the following information as a guide. Only cited text references are to be included. Each reference is to be referred to in the text by a number enclosed in a square bracket (i.e. [3] or [2] to [4] for more references; do not combine more than 3 references, explain each). No reference to the author is necessary. References must be numbered and ordered according to where they are first mentioned in the paper, not alphabetically. All references must be complete and accurate. Please add DOI code when available. Examples follow. Journal Papers: Surname 1, Initials, Surname 2, Initials (year). Title. Journal, volume, number, pages, DOI code. [1] Hackenschmidt, R., Alber-Laukant, B., Rieg, F. (2010). Simulating nonlinear materials under centrifugal forces by using intelligent cross-linked simulations. Strojniški vestnih - Journal of Mechanical Engineering, vol. 57, no. 7-8, p. 531-538, DOI:10.5545/sv-jme.2011.013. Journal titles should not be abbreviated. Note that journal title is set in italics. Books: Surname 1, Initials, Surname 2, Initials (year). Title. Publisher, place of publication. [2] Groover, M.P. (2007). Fundamentals of Modern Manufacturing. John Wiley & Sons, Hoboken. Note that the title of the book is italicized. Chapters in Books: Surname 1, Initials, Surname 2, Initials (year). Chapter title. Editor(s) of book, book title. Publisher, place of publication, pages. [3] Carbone, G., Ceccarelli, M. (2005). Legged robotic systems. Kordič, V., Lazinica, A., Merdan, M. (Eds.), Cutting Edge Robotics. Pro literatur Verlag, Mammendorf, p. 553576. Proceedings Papers: Surname 1, Initials, Surname 2, Initials (year). Paper title. Proceedings title, pages. [4] Štefanič, N., Martinčevič-Mikič, S., Tošanovič, N. (2009). Applied lean system in process industry. MOTSP Conference Proceedings, p. 422-427. Standards: Standard-Code (year). Title. Organisation. Place. [5] ISO/DIS 16000-6.2:2002. Indoor Air - Part 6: Determination of Volatile Organic Compounds in Indoor and Chamber Air by Active Sampling on TENAX TA Sorbent, Thermal Desorption and Gas Chromatography using MSD/FID. International Organization for Standardization. Geneva. WWW pages: Surname, Initials or Company name. Title, from http://address, date of access. [6] Rockwell Automation. Arena, from http://www.arenasimulation.com, accessed on 200909-07. EXTENDED ABSTRACT: When the paper is accepted for publishing, the authors will be requested to send an extended abstract (approx. one A4 page or 3500 to 4000 characters or approx. 600 words). The instruction for composing the extended abstract are published on-line: http://www.sv-jme. eu/information-for-authors/ . COPYRIGHT: Authors submitting a manuscript do so on the understanding that the work has not been published before, is not being considered for publication elsewhere and has been read and approved by all authors. The submission of the manuscript by the authors means that the authors automatically agree to transfer copyright to SV-JME when the manuscript is accepted for publication. All accepted manuscripts must be accompanied by a Copyright Transfer Agreement, which should be sent to the editor. The work should be original work by the authors and not be published elsewhere in any language without the written consent of the publisher. The proof will be sent to the author showing the final layout of the article. Proof correction must be minimal and executed quickly. Thus it is essential that manuscripts are accurate when submitted. Authors can track the status of their accepted articles on http://en.sv-jme.eu/. PUBLICATION FEE: Authors will be asked to pay a publication fee for each article prior to the article appearing in the journal. However, this fee only needs to be paid after the article has been accepted for publishing. The fee is 380 EUR (for articles with maximum of 6 pages), 470 EUR (for articles with maximum of 10 pages), plus 50 EUR for each additional page. The additional cost for a color page is 90.00 EUR (only for a journal hard copy; optional upon author's request). These fees do not include tax. Strojniški vestnik -Journal of Mechanical Engineering Aškerčeva 6, 1000 Ljubljana, Slovenia, e-mail: info@sv-jme.eu http://www.sv-jme.eu Contents Papers Andraz LipoLt, Brane Sirok, Marko Hocevar, Lovrenc Novak: Convective Drying of Sewage Sludge Layer in Through-flow Anna TreLka, Wojciech Zórawski, Anna GóraL: Microstructure and Property Modification of Cold Sprayed Coatings Using Different Grain Sizes of Cr3C2-25(Ni20Cr) Composite Powder Frank GoLdschmidtboeing, Uwe PeLz, Karen CLaire-Zimmet, Michael WoLf, RaLf GoerLach, Peter Woias: Bristle Motion, Forces, and Related Vertical Translation for a Novel Electric Toothbrush Design Xiangyang Xu, Xupeng Fan, Peitang Wei, Baojun Yang: Research on the Lubrication Characteristics of Harmonic Gear Transmission Meshing Areas Andres Lopez-Lopez, Jose BiLLerman RobLes-Ocampo, PerLa Yazmin SeviLLa-Camacho, OrLando Lastres-DanguiLLecourt, Jesús Muniz, Bianca Yadira Perez-Sariñana, Sergio de La Cruz: Dynamic Instability of a Wind Turbine Blade Due to Large Deflections: An Experimental Validation Yan Li, Lianhui Li: Enhancing the Optimization of the Selection of a Product Service System Scheme: A Digital Twin-Driven Framework Urban Močnik, Bogdan Blagojevič, Simon Muhič: Numerical Analysis with Experimental Validation of Single-Phase Fluid Flow in a Dimple Pattern Heat Exchanger Channel 481 494 505 513 523 544 9770039248001