>ince ^ „ . ...„ 1955 r Strojniški vestnik Journal of Mechanical Engineering no. 6 year 2017 volume 63 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: Papirografika Bori, printed in 300 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 Branko Širok University of Ljubljana, Faculty of Mechanical Engineering, Slovenia International Editorial Board Kamil Arslan, Karabuk University, Turkey Hafiz Muhammad Ali, University of Engineering and Technology, Pakistan Josep M. Bergada, Politechnical University of Catalonia, Spain Anton Bergant, Litostroj Power, Slovenia Miha Boltežar, UL, Faculty of Mechanical Engineering, Slovenia Franci Čuš, UM, Faculty of Mechanical Engineering, Slovenia Anselmo Eduardo Diniz, State University of Campinas, Brazil Igor Emri, UL, Faculty of Mechanical Engineering, Slovenia Imre Felde, Obuda University, Faculty of Informatics, Hungary Janez Grum, UL, Faculty of Mechanical Engineering, Slovenia Imre Horvath, Delft University of Technology, The Netherlands Aleš Hribernik, UM, Faculty of Mechanical Engineering, 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, UL, Faculty of Mechanical Engineering, Slovenia Milan Kljajin, J.J. Strossmayer University of Osijek, Croatia Peter Krajnik, Chalmers University of Technology, Sweden Janez Kušar, UL, Faculty of Mechanical Engineering, Slovenia Gorazd Lojen, UM, Faculty of Mechanical Engineering, Slovenia Thomas Lubben, University of Bremen, Germany Janez Možina, UL, Faculty of Mechanical Engineering, Slovenia George K. Nikas, KADMOS Engineering, UK José L. Ocana, Technical University of Madrid, Spain Miroslav Plančak, University of Novi Sad, Serbia Vladimir Popovic, University of Belgrade, Faculty of Mech. Eng., Serbia Franci Pušavec, UL, Faculty of Mechanical Engineering, Slovenia Bernd Sauer, University of Kaiserlautern, Germany Rudolph J. Scavuzzo, University of Akron, USA Arkady Voloshin, Lehigh University, Bethlehem, USA Vice-President of Publishing Council Jože Balič University of Maribor, Faculty of Mechanical Engineering, Slovenia Cover: Cover presents newly developed test rig for in-situ observation of the real contact area on a submicron scale. The sequence of frames captured during the experiment is visible on the right part of the image. It shows how the real contact area between the deformable rough surface and rigid sapphire window is developed under increasing static load, schematically presented on the left part of the image. Image Courtesy: Laboratory for Tribology and Interface Nanotechnology, Faculty of Mechanical Engineering, University of Ljubljana Bogisiceva ul. 8, 1000 Ljubljana, Slovenia https://www.tint.fs.uni-lj.si ISSN 0039-2480 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://en.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 peerreview process. © 2017 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 http://www.sv-jme.eu, where you access also to papers' supplements, such as simulations, etc. Contents Strojniški vestnik - Journal of Mechanical Engineering volume 63, (2017), number 6 Ljubljana, June 2017 ISSN 0039-2480 Published monthly Papers Blaž Brodnik Žugelj, Mitjan Kalin: In-situ Observations of a Multi-Asperity Real Contact Area on a Submicron Scale 349 Musibau O. Ogunlana, Esther T. Akinlabi, Mutiu F. Erinosho: Analysis of the Influence of Laser Power on the Microstructure and Properties of a Titanium Alloy-Reinforced Boron Carbide Matrix Composite (Ti6Al4V-B4C) 363 Ahmet Yildiz, Osman Kopmaz: Control-Oriented Modelling with Experimental Verification and Design of the Appropriate Gains of a PI Speed Ratio Controller of Chain CVTs 374 Yazan Taamneh, Kahled Bataineh: Mixed Convection Heat Transfer in a Square Lid-Driven Cavity Filled with Al2O3-Water Nanofluid 383 Xiaofeng Wang, Jun Liu, Weimin Ge: A Practical Method to Detect a Transverse Cracked Rotor Using Transient Response 394 Marko Peternelj, Benjamin Bizjan, Branko Širok: The Influence of Airflow Characteristics and Accumulation Grid Velocity on the Formation of a Stone Wool Primary Layer 405 Strojniški vestnik - Journal of Mechanical Engineering 63(2017)6, 351-362 © 2017 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2017.4366 Original Scientific Paper In-situ Observations of a Multi-Asperity Real Contact Area on a Submicron Scale Blaž Brodnik Žugelj - Mitjan Kalin* University of Ljubljana, Faculty of Mechanical Engineering, Slovenia We present apparatus that allows in-situ optical measurements of the evolving real contact area between a rigid glass and a deformable Al6026 surface with 700 nm of lateral and 20 nm of vertical resolution. In previous experimental studies of multi-asperity real contact area this was investigated either with much less accuracy or did not include the full (loaded) nominal contact area, which can hinder the relevant sub-micron deformation phenomena. During experiments involving the real contact area, the contact load and asperity deformations are simultaneously measured. To show the relevance of the developed experimental procedure measurements are compared to the results calculated with the Greenwood-Williamson (GW) and a modified Abbott-Firestone (AF(H)) models, which represent the two extreme deformation-regime models. The AF(H) model shows relatively good agreement between the real contact area and the asperity deformations (< 60 %), while the GW model deviates by up to 10 times, depending on the deformation value. In contrast, the GW model shows better agreement for the relationship between the contact load and the asperity deformation (< 20 %), while the AF(H) deviates by more, approximately 30 %. The results also indicate that the real contact area is a non-linear function of the contact load, while theoretical models predict their linearity. Finally, it is demonstrated that the real contact area reaches only up to 9 % of the nominal value in the loading range up to the material yield strength, as calculated for the nominal contact parameters. Keywords: test rig, in-situ experiment, optical technique, asperities, real contact area Highlights • A test rig for in-situ experimental investigations of a multi-asperity contact under static loading conditions is presented. • The contact between a deformable flat Al6026 specimen and a transparent rigid flat is investigated. • The contact load, the asperity deformations in vertical direction and Ar/An are analyzed. • Theoretical results of Greenwood-Williamson model and modified Abbott-Firestone model are compared to the experimental results. • The comparison between the theoretical and experimental results is discussed. 0 INTRODUCTION In a conventional engineering approach the contact conditions are normally calculated with the assumption of a nominal contact area between two surfaces, which, under certain conditions, can significantly over-estimate the size of the actual contact area. A consideration of such contact conditions leads to calculations of milder contact conditions than exist in reality. The problem of the real contact area has been investigated using several different approaches, where theoretical investigations tend to predominate. Probably the best-known real-contact-area model is the Greenwood-Williamson (GW) model from 1966 [1]. Because of several general assumptions, on which the GW model is based, a number of modifications to the model have appeared. Some of them show improvements in terms of mechanical characterization [2] to [5], while numerous studies have been made on geometrical treatments [6] to [8]. Based on investigations of the real contact area it is obvious that the problem is complex and involves many aspects. One of the most important issues with a theoretical investigation of the real contact area is the geometry of the engineering surfaces and their behaviour when in contact. Based on geometric characterization of contacting surfaces contact models can be divided into four groups, i.e., statistical [9] to [11], fractal [12] and [13], multi-scale [14] to [16] and deterministic [17] to [19]. While statistical models are the most widely used because of their simple adoption and good approximation to actual engineering surfaces, they are not a very powerful tool for a detailed contact analysis since the surface properties are determined with statistical functions. For a more detailed geometrical characterization of the surfaces, different measuring tools (contact profilometer, atomic-force microscope (AFM), scanning electron microscope (SEM), optical interferometer, etc.) can be used. The question here is the resolution of the different techniques, because the same surface can be characterized differently when different techniques are used. To overcome this problem, i.e., a non-uniform surface characterization, fractal methods of surface characterization can be adopted. The main advantage of fractal methods lies in a uniform surface characterization that is independent of the measuring technique. On the other hand, some problems are also connected to fractal contact theory. Namely, the distribution of contact areas is assumed geometrically without a consideration of the actual elasticity and it predicts that lighter loads cause greater percentage of plastically deformed contacts. The shortcomings of fractal models are addressed in multi-scale models, which similar as fractal models consider multi-scale effects of the surface but incorporate more accurate deformation mechanics [14] and [16]. Another type of geometric characterization of contacting surfaces presents deterministic method [20]. With this method, first, the real engineering surface is measured with one of the surface-measuring techniques. Then, in the next step, different criteria for identifying asperity peaks can be applied [21]. According to the implemented criterion, the number of detected asperity peaks, their locations and the individual asperity-tip radii can be accurately determined. Consequently, a surface topography that is the subject of the contact analysis is set based on actual surface measurements, which makes the analysis more accurate than in the case of using statistical methods. In any case, even the deterministic approach has some limitations. The main shortcoming of the deterministic approach is related to the application of the correct criterion for identifying the asperities. Thus far, we do not know which criterion is the most accurate for the appropriate detection and characterization of asperities [21]. Every theoretical model for determining a real contact area has some advantages and disadvantages. A common disadvantage of all theoretical approaches is simplification, which relates either to assumptions about asperity deformations or their geometrical properties. For an accurate investigation of the real contact area and understanding the actual contact, an experimental approach must be adopted that makes it possible to facilitate theoretical approaches with the actual measurements. In the literature [22] to [31] several experimental techniques are described for real-contact-area measurements. Probably the most practical in-situ methods for real-contact-area measurements are based on the electrical [24] and thermal [25] contact resistance, as well as ultrasonic [26] and [27] and optical [28] and [29] principles. However, the most common technique for a real-contact-area measurement is the optical method, the only critical limitation of which is the transparency of one of the contacting bodies so the light can pass through into the contact. The main advantage of optical methods is that the asperity junction can be "seen", which provides trust and an intuitive understanding of the physical processes occurring. Several attempts to experimentally analyse a real contact area with an optical method have been reported [32] to [36]. The vast majority of research has involved the contact between an ideally smooth, rigid flat and a ball, which simulates a single asperity peak [29], [35], [37] and [38]. Multi-asperity contacts and the real contact area between two realistic and rough flat surfaces have rarely been investigated experimentally [23], [28] and [30], and no detailed study with a submicron resolution of asperity growth has been presented, to the best of our knowledge. However, multi-asperity contact analyses are of great importance since the asperities interact with each other and influence how much load each asperity will carry. This means that their simultaneous interaction and growth should be understood for the proper development of a real-contact-area model. In this work an experimental investigation of the real contact area between "ideally" smooth and rough, nominally flat surfaces with an in-situ optical technique on a submicron scale is presented. The experimentally obtained results are compared with theoretical predictions, which are calculated with two fundamental, commonly applied contact models with the two extreme deformation situations; fully elastic and fully plastic. GW model is valid for solely elastically deformed surfaces and the Abbott-Firestone (AF) model, modified with hardness criteria (see Experimental section) AF(H), assumes only plastic deformation of the contacting surfaces. Namely, the main motivation for the present investigation is not to evaluate how accurate these two contact models predict the actual behaviour, but to compare actual contact behaviour with the two possible extreme situations of the contact-deformation regimes. In fact, it is well-known and expected that the actual contact behaviour is elastic-plastic, however, several elastic-plastic models, which all have many specifics and details require in-depth comments and analyses when presented, which greatly exceeds the scope of this work. A detailed comparison and evaluation of different elastic-plastic models with the same experimental technique is addressed in [39]. 1 EXPERIMENTAL 1.1 Test Apparatus A test apparatus was designed for in-situ measurements of the real contact area between two surfaces on the submicron scale (Fig. 1). The set-up allows normal loading of a deformable, multi-asperity, nominally flat specimen against a transparent rigid window where the changes to the contact area are captured by an optical microscope, simultaneously with the deformation in the vertical direction and the contact load. Fig. 1. Test apparatus for in-situ measurements of the real contact area Fig. 2. Schematic presentation of the test rig that was developed for investigating the real contact area between the transparent rigid window and deformable, multi-asperity, nominally flat specimen Fig. 2 shows a schematic of the test apparatus, which consists of the following components: a white-light optical microscope with a charge-coupled device (CCD) camera; a rigid, "ideally" smooth, flat sapphire window; a deformable specimen; a displacement sensor, a table that is movable in the z-direction; and an actuator. A very accurate actuator is able to press the specimen against the transparent sapphire window at a constant velocity as low as 10 nm/s. The contact load is detected with a resolution of 0.1 N using a compressive-force transducer (AEP transducers, Italy), which is installed between a lever and the movable table. An accurate capacitive displacement sensor (Micro-Epsilon, Germany) with a resolution of 20 nm measures the deformation of the specimen as it is pressed against a rigid flat. In order to magnify the contact between the deformable specimen and the transparent sapphire window, high-resolution (2,560 x 1,920 pixels) images are captured with an optical microscope (Eclipse LV 150, Nikon, Japan) using 200x magnification and equipped with a CCD camera (DS-fi1, Nikon, Japan) using a frame rate of 1 Hz and a lateral resolution of 700 nm. The output signals, detected by load and displacement sensors, are captured with a data-acquisition card. Commercial software (Labview 2013, National Instruments, USA) is used for the signal processing, full control of the actuator and synchronised storage of the measured displacement, the contact load and the captured images. 1.2 Specimens and Methods A commercial 6026 aluminium alloy was used in this investigation (Table 1). Young's modulus (E), Poisson's ratio (v) and yield strength (Y) were provided by the manufacturer while we measured Vickers hardness (HV) with a Leitz Miniload microhardness tester (Leitz Miniload, Wild Leitz GmbH, Germany). The measured dimensionless Vickers hardness was converted in MPa so it could be used in theoretical analysis (see Section 1.3). First, the specimens were cut out from a rod with a circular profile. The roughness, which in this study was Ra 0.6 ^m, was obtained with a sequence of grinding and polishing steps using a surface-grinding machine (RotoPol-21 with RotoForce-3 module, Struers, Denmark). The surface-roughness parameters were measured using a stylus-tip profilometer (T8000, Hommelwerke GmbH, Germany) with a TKE100/17 probe (90° angle and 5 ^m radius) according to the DIN 4768 standard at cut-off length 0.8 mm and traversing length 4.8 mm. After the surfaces were prepared, the samples were carefully milled to form a cone shape with a circular flat top and a pre-prepared roughness. The diameter of the remaining flat surfaces was approximately 200 ^m (Fig. 3), resulting in around 0.03 mm2 of nominal contact area (An). It is especially important to note that under all the testing conditions in this study, the whole nominal contact area was captured in the objective view. Consequently, all the asperities were detected when in contact. Prior to the experiment, each specimen was cleaned with acetone and ethanol so as to remove any contaminants from the surface. Table 1. Mechanical properties of the tested aluminium 6026 specimen and the sapphire window Material Aluminium 6026 Sapphire Young's modulus (E), [MPa] 69,000 335,000 Poisson's ratio (v), [/] 0.33 0.25 Hardness (H), [MPa] 1,370 27,130 Yield strength (Y), [MPa] 320 N/A Fig. 3. Final shape of the tested aluminium specimens and a detailed view of the top surface Each experiment began by mounting the specimen on a movable table. In the next step the actuator slowly lifted the table with the specimen into contact with the rigid flat at a constant velocity of 50 nm/s. While the sample was being pressed against the sapphire window, the displacement and load sensors were measuring the vertical deformation of the sample and the contact load, respectively. For the purposes of this investigation the experiments were controlled by the contact load. The actuator pressed the specimen against the sapphire window until a preset contact load was achieved. The maximum contact load was limited to the value where the resulting nominal contact pressure is equal to the yield strength of the material (Y = 320 MPa), which in this case corresponded to 12 N. While the captured load was used directly in the contact analyses, the obtained images and the vertical displacements were post-processed to precisely define the size of the real contact area and the asperity deformations. In order to accurately determine the asperity deformations, the vertical displacement had to be measured with nanoscale resolution. However, when the deformations are controlled on this scale, any small deviations (i.e., deformations) arising from the bulk material beneath the asperities and deformations of test rig have to be considered and excluded from the asperity-deformation measurements. For this reason a polished reference sample with a roughness less than Ra 0.01 pm was prepared and tested in exactly the same way as described in the previous section. We assumed that the measured deformations of a smooth sample were only a consequence of the bulk deformations of the specimen and the deformations of the test rig (Eq. (1)). The measured deformations of the reference specimen were subtracted from the deformations of the rough sample (Eq. (2)), measured under the same loads, to obtain the actual asperity deformations (Eq. (3)) [30]. dreference d bulk + dt test rig ; dsample db bulk + dtest rig + dasperities , dasperities dsample dreference . (1) (2) (3) During the experiments the images of the microcontacts between the specimen and the sapphire were captured with a CCD camera at a frame rate of 1 Hz. When the specimen was moved into contact with the transparent window, micro-contacts between the asperities and the transparent window were observed. Moreover, because of the interference between the light reflected from the interface between the glass and air and the light reflected from the specimen, an optical fringe-pattern phenomenon was detected as well. If a threshold was applied to the obtained unprocessed image, fringes surrounding the actual contacts would produce false contacts, which would result in an overestimated contact area. To reduce this error and assess the actual contact area, each image was post-processed with the image-processing technique introduced by Krick et al. [29]. In Fig. 4 the contact between the specimen and the sapphire window is presented before and after the applied image-processing technique. Fig. 4. Image of the contact between the deformable sample and the rigid sapphire window under a load of 10 N; a) captured during the test and b) the resulting real contact area which was set by applying the image-processing technique [29] Several experiments were performed to obtain statistically relevant results. The variation among different samples at the same contact load in asperity deformations and real contact area was 8 % and 18 %, respectively. However, for the consistency of surface features, data and analyses, results from only one sample are used throughout the analyses; namely, the one that is the closest to the mean values of analysed parameters. 1.3 Theoretical Models for the Contact of Rough Surfaces In this work the results of the experiments were compared with theoretical calculations obtained using the two fundamental contact models, i.e., the GW model [1] considering only elastic deformations and the AF model [40] considering only plastic deformation, but modified with hardness criteria, as explained later on, namely the AF(H) model. Fig. 5. Determination of the real contact area at certain degree of asperity deformations by slicing the 2-D surface profile according to the definition of the AF model That is to say, in accordance with the definition of the AF model, the dependence between the real contact area Ar and the asperity deformations ra for the contact of a rigid and deformable flat was determined by truncation of the un-deformed surface (Fig. 5). Therefore, at any level of separation between the deformable and rigid flat d = z - ra, the size of the real contact area was calculated by summing the intersections between two surfaces (A,) (Eq. (4)) [40]. Ar( AF) (d) = A(d) + ■■■ + An (d) = (4) In the original work [40] authors did not determine any expression for calculating the load in fully plastic contacts. Therefore, in this work we introduce a "modified Abbott-Firestone" model, denoted as "AF(H)", which assumes that for plastically deformed contacts the mean contact pressure remains constant and is equal to the material hardness H, as suggested by Bowden and Tabor [41]. The abbreviation AF(H) is used to explicitly show that the contact model for fully plastic deformation regime is not based only on AF model, which is usually misquoted by other researchers, but considers hardness criteria for deformation, as well. Therefore, the contact loads for the fully plastic deformation regime can be calculated with the following equation: r (d) = Ar(af)(d) •H, (5) where H denotes hardness measured with microhardness tester. To precisely assess the real contact area based on the AF model, the 3-D topography of the tested surface was measured prior to the actual testing. The topography was obtained with a 3-D optical interferometer (Contour GT-K0, Bruker, USA). A 10* magnification lens was used for the measurements, which resulted in the same lateral resolution in x and y directions (Ax = Ay = 0.334 ^m) and a vertical resolution better than 0.1 nm. Additionally, median filter [42] was used to eliminate any captured "artificial asperities" due to measurement errors, which could influence the theoretical results. Fig. 6 shows the surface topography of the tested sample measured with the interferometer. A quite good agreement between nominally flat surface and Gaussian distribution is also shown in Fig. 6, which is assumed in most of statistical contact models including GW model [43]. The second contact model used in this study was the GW model. It is probably the best-known statistical model for an elastic contact between a deformable engineering surface and an ideally flat, rigid counter surface. The contact as predicted by GW model [1] is schematically shown in Fig. 7. Based on the definition of the GW model, the real contact area and the contact load were calculated with the following equations: (d ) = nNRAnjm^)(z )dz (6) 4 Pr(GW) (d) = - NERm Anfa>3/20(z )dz, (7) where 1_1-v' 1-v22 (8) In Eqs. (6) and (7), N and R denote the asperity density and the average radius of the detected asperity tips, respectively; An is the nominal contact Surface Height [nm] Fig. 6. Topography of the tested aluminium specimen, measured with the 3-D optical interferometer Deformable surface Fig. 7. Contact between the simplified deformable surface and rigid flat in accordance with the GW model [1] area between two surfaces; $(z) denotes the normal distribution function of the asperity heights; and z represents the asperity height measured from the mean line of the summit heights. Eq. (8) presents the expression for calculating the equivalent Young's modulus E, where Ex and E2 are the Young's moduli and Vj and v2 are the Poisson's ratios of the contacting surfaces. and characterize the asperities. Besides N and R, the maximum asperity peak height (Zmax), the distance between the mean of asperity heights and the mean of the surface heights (ys), the standard deviation of the asperity heights (os) and the standard deviation of the surface heights (o) were also calculated (Table 2) for a unified comparison of the experimental and theoretical results. Table 2. Surface parameters of the analysed aluminium surface with a roughness Ra 0.6 pm calculated using the 9PP criterion Asperity density (N), [1/m2] 5.02e+10 Average radius of detected asperity tips (R), [mm] 1.70e-03 The maximum asperity peak height (Zmax), [mm] 2.48e-03 Distance between the mean of asperity heights and 0.83e-03 the mean of the surface heights (ys), [mm] Standard deviation of asperity heights (omax), [mm] 0.49e-03 Standard deviation of surface heights (o), [mm] 0.85e-03 The material properties of the aluminium samples and the sapphire window were provided by the distributer and are specified in Table 1. The surface parameters for the deformable specimen were determined from the same 3-D topography as presented in Fig. 6. The asperities on the surface needed to be determined arbitrarily. Therefore, the 9-point-peak (9PP) criterion [21] was adopted to detect 2 RESULTS 2.1 Relationship between the Contact Load and the Asperity Deformations Fig. 8 shows the asperity deformations as a function of the contact load for the experimental measurements and the theoretical results obtained using the GW and AF(H) contact model. The calculated asperity deformations for the same loads were the largest for the AF(H) model, where fully plastic deformations of the contacting asperities were assumed. The smallest deformations were calculated for the surfaces that were only deformed elastically, according to the GW model. It is clear that the experimental results are in better agreement with the GW model than with the AF(H) model over the whole region. Moreover, up to 2 N, the deviations between the experimentally obtained data and the results predicted by the GW model are negligible (up 4 5 6 7 Contact load, [N] Fig. 8. Relationship between the contact load and asperity deformations, obtained experimentally and predicted using the AF(H) and GW models to 3 % at 2 N). For loads above 2 N the experimental results start to deviate more significantly, also from the results of GW model and, in terms of the trend, i.e., slope, are more similar to the predictions of the AF(H) model. At a contact load of 12 N the measured asperity deformations are approximately 20 % greater than those predicted by the GW model. In contrast, the largest discrepancy between the experimental data and the results of the AF(H) model is about 31 %, at 12 N. According to the GW model the deformation regime of a random surface can be predicted using the plasticity-index criterion [1], which is defined as: Y= (E / H) • (as / R)1/2. (9) By considering the material (Table 1) and topographic (Table 2) properties, the plasticity index W for the tested aluminium surface with a roughness Ra 0.6 ^m was 25. Therefore, our surface should undergo plastic deformation and should, consequently, be in better agreement with the AF(H) model. However, as seen in Fig 8, the experimental results deviate significantly from these theoretical predictions. Furthermore, while the AF(H) model predicts that the asperity deformations will match the mean value of the asperity heights at 12 N, this level of deformation is not achieved within the loading range for the case of experimental measurements and the GW model. 2.2 Relationship between the Asperity Deformations and Ar/An Fig. 9 shows Ar/An as a function of the asperity deformations. It is clear that the theoretical results have a similar trend to the experimentally obtained data. Three regions can be seen in Fig. 9. When the asperity deformations are below 890 nm the measured Ar/An values are larger than the theoretical predictions of the AF(H) model. The largest deviation between the experimental and theoretical results is at 670 nm, where the experimentally measured Ar/An value is 2.2 % and the Ar/An value predicted by the AF(H) model is 1.2 %. For the asperity deformations between 890 nm and 1,100 nm the differences between the experimental results and the predictions of the AF(H) model are negligible. With a further increase in the asperity deformations (i.e., > 1,100 nm) the measured Ar/An value is as much as 17 % (at 1,280 nm) smaller than the Ar/An value predicted by the AF(H) model. On the other hand, the GW model predicts roughly up to 10-times smaller Ar/An values than the experimental results across the deformation range. Based on the dependency between Ar/An and the asperity deformations we can claim that the results obtained using the AF(H) model are in better agreement with the experimental data than the predictions of the GW model, which is a contrast to the relationship between the contact load and the asperity deformations in Fig. 8. 2.3 Relationship between the Contact Load and the Ar/An Fig 10a compares the experimentally obtained relationship between Ar/An and the contact load with the results predicted by the AF(H) and GW models. In the case of the theoretical GW and AF(H) models the relationship between the contact load and Ar/An is linear over the whole range of loads. It is clear that the experimental values lie in between the results obtained 400 600 800 1000 Asperity deformations, [nm] Fig. 9. Relationship between the asperity deformations and Ar/An, obtained experimentally and predicted by the AF(H) and GW models 25 n (a) Contact load, [N] (b) 5 4 ] ^ 3 K ^ 2 1 0 7 3 4 Contact load, [N] Fig. 10. a) The results for Ar/An depending on the contact load, obtained experimentally and predicted by the AF(H) and GW models; and b) a detailed view from 10a with two separate regions for the relationship between the contact load and Ar/An using the GW and AF(H) models over the whole loading region. The experimentally obtained Ar/An is up to 5-times higher than that predicted using the GW model, while it is up to 2 or 3 times lower compared to the predictions of the AF(H) model, depending on the load. The deviations increase as the load increases. The experimentally measured Ar/An at the maximum contact load, i.e., 12 N, is 9 % (Fig. 10a). However, a closer look (see Fig 10b) at the experimental results shows two separate regions for the experimental relationships between Ar/An and the contact load. In both regions the relationship between Ar/An and the contact load is fairly linear for the GW and AF(H) models. However, up to 2 N, which corresponds to 600 nm of asperity deformation, the slope coefficient (k) for the experimental data is 1.2, and is in better agreement with the AF(H) model. In this region the GW model deviates significantly, i.e., a real contact area that is as much as a 10-times smaller. For the contact loads higher than 2 N (corresponding to an asperity deformation from 670 nm to 1,280 nm), the slope coefficient for the experimental results is reduced to about 0.7, and it becomes a little more like the GW model. This indicates that the actual real contact is not changing linearly with load, as the theoretical models would suggest, but it increases faster and becomes more plastic in nature at low loads, while at higher loads the increase in the real contact area is smaller and it is closer to elastic behaviour, which the GW model predicts. 3 DISCUSSION This study looks at the contact between two flat surfaces. The results were obtained using a test rig that allows measurements of the asperity deformations and the real contact area with a submicron resolution, i.e., 20 nm and 700 nm, respectively, and is controlled by the contact load with a resolution of 0.1 N. Moreover, the test rig, together with the test specimens, is designed to monitor the full nominal contact area, thus including every possible asperity that can carry the load. This is very important in order to control all the asperities and their deformations together with the contact load from the moment when the first contact occurs, i.e., at a nanoscale resolution. In the multi-asperity contact analyses of Azushima et al. [28], based on an optical method in a similar study to ours, the resolution of the test rig was lower and only a portion of the nominal contact area was captured during the experiment where the "unseen" asperity contacts can affect the results significantly. Therefore, we believe that a complete set of experimental contact parameters must be captured for accurate measurements, such as those developed in this test rig. The analysed specimen with a roughness Ra 0.6 ^m was tested for the full range of possible engineering loads, i.e., from almost negligible up to a macro yield stress, which was calculated to initiate at a nominal contact load of 12 N. Moreover, it should be noted that for a selected material and roughness, for the maximum contact load at yield (12 N), the asperities only deformed by 1,280 nm, while their deformations at more typical engineering loads are much less, around several hundred nm, indicating that the major part of the contact deformations is indeed on the nanoscale for the surfaces we used. This further shows the necessity for high precision, i.e., nanoscale, asperity deformation measurements. Accordingly, at this experimental precision, we clearly show that the real contact area under such static loads can actually be very small. Only 9 % of the nominal contact area was in contact at the maximum contact load at the yields stress, while for more typical loads around 0.5Y, i.e., 10 N, this was 7.5 %. So far, this has not been experimentally shown at the precision we used or for the material and topographic properties similar to ours. In [28] the Ar/An measured at the yield initiation was around 30 %; however, rougher (Ra 5.22 ^m) and much softer (Y = 120 MPa) material was analysed in that study. In this study the predictions of the fundamental GW and AF(H) contact models were also compared to the experimental results. There are two different aspects when comparing the theoretical and experimental results. Namely, when Ar/An is analysed as a function of the asperity deformations, the relationship calculated with the predominantly plastic AF(H) model is very close to the experimental results, Fig. 9. Moreover, if we look at the calculated plasticity index (Eq. (9)), which for our tested surface predicts a fully plastic deformation at a value of ¥ = 25, this relationship shows a very good agreement between this, very often used plasticity criterion [1], and actual contact behaviour in this work. Therefore, for the relationship between Ar/An and the asperity deformations, the GW model was not appropriate. The results at this point appear conventional and predictable, following the plasticity index criterion, which is indeed typically used. However, Fig. 8 clearly shows that the contact loads to achieve a certain micro-asperity deformation are much higher than those predicted by the AF(H) model. Accordingly, just the opposite of the above, when the relationships between the asperity deformations and the contact loads are analysed, the results of the GW model [1] for predominantly elastic deformations are in better agreement with the experimentally obtained relationship, Fig. 8, than the AF(H) model. This indicates that the contact loads predicted by the AF(H) model are underestimated for actual engineering contacts. Accordingly, the two models are contradictory for the parameters analysed within the same contact. Moreover, even the very often used plasticity index that predicts the asperity deformation - irrespective of the real contact area and the load - shows a contradictory prediction, since when it is correlated with the real contact area it is in agreement with the AF(H) model and "properly predicts" the plastic deformation [44]; while when correlated with the load, the prediction is "wrong", since the experimental results show a similarity with the GW model, which predicts elastic behaviour. These findings indicate that the actual contact behaviour within a multi-asperity contact cannot be accurately analysed with either of the two most commonly adopted theoretical models, i.e., the AF(H) or GW model, neither can it be properly anticipated using the plasticity index criterion. The last observation, which can be made based on the obtained results, refers to the relationship between the contact load and Ar/An. While both theoretical contact models predict an almost linear relationship between the contact load and Ar/An, the experimental results show that the asperities have two different behaviours. At lower loads and, consequently, for smaller deformations, Ar/An grows faster, Fig 10b. Therefore, the behaviour of Ar/An with respect to the contact load is up to 0.2Y (i.e. 2 N) more similar to the predictions of the AF(H) model than GW, Fig 10b. However, at higher loads and deformations, the growth of Ar/An slows down, so the trend becomes more like that of the GW model and thus more elastic. A similar trend was also observed by Jackson et al. [33] using a gold-coated rough ball. They suggested that the most plausible reason for such behaviour lies in strain hardening. From our results it seems that the material properties change with load, probably due to strain-hardening [33] and [45] or induced hydrostatic stresses [46]; however, the actual geometrical properties of the asperity [21] and [30] and interactions between neighbouring asperities [14] cannot be excluded from any influence since they also change simultaneously. This particular issue should be further analysed in detail to properly understand the mechanisms behind multi-asperity contact behaviour. 4 CONCLUSIONS 1. A test rig for in-situ experimental investigations of a multi-asperity contact under static loading conditions is presented. The apparatus allows measurements of the asperity deformations, the real contact area with a submicron resolution, 20 nm and 700 nm, respectively, and the contact load with a resolution of 0.1 N. 2. A careful specimen-preparation procedure and the optical specifications of the test rig make it possible to monitor the whole nominal contact area of the tested specimen, where each asperity can be detected when in contact. Therefore, no mechanical or topographic simplifications or assumptions are considered with the experimental approach presented in this work. 3. The comparison between the theoretical and experimental results indicates that comprehensive contact behaviour within a multi-asperity contact cannot be accurately predicted with only one of the two most commonly adopted theoretical models, e.g., AF(H) and GW, since the predictions of the AF(H) model are in better agreement with the experimental results for the relationship between the asperity deformations and Ar/An, while the GW model better predicts the relationship between the asperity deformations and the contact load. Moreover, we experimentally showed that the actual contact behaviour cannot be properly anticipated using the plasticity-index criterion either. 4. The experimental results for the aluminium 6026 specimen with a roughness Ra 0.6 ^m show that at macro yield initiation only about 9 % of the nominal contact area was in contact. Moreover, two distinctive regions for Ar/An growth with respect to the contact load were observed from the experimental results. In contrast, the GW and AF(H) models predict a nearly linear relation. Therefore, an experimental analysis of the real contact area on the submicron level is crucial for any in-depth investigation of the actual contact behaviour. A comparison between the experimental and theoretical results indicates discrepancies between the actual contact behaviour and the theoretical analysis. We believe that both the actual material and the topographic properties cause the deviations between theory and actual behaviour. Therefore, a further in-depth investigation of the actual contact behaviour is necessary in order to properly understand the mechanisms in multi-asperity contacts. 5 ACKNOWLEDGEMENTS The authors wish to acknowledge Slovenian Research Agency of the Republic of Slovenia for financial support. 6 REFERENCES [1] Greenwood, J., Williamson, J. (1966). Contact of nominally flat surfaces. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, vol. 295, no. 1442, p. 300-319, D0l:10.1098/rspa.1966.0242. [2] Zhao, Y., Maietta, D.M., Chang, L. (1999). An asperity mlcrocontact model incorporating the transition from elastic deformation to fully plastic flow. Journal of Tribology, vol. 122, no. 1, p. 86-93, D0I:10.1115/1.555332. [3] Jackson, R.L., Green, I. (2006). A statistical model of elasto-plastic asperity contact between rough surfaces. Tribology International, vol. 39, no. 9, p. 906-914, D0I:10.1016/j. triboint.2005.09.001. [4] Kogut, L., Etsion, I. (2002). Elastic-plastic contact analysis of a sphere and a rigid flat. Journal of Applied Mechanics, vol. 69, no. 5, p. 657-662, D0I:10.1115/1.1490373. [5] Liu, M., Proudhon, H. (2016). Finite element analysis of contact deformation regimes of an elastic-power plastic hardening sinusoidal asperity. Mechanics of Materials, vol. 103, p. 78-86, D0I:10.1016/j.mechmat.2016.08.015. [6] Jeng, Y.-R., Peng, S.-R. (2005). Elastic-plastic contact behavior considering asperity interactions for surfaces with various height distributions. Journal of Tribology, vol. 128, no. 2, p. 245-251, D0I:10.1115/1.2162557. [7] Rostami, A., Jackson, R.L. (2013). Predictions of the average surface separation and stiffness between contacting elastic and elastic-plastic sinusoidal surfaces. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 227, no. 12, p. 1376-1385, D0I:10.1177/1350650113495188. [8] Zhao, B., Zhang, S., Qiu, Z. (2015). Analytical asperity interaction model and numerical model of multi-asperity contact for power hardening materials. Tribology International, vol. 92, p. 57-66, D0I:10.1016/j.triboint.2015.05.027. [9] Beheshti, A., Khonsari, M. (2012). Asperity micro-contact models as applied to the deformation of rough line contact. Tribology International, vol. 52, p. 61-74, D0I:10.1016/j. triboint.2012.02.026. [10] Buczkowski, R., Kleiber, M. (2009). Statistical models of rough surfaces for finite element 3D-contact analysis. Archives of Computational Methods in Engineering, vol. 16, no. 4, p. 399424, D0I:10.1007/s11831-009-9037-2. [11] Xu, Y., Jackson, R.L. (2017). Statistical models of nearly complete elastic rough surface contact-comparison with numerical solutions. Tribology International, vol. 105, p. 274291, D0I:10.1016/j.triboint.2016.10.003. [12] Majumdar, A., Bhushan, B. (1991). Fractal model of elastic-plastic contact between rough surfaces. Journal of Tribology, vol. 113, no. 1, p. 1-11, D0I:10.1115/1.2920588. [13] Ciavarella, M., Delfine, V., Demelio, G. (2006). A "re-vitalized" Greenwood and Williamson model of elastic contact between fractal surfaces. Journal of the Mechanics and Physics of Solids, vol. 54, no. 12, p. 2569-2591, D0I:10.1016/j. jmps.2006.05.006. [14] Jackson, R.L., Streator, J.L. (2006). A multi-scale model for contact between rough surfaces. Wear, vol. 261, no. 11-12, p. 1337-1347, D0I:10.1016/j.wear.2006.03.015. [15] Wilson, W.E., Angadi, S.V., Jackson, R.L. (2010). Surface separation and contact resistance considering sinusoidal elastic-plastic multi-scale rough surface contact. Wear, vol. 268, no. 1-2, p. 190-201, D0I:10.1016/j.wear.2009.07.012. [16] Gao, Y.-F., Bower, A.F. (2006). Elastic-plastic contact of a rough surface with Weierstrass profile. Proceedings of the Royal Society A, vol. 462, no. 2065, p. 319-348, D0l:10.1098/ rspa.2005.1563. [17] Belghith, S., Mezlini, S., BelhadjSalah, H., Ligier, J.-L. (2010). Modeling of contact between rough surfaces using homogenisation technique. Comptes Rendus Mécanique, vol. 338, no. 1, p. 48-61, D0I:10.1016/j.crme.2009.12.003. [18] Jackson, R.L., Green, I. (2011). On the modeling of elastic contact between rough surfaces. Tribology Transactions, vol. 54, no. 2, p. 300-314, D0I:10.1080/10402004.2010.542277. [19] Buchner, B., Buchner, M., Buchmayr, B. (2009). Determination of the real contact area for numerical simulation. Tribology International, vol. 42, no. 6, p. 897-901, D0I:10.1016/j. triboint.2008.12.009. [20] Pasaribu, H., Schipper, D. (2005). Application of a deterministic contact model to analyze the contact of a rough surface against a flat layered surface. Journal of Tribology, vol. 127, no. 2, p. 451-455, D0I:10.1115/1.1866163. [21] Kalin, M., Pogacnik, A. (2013). Criteria and properties of the asperity peaks on 3D engineering surfaces. Wear, vol. 308, no. 1-2, p. 95-104, D0I:10.1016/j.wear.2013.09.010. [22] Ovcharenko, A., Halperin, G., Etsion, I., Varenberg, M. (2006). A novel test rig for in situ and real time optical measurement of the contact area evolution during pre-sliding of a spherical contact. Tribology Letters, vol. 23, no. 1, p. 55-63, D0I:10.1007/s11249-006-9113-9. [23] Kucharski, S., Starzynski, G. (2014). Study of contact of rough surfaces: Modeling and experiment. Wear, vol. 311, no. 1-2, p. 167-179, D0I:10.1016/j.wear.2014.01.009. [24] Sick, J.-H., Ostermeyer, G.-P. (2007). In situ measurement of contact area in coated surfaces. Computer Methods and Experimental Measurements for Surface Effects and Contact Mechanics VIII, (De Hosson, J.T.M., Brebbia, C.A., Nishida, S.I., eds.) vol. 55, p. 259-270, D0I:10.2495/SECM070251. [25] Bahrami, M., Yovanovich, M.M., Culham, J.R. (2005). Thermal contact resistance at low contact pressure: Effect of elastic deformation. International Journal of Heat and Mass Transfer, vol. 48, no. 16, p. 3284-3293, D0I:10.1016/j. ijheatmasstransfer.2005.02.033. [26] Baltazar, A., Kim, J.Y., Rokhlin, S.I. (2006). Ultrasonic determination of real contact area of randomly rough surfaces in elastoplastic contact. Revista Mexicana de Física, vol. 52, no. 1, p. 37-47. [27] Aymerich, F., Pau, M., Ginesu, F. (2003). Evaluation of nominal contact area and contact pressure distribution in a steel-steel interface by means of ultrasonic techniques. JSME International Journal Series C Mechanical Systems, Machine Elements and Manufacturing, vol. 46, no. 1, p. 297-305, D0I:10.1299/jsmec.46.297. [28] Azushima, A., Kuba, S., Tani, S., Olsson, D.D. (2006). Direct observation of asperity deformation of specimens with random rough surfaces in upsetting and indentation processes. Wear, vol. 260, no. 3, p. 258-264, D0I:10.1016/j.wear.2005.04.022. [29] Krick, B.A., Vail, J.R., Persson, B.N.J., Sawyer, W.G. (2012). Optical in situ micro tribometer for analysis of real contact area for contact mechanics, adhesion, and sliding experiments. Tribology Letters, vol. 45, no. 1, p. 185-194, D0l:10.1007/ s11249-011-9870-y. [30] Pogačnik, A., Požar, T., Kalin, M., Možina, J. (2013). A homodyne quadrature laser Interferometer for micro-asperity deformation analysis. Sensors, vol. 13, no. 1, p. 703-720, D0l:10.3390/s130100703. [31] Matsuda, K., Hashimoto, D., Nakamura, K. (2016). Real contact area and friction property of rubber with two-dimensional regular wavy surface. Tribology International, vol. 93, part B, p. 523-529, D0I:10.1016/j.triboint.2014.11.011. [32] Ovcharenko, A., Halperin, G., Etsion, I. (2008). In situ and real-time optical investigation of junction growth in spherical elastic-plastic contact. Wear, vol. 264, no. 11-12, p. 10431050, D0I:10.1016/j.wear.2007.08.009. [33] Jackson, R.L., Down, M.P., Liu, H., McBride, J.W. (2014). A comparison of the predictions of a multiscale model and optical real area of contact measurements. IEEE 60th Holm Conference on Electrical Contacts, p. 79-86, D0I:10.1109/ H0LM.2014.7031026. [34] McBride, J.W., Cross, K.C. (2008). An experimental investigation of the contact area between a glass plane and both metallic and carbon-nano-tube electrical contacts. Proceedings of the 54th IEEE Holm Conference on Electrical Contacts, p. 325-331, D0I:10.1109/H0LM.2008.ECP.63. [35] Nitta, I., Tsukiyama, Y., Tsukada, T., Terao, H. (2014). Measurement of real contact area on thermal print head using a laser microscope with a wide field of view. Tribology International, vol. 79, p. 162-173, D0I:10.1016/j. triboint.2014.06.005. [36] Eguchi, M., Shibamiya, T., Yamamoto, T. (2009). Measurement of real contact area and analysis of stick/slip region. Tribology International, vol. 42, no. 11-12, p. 1781-1791, D0I:10.1016/j. triboint.2009.04.046. [37] Ovcharenko, A., Halperin, G., Verberne, G., Etsion, I. (2007). In situ investigation of the contact area in elastic-plastic spherical contact during loading-unloading. Tribology Letters, vol. 25, no. 2, p. 153-160, D0l:10.1007/s11249-006-9164-y. [38] Jamari, J., Schipper, D.J. (2007). Plastic deformation and contact area of an elastic-plastic contact of ellipsoid bodies after unloading. Tribology International, vol. 40, no. 8, p. 13111318, D0I:10.1016/j.triboint.2007.02.015. [39] Brodnik Zugelj, B., Kalin, M. (2017). Submicron-scale experimental and theoretical analyses of multi-asperity contacts with different roughnesses. Tribology International, (communicated). [40] Abbott, E., Firestone, F. (1933). Specifying surface quality: a method based on accurate measurement and comparison. Mechanical Engineering, vol. 55, p. 569-572. [41] Bowden, F., Tabor, D. (1967). Friction and Lubrication, Methuen & Co. Ltd., London. [42] Rajni, R., Anutam, A. (2014). Image denoising techniques - An overview. International Journal of Computer Applications, vol. 86, no. 16, p. 13-17, D0I:10.5120/15069-3436. [43] Lee, C.-H., Eriten, M., Polycarpou, A.A. (2010). Application of elastic-plastic static friction models to rough surfaces with asymmetric asperity distribution. Journal of Tribology, vol. 132, no. 3, p. 031602, D0I:10.1115/1.4001547. [44] Jamari, J., Schipper, D.J. (2006). Experimental investigation of fully plastic contact of a sphere against a hard flat. Journal of Tribology, vol. 128, no. 2, p. 230-235, D0I:10.1115/1.2164470. [45] Chatterjee, B., Sahoo, P. (2012). Effect of strain hardening on elastic-plastic contact of a deformable sphere against a rigid flat under full stick contact condition. Advances in Tribology, vol. 2012, p. 1-8, D0I:10.1155/2012/472794. [46] Krithivasan, V., Jackson, R.L. (2007). An analysis of three-dimensional elasto-plastic sinusoidal contact. Tribology Letters, vol. 27, no. 1, p. 31-43, D0I:10.1007/s11249-007-9200-6. Strojniški vestnik - Journal of Mechanical Engineering 63(2017)6, 363-373 © 2017 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2016.4159 Original Scientific Paper Analysis of the Influence of Laser Power on the Microstructure and Properties of a Titanium Alloy-Reinforced Boron Carbide Matrix Composite (Ti6Al4V-B4C) Musibau O. Ogunlana* - Esther T. Akinlabi - Mutiu F. Erinosho University of Johannesburg, Department of Mechanical Engineering Science, South Africa Laser Metal Deposition (LMD) process Is a means of producing metal composites with the aid of a laser beam, ejected onto the substrate with the participating powder and fused together after solidification. In this research work, Ti6Al4V alloy is fused with 20 wt % of B4C in order to form metal matrix composites (MMCs). Using the Ytterbium Fibre Laser System powdered at 3000 W, the laser powers were varied between 800 W and 2400 W while all other supporting process parameters were kept constant. The deposited Ti6Al4V-B4C composites were characterized through the surfacing microstructure, microhardness and dry sliding wear. The microstructural properties of the deposited samples were profound, with a Widmanstatten structure of a-Ti, P-TI and (a+P) Ti phases. The microhardness tests revealed that the composites deposited with a laser power of 2000 W exhibited the highest hardness value and standard deviation of HV 445 ± 61. Furthermore, characterisation revealed that the sample produced with the laser power of 800 W had the lowest wear loss and wear rate of 35.2 x 10-3 mm3 and 6.42 x 10-4 mm3/Nm. However, the motivation for this work is to improve the material properties of the Ti6Al4V alloy for surface engineering applications. Keywords: Dry sliding wear, LMD, microhardness, microstructure, Ti6Al4V-B4C composites, X-Ray Diffraction Highlights • The presence of TiB and TiB2 compounds was detected between the interface of the Ti6Al4V-B4C composite and substrate. • The microstructures were mostly characterised with Widmanstatten structures of a-Ti, P-Ti and (a+P) Ti phases. • The wear loss and wear rate of 35.2 x 10-3 mm3 and 6.42 x 10-4 mm3/Nm were achieved. • The optimized process parameters were established between the laser power of 2000 W and 2400 W due to the defect-free deposition output. 0 INTRODUCTION Titanium and its alloys are used for highly demanding applications, such as in the fabrication of some of the most critical and highly stressed civilian and military aircraft parts, chemical processing, automobile industries, nuclear power plants, food processing plants and oil refinery heat exchangers. This wide range of applications has been attributed to its excellent properties, such as low density, high specific strength, heat resistance, corrosion resistance, low temperature resistance, and excellent biocompatibility. Several suggestions have been made that the physical and mechanical properties of titanium can be improved through the integration of reinforcing compounds using the principle of metal matrix composites (MMCs), since they combine the properties of ceramics and metals to produce good shear strength, high temperature strength and elevated hardness composites [1]. Titanium carbide has been agglomerated with Ti6Al4V alloy to improve its wear properties [2]. Different coatings such as plasma-sprayed Al-bronze have been applied to the surface of titanium to improve the operational life [3]. The microstructural behaviour of titanium alloy is an important aspect to be studied with the establishment of new phases [4]. In contrast, laser metal deposition (LMD) is an additive manufacturing (AM) technique that serves as a recommended technique for processing titanium and its alloy, since it addresses most of the problems of the traditional manufacturing methods [5]. This AM technology is also a promising aerospace manufacturing technique due to its potential of reducing the buy-to-fly ratio and repairing high valued parts [6]. Among the various methods of coating the surface of materials, which includes chemical vapour deposition, physical vapour deposition, spraying, etc., the LMD process is believed to have a greater advantage over other methods of coating [7]. Some of the advantages of using the LMD process include the ability to produce parts directly from a 3-dimensional computer aided design (CAD) model of the part [8] with the required surface coating in one step [9]; and its ability to be used for repair of existing worn out parts that were not repairable in the past [10]. Despite the significant progress in this field, there are still some of technical challenges that need improvement, which includes material characterization and availability [11]. The adopted approach of layer-by - layer of building a component has been attributed to the laser powder deposition technique. However, homogeneous structure, production of complex geometries and improved mechanical properties are the possible features attributed to the CAD systems [12]. In the same vein, complex parts can be produced as one single component through the LMD process, without requiring later assembly [13]. The mechanical properties of titanium carbide and titanium alloy (TiC/Ti6Al4V) composites were improved with the addition of different percentages of boron. However, the refinement mechanism of boron was ascribed to the combined effects of the nucleation growth at the super-cooled zone [14]. In addition, this TiC is regarded as the refractory ceramic material that has been widely used to improve the mechanical properties of titanium and its alloy. The thermal oxidation treatment carried out on Ti6Al4V/10 vol.% TiC composites have significantly improved its wear resistance properties due to the formation of hard oxide layer [15]. A different approach can also be used for the deposition process such as the direct laser fabrication techniques and still give good bonding [16], or by using a vacuum induction melting furnace for the fabrication process [17]. Thus, generally, laser material processing has developed a significant advantage in the field of engineering [18]. Furthermore, research has been conducted on laser surface re-melting, laser alloying, and laser cladding to improve the surface properties of many kinds of metals. It is noteworthy that the coatings formed by laser cladding process show dense microstructures and exhibit strong metallurgical bonding with the associated substrate [19], and these would always be achieved if the process parameters employed are well understood. However, the clad width, height, and geometry need to be put in place for complete deposition [20]. In the group of the most important non-metallic hard materials, boron carbide, B4C has played a role in the MMC of the primary alloy. It is currently used in high-technology industries fast-breeders, lightweight armours and high-temperature thermoelectric conversion due to such properties as high melting point, outstanding hardness, good mechanical properties, low specific weight and great resistance to chemical agents [21]. It is an extremely promising material for a variety of applications that require elevated hardness, good wear and corrosion resistance [22] and [23], and the highlighted areas of application are due to its excellent properties [24]. B4C also finds application in the nuclear industry and high-temperature thermo-electric conversion, but there is a restriction to its wide industrial application because of low strength (about 200 MPa to 400 MPa), low fracture toughness (2 MPaWm to 3 MPaWm), as well as poor sinterability that results from its low self-diffusion coefficient. Grain-refined boron-modified Ti6Al4V alloy has shown significant improvement in strength, stiffness, fatigue resistance and fracture toughness when fused together. However, the enhanced formability of the boron-modified alloy during large deformation without cracking as opposed to the normal alloy was reviewed. The microstructural evolutions of the specimen during the hot deformation were also revealed to have significant influence on the mechanical properties of the final product. However, the understanding of the processing-microstructural relationship during the thermo-mechanical processing of boron-modified Ti6Al4V alloy was discovered to be of great significant [25] and [26]. The failures of materials vary, and it is more noticeable in complex laminates and composites than in the homogeneous materials [27]. The wear properties of laser-deposited titanium boride-and carbide-reinforced composite were investigated. Excellent abrasive and adhesive wear resistance under sliding wear test conditions were achieved. In addition, the strength and hardness of the surface of samples were significantly enhanced by the in situ formed compounds [28]. Premixed titanium and boron were concocted in the weight ratio of 8:1 with a sodium silicate solution. The mixtures were laser deposited onto a clean Ti6Al4V alloy substrate. Hard ceramic compounds of TiB and TiB2 were formed in the deposited composite. In addition, the excellent wear resistance of the laser-deposited samples was also discerned. However, it was finally reported that the load-bearing capability of the substrate has a strong support for the deposit [29]. Despite the fascinating properties Ti6Al4V alloy has, its wear resistance properties need to be improved. This poor wear characteristic of the alloy has inspired the addition of 20 weight per cent (wt %) of B4C to upgrade its wear resistance. The laser powers were varied between 800 W and 2400 W respectively while keeping the scanning speed, the powder flow rate and gas flow rate constant throughout the experiments. The microstructures, the microhardness as well as the wear properties of the laser deposited Ti6Al4V-B4C have been characterized. 1 EXPERIMENTAL PROCEDURE The LMD process was achieved with a 3.0 kW Ytterbium fibre laser system at the Council for Scientific and Industrial Research, (CSIR), South Africa. A Kuka robot attached with a nozzle are used to carry out the deposition process. A rectangular Ti6Al4V alloy substrate block with dimensions of 102 mm x 102 mm x 7 mm is prepared for the LMD of the Ti6Al4V-B4C. The substrate was sandblasted and cleaned under tap water prior to the LMD coating. The Kuka robot and the powder feeder are in connection with the laser system. The powder feeder has two cylindrical glass jars (hopper) where the powders are kept. A three-way nozzle is attached to the end of the Kuka arm. It is through the nozzle that the powders and the laser beam were ejected onto the substrate. The powders are sucked (with aid of gas) into the nozzle passing through hose connections; and were delivered into the melt pool that was created by the laser beam. Ti6Al4V alloy and B4C powders are the two powders used for the experiment. Both powders were supplied by the Alfa Aesar Company in Germany. The particle size of the boron carbide powder is between 22 ^m to 59 ^m, while the titanium alloy powder's particle size ranges between 45 ^m to 90 ^m. Fig. 1 shows the scanning electron microscope (SEM) morphology and the energy dispersion spectroscopy (EDS) analysis of the Ti6Al4V alloy and B4C powders. Fig. 1. SEM morphology and EDS analysis; a) Ti6Al4V powder; and b) B4C powder The SEM morphologies of the Ti6Al4V powder are spherical with dissimilar sizes. The major peak observed in the spectrum is Ti. The SEM morphologies of the B4C powder appear amorphous and inform of stone pebbles. The major peaks observed in the B4C spectrum are boron (B) and carbide (C). Table 1 shows the experimental matrix used for the deposition of the Ti6Al4V-B4C composite. Table 1. Experimental matrix Sample Laser power Scanning Powder flow rate [rpm] name [W] speed [m/min] Ti6Al4V B4C S1 800 1.0 3.2 0.8 S2 1000 1.0 3.2 0.8 S3 1200 1.0 3.2 0.8 S4 1400 1.0 3.2 0.8 S5 1600 1.0 3.2 0.8 S6 1800 1.0 3.2 0.8 S7 2000 1.0 3.2 0.8 S8 2200 1.0 3.2 0.8 S9 2400 1.0 3.2 0.8 The depositions were made in the ratio of 4 to 1 based on the 80 wt.% for titanium alloy (Ti6Al4V) and 20 wt.% for boron carbide (B4C) as shown in Table 1. The laser powers used were varied between 800 W and 2400 W at an interval of 200 W. The scanning speed of 1 m/min, gas flow rate of 2 l/min, powder flow rate of 3.2 rpm for Ti6Al4V alloy and 0.2 rpm for B4C were kept throughout the experiment. The beam diameter or spot size of 4 mm was employed to create the melt pool on the substrate. A focal distance of 12 mm is maintained between the nozzle tip and the substrate. During the laser surface melting process, the powders were dissolved into the melted pool, leading to the alloying of the sample's surface. To protect the melt pool as well as the deposit from oxidation during laser deposition process, argon gas was initiated at a pressure of 2 MPa to provide shielding. Nine single deposits were made on the substrate and labelled from S1 to S9. Samples were cut laterally from the deposit for further characterization. 1.1 Microstructure The laterally cut samples from S1 to S9 were mounted in resin prior to micro structural investigation. The mounted samples were ground, polished, and etched according to the Struers note of standard metallographic preparation of titanium [30]. The etchant was prepared using the Kroll's reagent with 100 ml H2O; 1 ml to 3 ml HF and 4 ml to 6 ml HNO3. The samples were studied under the Optical Microscope (OM), using the Olympus BX51M and the SEM using the TESCAN instrument with Vega TC software. where, hf = wear depth, Rf = radius of the two spherical ends, W = width, Vf = wear volume, Ls = stroke length. The three highlighted equations are used simultaneously for the wear volume calculation. 1.2 Microhardness Test The microhardness characterization was performed on a digital Vickers hardness tester. The etched samples were polished again for the hardness test. Indentations were carried out according to ASTM E384-11e1 [31], using a load setting of 500 g and a dwell time of 15 seconds. Ten indentations were taken on the samples laterally and the distance of 10 ^m between each indentation was taken into consideration. 1.3 Dry Sliding Wear Test The dry sliding wear tests were conducted using ballon-disc tribometer equipment: a Universal Micro Materials Tester (UMT-2) which operates with a linear reciprocating motion drive. The wear tester was produced by the Centre for Tribology (CETR) Incorporated, USA. A tungsten carbide ball of 10 mm in diameter under the spindle is designed to rub over the surface of the samples to create the wear. A normal load of 25 N and a constant stroke length of 2 mm were used. The frequency and the speed of the reciprocating spindle were maintained at 5 Hz and 5 mm/s respectively. Each deposited Ti6Al4V-B4C sample was affixed to a non-moving table for the sliding operation. This was carried out according to the standard, ASTM G133-05 [32]. Evaluating the wear loss or volume, the measurements from the wear tester as well as the width measurement on the wear track from the SEM are taken into consideration. However, the wear volumes as well as the wear rates of the linearly operated deposited composite under the normal load application and constantly controlled speed were calculated according to Archard's wear model of Eqs. (1) to (3), respectively [33]. h = R R2 - Rf =- 4hf + W2 (1) (2) Vf = Ls f ( TTr \ Rf arcsin 2 Rf v f J - W(Rf - Af ) -J Af (3Rf - hf ) ' (3) 2 RESULTS AND DISCUSSION 2.1 Microstructural Evaluation The microstructural evaluation of the Ti6Al4V alloy substrate and the laser deposited Ti6Al4V-B4C composites are presented in this section. The microstructure of the substrate is characterized by the lighter structures called the alpha phase (a-phase) and the darker structures called the beta phase (P-phase). The micrographs of the laser-deposited Ti6Al4V-B4C composites are characterized by three different zones: The deposit zone, the fusion zone, and the heat affected zone. Fig. 2 shows the micrographs of samples S4, S5, S6, S7, S8 and S9 deposited with the laser powers of 1400 W, 1600 W, 1800 W, 2000 W, 2200 W and 2400 W, respectively. The micrographs of Figs. 2a to c of samples S4 to S6 are characterized by poor bonding which is observed at the right side of the deposit. The regions showing red ellipses were delaminated due to the embedded non-melted B4C in those regions. These occur between the deposit zone and the fusion zone at the laser powers of 1400 W, 1600 W, and 1800 W. Obviously in the spotted region, larger particle sizes and the higher melting temperature of B4C could also create the void at the interface. The entrapped B4C were not fully melted before solidification with the laser powers used. A further increase in the laser power gives good laminate as observed in Figs. 2d to f. Significantly, basket wave-like of alpha titanium (a-Ti) lamella and beta phases were formed on the deposited composites. They were found growing epitaxial and elongated towards the fusion zone and the heat-affected zone (HAZ). The HAZ, in contrast, heated up during the deposition process and acts as a heat sink. The energy density of the laser power and the density of the B4C powder initiated the elongated grains of the a-Ti lamella and the P-phase respectively. These have made the grains to extend below the fusion zone. The region under the main deposit is known as the laser melt zone, and this increases with greater energy input. However, this has a significant effect on the cooling rate during solidification [34]. Fig. 3 shows the SEM micrographs of sample S8 deposited at a laser power of 2200 W. Fig. 3b shows the enlarged region of the circle in Fig. 3 a, and this was observed to be non-melted B4C. Most of these have 2 4 Fig. 2. Micrographs of laser deposited samples S4, S5, S6, S7, S8 and S9 deposited with the laser powers of a) 1400 W, b) 1600 W, c) 1800 W, d) 2000 W, e) 2200 W and f) 2400 W created voids in the deposited samples towards the dilution zone. TiB and TiB2 are very much observed in this region. Moreover, this tends to improve the martensitic phase of the deposit. Macroscopic banding was observed in the deposit of sample S7 deposited with a laser power of 2000 W and scanning speed of 1.0 m/min. The 20 wt. % of B4C introduced into the deposit has enhanced the properties of the primary alloy. The level of homogeneity occurs at a higher laser power. Thus, this property is characterized by a martensitic microstructure that shows a significant microhardness increase and corrosion resistance improvement [19]. All the samples deposited between the laser powers of 800 W and 1800 W were characterized with significance of pores and poor bonding, which is due to the lower laser power used. This behaviour has been traced to the presence of non-melted B4C powder that was observed after Fig. 3. SEM micrographs of sample S8 deposited at laser power of2200 W the solidification of the deposits. Basket weave-a Ti has, however, been known to benefit the mechanical properties of titanium alloy. The deposition width as well as the HAZ width increased as the laser powers was increased. The beam diameter of 4 mm used has also enhanced the volume of deposited composites. The EDS analyses were conducted on the defect free composites in order to show the elements that are present in the analysed region. From the spectra observed, titanium exhibits the highest peak. The presence of Aluminium, Vanadium and Carbon was also indicated in the spectra of both samples. The microstructures of the deposited composites were characterized with macroscopic banding, martensitic structures and intermetallic (a+P)-phase of titanium alloy. Sample S1 to S6 were characterized by a lack -of -fusion and poor bonding, which create a void of non-melted B4C powder after the solidification process and thereby reduces the ductile properties of the grain structures. Samples S7, S8 and S9 deposited at laser powers of 2000 W, 2200 W and 2400 W, respectively, were selected to be the best with elongated columnar grains and Widmanstatten grain structural properties. The analyses of the x-ray diffraction (XRD) for the laser-deposited composites are presented in the Fig. 4 from samples S1 to S3, respectively. 300 3BO J (a) 1 I teta [dig] Fig. 4. XRD spectra of the laser deposited Ti6Al4V-B4C composites from samples a) S1, b) S2 and c) S3 and at varying laser power between 800 W and 1200 W The XRD patterns presented in Fig. 4 show that the phases were characterized with the increase in the intensity of the diffraction peak patterns of the coatings as the laser power increases. These results, however, suggest that the coatings were composed of slightly low traces of impurities or no irregular impurity of crystalline. Thus, the crystalline phases of the XRD patterns are marked as follow: titanium vanadium carbide, Ti0.33 V1.67 C (★); titanium boride, TiB (■); titanium, Ti (•); titanium diboride, TiB2 (▲); vanadium carbide, V2C (♦); aluminium titanium vanadium, A10.5 Ti0.5 V (▼). From the identified peaks, there was no specific shift in the phases as regards to the increase in the laser powers. Eqs. (4) to (8) show the balanced chemical equations of titanium and boron carbide. Ti + B4C ^ TiC + 4B, (4) 5Ti + B4C ^ 4TiB + TiC, (5) 3Ti + B4C ^ 2TiB2 + TiC, (6) Ti + i/2B4C ^ TiB2 + i/2C, (7) Ti + i/3B4C ^ 2/3TiB2 + !/3TiC. (8) However, the possible overall stoichiometric equation between the specific grade 5 of Ti6Al4V alloy and the ceramic, B4C used in this experiment is given in Eq. (9): 4Ti6Al4V(s) + B4C(s) ^ 2T®2(s) + + 5V2V5C(s) + 6Al4i/3TiV(s). (9) 2.2 Microhardness Profiling The microhardness results of the laser deposited composites from samples S1 to S9 are presented in Fig. 5. Fourteen indentations have been put into consideration from the top of the deposit to the substrate. From the microhardness profiling as indicated, there is an increase in the hardness values (HV) from samples S1 to S9 as the laser power increases, and at a certain instance, the hardness values show a decrease as the laser power continues to increase. There is also a rise in the hardness value of sample S7 and falls again at samples S8 and S9. From the plot, sample S1 deposited with a laser power of 800 W exhibits the lowest hardness value and standard deviation of HV 339 ± 29, while sample S7 deposited at a laser power of 2000 W displays the highest hardness value and standard deviation of HV 445 ± 61. However, samples S7 and S9 show some high values of hardness of between HV 530 and HV 545 in the main deposit, and this can be due to the indentation made on the non-melted B4C. From all indications, the substrate shows the lowest hardness values which falls below HV 339 ± 29. The improvement in the hardness values of the laser-deposited samples can be attributed to the B4C added. Thus, the variation in hardness values is set as a criterion to improve the ductility property of the composites. It is, however, observed that, an increase in the laser power leads to an increase the hardness values of the deposits. An increase-decrease-increase phenomenon was achieved from the trend of the hardness plot. 2.3 Wear characterization The SEM images of the worn scar for the respective coatings as shown in Figs. 6 a to e from samples S1 to S6 at low magnification are presented in this section. The rubbing together of the tungsten ball against the surface of the laser-deposited samples has initiated the wear scar. The worn surfaces are characterized by severe wear with plough and ridge; and wear debris 10 12 Number of indentations Fig. 5. Microhardness profiling of the laser-deposited Ti6Al4V-B4C composites as observed is apparently in all the wear samples. However, the worn surfaces are pronounced by elliptical groove shapes. Samples S1 and S2 deposited with laser powers of 800 W and 1000 W show mildly worn surfaces in comparison to other samples. The values of the wear track radius, wear width, stroke length are measured under the SEM in order for the wear volume to be evaluated. The sliding friction exists between the surfaces of the Ti6Al4V-B4C composites and the tungsten carbide ball, which was initiated by the force in action. Wear occurs due to the mechanical action and created an oval shaped groove on the surface of the composites. The wear volume was instigated as a result of the frictional force and the abrasive wear between the tungsten ball and the composite's surfaces. A thicker ridge produced at the plough way shows the extent of the wear volume. Fig. 7 shows the plots of the wear depth (hf), wear width (w) and stroke length (Ls) for all the samples, while Fig. 7 shows the histogram plot of the wear volume (Vf) and the wear rate (K). Figs. 7 and 8 clearly present the plots of the wear depth, the wear width, the stroke length, the wear Fig. 6. SEM images of the wear scar on the deposited Ti6Al4V- B4C composites from samples from S1 to S6 Samples Fig. 7. Plots of the wear depth, wear width and stroke length of the laser-deposited Ti6Al4V-B4C composites Samples Fig. 8. Graph of the wear volume and wear rate of the laser-deposited Ti6Al4V-B4C composites volume and the wear rates of all the laser deposited samples. Among the deposited composites, sample S6 deposited at a laser power of 1800 W has the lowest wear depth with 74.6 ^m while sample S8 deposited at a laser power of 2200 W exhibits the highest wear depth with 646.3 ^m. Thus, sample S6 showed profound ductile property than sample S8, which could be due to more heat input associated with the deposition of sample S8 when compared to sample S6. In contrast, more non-melted particles of B4C in the melt pool during deposition may have caused the wear depth of sample S6 to have the lowest wear depth value among the rest of the samples. Furthermore, from the histogram chart, both the wear volume and the wear rate were characterized by irregular behaviour, as observed in the trend. From all the deposited composites observed, sample S1 deposited at a laser power of 800 W has the lowest wear volume and wear rate of 35.2 x 10-3 mm3 and 0.642 x 10 3 mm3/Nm, while sample S7 deposited at a laser power of 2000 W possesses the highest wear volume and wear rate of 93.3 x 10-3 mm3 and 2.62 x 10 3 mm3/Nm, respectively. The depth and width of wear have a significant role to play in determining the wear volume. However, the major role is attributed to the width and length of the wear track. From the histogram plot, the substrate shows the highest wear volume of 569.9 x 10-3 mm3 and the occurrence was as a result of the larger wear width of 1844 ^m exhibited by the substrate. In other words, all the deposited Ti6Al4V-B4C samples show good wear volume results in comparison with the substrate. 3 CONCLUSIONS The improvement in the surface properties of the Ti6Al4V-B4C composites deposited through the LMD process was significantly achieved in this study. The surface effects of varying the laser power on the microstructure, the microhardness, and the wear were extensively studied and concluded as follows: • From the micro structural analyses, sample S8 deposited with the laser power of 2200 W and scanning speed of 1 m/min revealed a defect-free surface and good intermetallic phase of (a+P) alloy with boron carbide (B4C) particulates. • Fine globular primary alpha and martensite structures changed to thick and coarse structures as the laser power was increased. • Sample S1 deposited with a laser power of 800 W exhibited the lowest hardness value and standard deviation of HV 373 ± 48 while sample S7 deposited at a laser power of 2000 W displayed the highest hardness value and standard deviation of HV 445 ± 61. • The XRD patterns of the laser deposited composites revealed the presence of intermetallic compounds of TiB and TiB2, which were formed at the interface. • Both wear volume and wear rate for the deposited composites were calculated using the proposed Archard's wear equation. Sample S7 deposited with a laser power of 2000 W demonstrated the highest wear volume and wear rate of 93.3 x 10-3 mm3 and 2.62 x 10-3 mm3/Nm. In the same vein, they both experienced irregular increase-decrease-increase phenomenon as the laser power was increased. Based on the range of experimental matrix or process parameters adopted, this research study has been able to establish the optimized process parameters between the laser power of 2000 W and 2400 W due to the defect-free deposition output. However, this research work can be recommended for the industries dealing with aerospace parts and nuclear device applications. 4 REFERENCES [1] Jimoh, A. (2010). The Particulate-Reinforcement of Titanium Matrix Composites with Borides. PhD Thesis, The Faculty of Engineering and the Built Environment, University of the Witwatersrand, Johannesburg. [2] Mahamood, R.M., Akinlabi, E.T. (2015). Effect of laser power and powder flow rate on the wear resistance behaviour of laser metal deposited TiC/Ti6Al4V composites. 4th International Conference on Materials Processing and Characterization. Journal of Materials Today, vol. 2, no 4-5, p. 2679-2686, D0l:10.1016/j.matpr.2015.07.233. [3] Hager, C.H.Jr., Sanders, J.H., Sharma, S. (2008). Unlubricated gross slip fretting wear of metallic plasma-sprayed coatings for Ti6Al4V. Surfaces, Wear, vol. 2653, no. 3-4, p. 439-451, D0I:10.1016/j.wear.2007.11.026. [4] Lu, Y., Tang, H.B., Fang, Y.L., Liu, D., Wang, H.M. (2012). Microstructure evolution of sub-critical annealed laser deposited Ti6Al4V alloy. Journal of Materials and Design, vol. 37, p. 56-63, D0I:10.1016/j.matdes.2011.12.016. [5] Akinlabi, E.T., Akinlabi, S.A. (2015) Characterization of Functionally Graded Commercially Pure Titanium (CPTI) and Titanium Carbide (TiC) Powders. Proceedings of the World Congress on Engineering, Vol. II, London. [6] Mahamood, R.M., Akinlabi, E.T., Shukla, M., Pityana, S. (2014). Characterization of Laser Deposited Ti6A14V/TiC Composite Powders on a Ti6A14V Substrate. Lasers in Engineering. vol. 29, no. 3-4, p. 197-213. [7] Matilainen, V. (2012). Benchmarking of Laser Additive Manufacturing Process. BSc Thesis. Lappeenranta University of Technology, Faculty of Technology and Metal Technology, Lappeenranta. [8] Mahamood, R. M., Akinlabi, E.T., Shukla, M., Pityana, S. (2014). Effect of processing parameters on the properties of laser metal deposited Ti6Al4V using design of experiment. IAENG Transactions on Engineering Sciences, p. 331-339, D0I:10.1201/b16763-37. [9] Mahamood, R.M., Akinlabi, E.T., Shukla, M., Pityana, S. (2013). Material Efficiency of Laser Metal Deposited Ti6Al4V: Effect of Laser Power. Engineering Letters, vol 21, no. 1. [10] Mahamood, R.M., Akinlabi, E.T., Shukla, M., Pityana, S. (2013). Scanning velocity influence on microstructure, microhardness and wear resistance performance of laser deposited Ti6Al4V/ TiC composite. Materials and Design, vol. 50, p. 656-666, D0I:10.1016/j.matdes.2013.03.049. [11] Scott, J., Gupta, N., Wember, C., Newsom, S., Wohlers, T., Caffrey, T. (2013). Additive Manufacturing: Status and Opportunities. Science and Technology Policy Institute, Washington D.C. [12] Toyserkani, E., Khajepour, A. (2006). A mechatronics approach to laser powder deposition process. Mechatronics, vol. 16, no. 10, p. 631-641, D0I:10.1016/j.mechatronics.2006.05.002. [13] Choe, H., Abkowitz, S., Abkowitz, S.M. (2008). Influence of processing on the mechanical properties of Ti-6Al-4V-based composites reinforced with 7.5 mass% TiC and 7.5 mass% W. Materials Transactions, vol. 49, no. 9, p. 2153-2158, D0I:10.2320/matertrans.MER2008049. [14] Xiang, W., Xuliang, M., Xinlin, L., Lihua, D., Mingiia, W. (2012). Effect of boron addition on microstructure and mechanical properties of TiC/Ti6Al4V composites. Materials & Design, vol. 36, p. 41-46, D0l:10.1016/j.matdes.2011.10.040. [15] Dalili, N., Edrisy, A., Farokhzadeh, K., Li, J., Lo, J., Riahi, A.R. (2010). Improving the wear resistance of Ti-6Al-4V/TiC composites through thermal oxidation (TO). Wear, vol. 269, no. 7-8, p. 590-601, D0I:10.1016/j.wear.2010.06.006. [16] Wang, F., Mei, J., Jiang, H., Wu, X. (2007). Microstructure and processing properties of Laser fabrication of Ti6Al4V/ TiC composites using simultaneous powder and wire feed. Materials Science and Engineering. A, vol. 445-446, p. 461466, D0I:10.1016/j.msea.2006.09.093. [17] Rastegari, H.A., Asgari, S., Abbasi, S.M. (2011). Producing Ti-6Al-4V/TiC composite with good ductility by vacuum induction melting furnace and hot rolling process. Materials & Design, vol. 32, no. 10, p. 5010-5014, D0I:10.1016/j. matdes.2011.06.009. [18] Uhlmann, E., Hollan, R., El Mernissi, A. (2006). Hybrid Dry-Ice Blasting Laser Processing: Nd-YAG-Laser-assisted Dry-Ice Blasting for De-Coating. Strojniški vestnik - Journal of Mechanical Engineering, vol. 52, no. 7-8, p. 458-462. [19] Weng, F., Chen, C., Yu, H. (2014). Research status of laser cladding on titanium and its alloys: A review. Materials and Design, vol. 58, p. 412-425, D0I:10.1016/j. matdes.2014.01.077. [20] Grigoriev, S.N., Tarasova, T.V., Gvozdeva1, G.O., Nowotny, S. (2014). Structure Formation of Hypereutectic Al-Si-Alloys Produced by Laser Surface Treatment. Strojniški vestnik -Journal of Mechanical Engineering, vol. 60, no. 6, p. 389-394, D0I:10.5545/sv-jme.2013.1211. [21] Thévenot, F. (1990). Boron Carbide - A comprehensive review: Journal of the European Ceramic Society, vol. 6, no. 4, p. 205225, D0I:10.1016/0955-2219(90)90048-K. [22] Levin, L., Frage, N., Dariel, M.P. (2000). A novel approach for the preparation of B4C-based cermets. International Journal of Refractory Metals and Hard Materials, vol. 18, no. 1-2, p. 131135, D0I:10.1016/S0263-4368(00)00012-3. [23] Zhang, G.J, Ando, M., Yang, J.F., Ohji, T., Kanzaki, S. (2004). Boron carbide and nitride as reactants for in situ synthesis of boride-containing ceramic composites. Journal of the European Ceramic Society, vol. 24, no. 2, p. 171-178, D0I:10.1016/S0955-2219(03)00607-1. [24] Dudina, D. V., Hulbert, D. M., Jiang, D., Unuvar, C., Cytron, S.J., Mukherjee, A. K. (2008). In situ boron carbide-titanium diboride composites prepared by mechanical milling and subsequent Spark Plasma Sintering. Journal of Material Science, vol. 43, no. 10, p. 3569-3576, D0I:10.1007/s10853- 008-2563-8. [25] Roy, S., Sarkar, A., Suwas, S. (2010). Characterization of deformation microstructure in Boron modified Ti-6Al-4V alloy. Materials Science and Engineering: A, vol. 528, no. 1, p. 449458, D0I:10.1016/j.msea.2010.09.026. [26] Arrazola, P.J., Garay, A., Iriarte, L.M., Armendia, M., Marya, S., Le Maitre, F. (2009). Machinability of titanium alloys (Ti6Al4V and Ti555.3). A review. Journal of Materials Processing Technology, vol. 209, no. 5, p. 2223-2230, D0I:10.1016/j. jmatprotec.2008.06.020. [27] Hadâr, A., Nica, M. N., Constantinescu, I. N, Pastramâ, S. D. (2006). The constructive and geometrical optimization of the junctions in structures made from laminated composite materials. Strojniški vestnik - Journal of Mechanical Engineering, vol. 52, no.7-8, p. 546-551. [28] Tian, Y.S., Chen, C.Z., Chen, L.B., Liu, J.H., Lei, T.Q. (2005). Wear properties of alloyed layers produced by laser surface alloying of pure titanium with B4C and T{i} mixed powders. Journal of Materials Science, vol. 40 p. 4387-4390, DOI:10.1007/s10853-005-0736-2. [29] Tian, Y.S., Zhang, Q.Y., Wang, D.Y. (2009). Study on the microstructures and properties of the boride layers laser fabricated on Ti-6Al-4V alloy. Journal of Materials Processing Technology, vol. 209, no. 6, p. 2887-2891, DOI:10.1016/j. jmatprotec.2008.06.043. [30] Struers application Note on titanium. (2015). From http:// www.struers.com/resources/elements/12/104827/ Application_Note_Titanium_English.pdf, accessed on 201515-04. [31] ASTM E384-11e1. (2011). Standard Test Method for Knoop and Vickers Hardness of Materials. ASTM International West Conshohocken, DOI:10.1520/E0384-11E01. [32] ASTM Standard G133-05. (2005). Standard Test Method for Linearly Reciprocating Ball-on-Flat Sliding Wear, ASTM International. West Conshohocken, DOI:10.1520/G0133-02. [33] Shen, X., Cao, L., Li, R. (2010). Numerical simulation of sliding wear based on Archard's model. International Conference on Mechanic Automation and Control Engineering, DOI:10.1109/ MACE.2010.5535855. [34] Ravnikar, D., Mrvar, P., Medved, J., Grum, J. (2013). Microstructural Analysis of Laser Coated Ceramic Components TiB2 and TiC on Aluminium Alloy EN AW-6082-T651. Strojniški vestnik - Journal of Mechanical Engineering, vol. 59, no. 5, p. 281-290, DOI:10.5545/sv-jme.2012.904. Strojniški vestnik - Journal of Mechanical Engineering 63(2017)6, 374-382 © 2017 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2016.4184 Original Scientific Paper Control-Oriented Modelling with Experimental Verification and Design of the Appropriate Gains of a PI Speed Ratio Controller of Chain CVTs Ahmet Yildiz* - Osman Kopmaz Uludag University, Engineering Faculty, Mechanical Engineering Department, Turkey The continuously variable transmission (CVT) system represents one of the best solutions to minimize the fuel consumption for vehicles that are driven by an internal combustion engine or an electrical motor Hence, the theoretical analysis of a chain CVT is crucial in the optimization of the design and control strategy. This paper is concerned with the transient dynamics of a chain CVT with a speed ratio controller adopted from the Carbone-Mangialardi-Mantriota model and its experimental verification. To this end, a theoretical model is developed with a PI speed ratio controller considering the momentary solution of a first order differential equation governing the shifting speed, which represents the whole dynamics of the CVT. To verify the developed model, the experiments are carried out on a chain CVT test rig. It is observed that the numerical and experimental results are in good agreement, which implies that the developed model embedded in a speed ratio controller is appropriate to predict the shifting dynamics of the chain CVT. Afterwards, the developed model is used to design the appropriate PI gains by numerical experiments in order to obtain the same slope of the time response of speed ratio for different input angular velocities. Therefore, the same shifting speed can be secured for the different working conditions in this manner. This paper contains very important results for manufacturers about the control parameter and gain effect on the shifting dynamics of chain CVT, and the developed model can be used for different control algorithms in automotive applications. Keywords: chain CVT, speed ratio control, PID, shifting dynamics Highlights • Obtaining a lookup table of the clamping force ratio under the steady-state condition from the theoretical model with experimental verification. • Implementing this lookup table in the first order differential equation and feedbacking the speed ratio with a PI controller. • Comparison of the theoretical and experimental results of the time response of chain CVT with a controller. • Design a PI gain schedule to ensure the same shifting speed by numerical experiment method. 0 INTRODUCTION The environmental problems caused by harmful gas emissions of the vehicles remain significant issues requiring further investigation. To minimize this pollution, research on the efficiency of the vehicle in terms of the engine and transmission systems is increasingly common. The continuously variable transmission (CVT), shown in Fig. 1, is one of the best solutions to reduce the fuel consumption for vehicles that are driven by internal combustion engines [1], electric motors [2], and hybrid systems [3]. This is a result of the CVT systems providing smoothness in changing of speed ratio, resulting in reduced fuel consumption for the vehicles. Moreover, CVT systems are significant and advantageous alternatives in the machinery industry when continuous speed variation is required. Mangialardi and Mantriota [4] analyse the dynamic behaviour of a wind power system equipped with automatically regulated, continuously variable transmission. Yildiz and Kopmaz [5] propose a novel mechanical press concept that consists of a half- toroidal continuously variable transmission (CVT) to obtain more speed flexibility than the classical one. Yildiz et al. [6] analyse the dynamic behaviour of a four-bar mechanism coupled with a CVT system. Thus, the speed variability which is crucial for obtaining different manufacturing and transportation processes can be provided by mechanically. To obtain better transmission performance, the dynamics of the transmission must be fully understood. Determining the dynamics of CVT is crucial to obtaining and controlling the desired the speed ratio accurately. In this type of variator, the positions of the pulleys are varied by a hydraulic system, which results in a change in the contact radii of chain at the drive and driven sides; thus, the desired value of the transmission ratio is obtained. All these dynamic phenomena have been studied for decades. In particular, the theoretical analyses for both steady-state and transient cases have been studied in depth by Carbone, Mangialardi and Mantriota [7] and Yildiz et al. [8]. A theoretical model of the shifting dynamics of a chain CVT during rapid ratio speed changes is developed by Carbone et al. [7]. The clamping force ratio as a function of torque load and traction coefficient is derived in detail. Fig. 1. Chain CVT system Carbone et al. [9] demonstrate that the pulley bending has an important role in the dynamics of the variator. The sliding angle that is one of the main parameters in equilibrium equations varies if the pulley deformation is considered. Therefore, it is proved that neglecting the pulley bending causes errors. It is also shown that the shifting speed of the variator increases as the stiffness of the pulley decreases. The dynamic behaviour of the CVT is investigated at a steady state and in transient conditions by Carbone et al. [10]. However, in their study, the shifting dynamics are examined under no-load conditions. To verify the proposed theoretical model, experimental investigations are performed. Furthermore, the slip of the variator that is one of the crucial issues of the mechanical behaviour is analysed and controlled for efficiency optimization in [11]. Controlling the speed ratio is one of the crucial issues of the CVT systems, and there are several methods to do so. The most common technique is the proportional-integrate-derivative (PID) controller. The controller design with PI(D) and their applications are examined in [12] and [13]. Moreover, PI and PID controllers tuning is investigated by Precup and Preitl in [14] for integral-type servo-systems to ensure strong stability and controller robustness. Although extensive works have been performed on the steady-state performance and the mechanical behaviours of the variator, such as slip predictions, traction capacities, the efficiency [15] to [17], only a few studies focus on the shifting dynamics of the CVT [18] and [19] and the speed ratio control in real time. Furthermore, most researchers use semi-empirical formulations in their model that can cause some errors in different applications. This paper differs from the previous works in that the shifting dynamics of the CVT with a PI controller is fully modelled based on the Carbone-Mangialardi-Mantriota model, and its experimental verification is carried out. From the theoretical model, a lookup table of the clamping force ratios as a function of speed ratio and traction coefficient is obtained for fast implementation in the control algorithm. After that, a PI feedback control algorithm is designed considering the momentary solution of a first order differential equation governing the shifting speed, which represents the whole dynamics of the CVT. To verify the developed model, several experiments are performed by setting the clamping forces of the drive pulleys with a feedback PI controller in a real time Labview environment. The experimental results show good agreement with the theoretical ones obtained from the control-oriented model based on the Carbone-Mangialardi-Mantriota approach. After verifying the model, it is used to obtain the appropriate PI gains via numerical experiments. A PI gain schedule is proposed for the different angular velocity of the drive pulley, which affects the shifting dynamics sharply. This model may constitute the basis of the different CVT control strategies, in particular, the feedforward control algorithm. 1 MODELLING THE CHAIN CVT In this part, the theoretical model of the speed ratio control based on the Carbone-Mangialardi-Mantriota model is presented. First, the main equations of this model are presented. After that, a control algorithm for the speed ratio control of the CVT is proposed. As presented in [7], the basic kinematic quantities and their relations are presented as follows. The mass conservation law yields the following relation as the longitudinal elongation of the chain is neglected: v +- dvg de = 0, (1) where vr and ve are radial and tangential sliding velocity of the chain respectively. The sliding angle ^ is defined as; tan^ = —-. Vr The radial sliding velocity is calculated as: v = R + AmR (1 + cos2 ßo ) . sin2ß0 sin(0-0c), (2) (3) where R is the pitch radius of the chain, fi0 is the non-deformed half-groove angle of the pulley, Qc is the centre of wedge expansion, m is the angular velocity of the pulley, and A is the amplitude of the deformed groove angle given in [15]. The deformed half groove angle is determined as a function of local angular coordinate Q, which considers a uniform deformation r of the pulley, given in [15]: 1 d p-p0 =r + 0.5 A cos(Q-Qc ). (4) The centre of wedge expansion can be estimated as follows: tan Qc = - u. J p(d) sin ddd _0_ a J p(d)cosddd (5) where a is the contact arc extension and p is the force per unit longitudinal chain length (linear pressure). aj Fig. 2. The kinematic representation of a) chain CVT and b) the forces acting on the chain All forces acting on the chain are represented in Fig. 2. The equilibrium equations in the radial and tangential direction yield the following equations: F-gR2o>2 d0 U cos Ps sin^ (F -oRO>2 ) = P = - (sinP0 — jucosPs cosy) F-gR2G)2 (7) 2R [sin P0 - p cos Ps cos^ ] where F is the tension of the chain, and a is the mass per unit length of the chain. The angle satisfies the relation: tan PS = tan Pcosy. (8) In the theoretical model, a friction coefficient / = 0.09 is implemented according to Coulomb friction. Once the local pressure and tension distributions are known, the axial clamping force S and torque T acting on each pulley are calculated as below: a S = J (cos P+,u sin PS) pRdd, (9) 0 T = (F -F1)R, (10) where F1 and F2 are the tensile belt forces at the entry and exit points of the pulley groove. In Eq. (9), the extension arc of the chain a can be calculated with the following relations: smffl = R2 - R1 d aj = n- 20, a2 = n + 20, (11) (12) (13) where d is the distance of the pulleys' centres. The Carbone-Mangialardi-Mantriota model shows that the whole dynamics of the CVT can be represented as a first order non-linear differential equation as follows: T =A DRk^DRg ln I (14) v eq J where k = (1 + cos2 P0) / sin (2P0), % = SDR / SDN and is the clamping force ratio at the steady state condition. As can be seen from Eq. (14), the time derivative of the speed ratio is a function of the input angular velocity, the elastic displacement of the input pulley, the non-deformed value of the groove angle, the actual value of the transmission ratio through the function 0 g(r) and the logarithmic ratio of £ and £eq. The term g(z) is as follows: « T)= TC 1 -T \n + 2 arcsin r (i l p| 1 \ 1 _ V T )_ \n- 2 arcsin r (i l p| i 1 V T ) J (15) where p(r)=R/ d and c ~ 0.81. Since the SDR/SDN is a function of the traction coefficient and speed ratio, it is necessary to yield the traction coefficient futr\ Vtr = cos ß0 T 2 RS ' (16) In order to calculate the clamping forces ratio, the numerical steps are followed as shown in Fig. 3. A clamping forces ratio matrix is obtained by the same way for use in the control algorithm. This matrix is a function of the traction coefficient and speed ratio. Fig. 3. Numerical flowchart for calculating the necessary clamping force ratio 2 SPEED RATIO CONTROLLER DESIGN In order to control the speed ratio of the chain CVT, a PI controller is used, and the transient dynamics of the variator are investigated. In this controller, the initial condition of the angular velocity and the torque load of the drive pulley are given, and then the initial value of the clamping force ratio calculated by the theoretical model is implemented. For a constant clamping force of the driven pulley, the speed ratio is fed backed. The program controls the speed ratio by changing the clamping force of the drive pulley while the new conditions calculated in CVT box spontaneously until it reaches the desired value of the speed ratio. Since the proportional pressure valves are much faster than the variator in terms of the time constant, the dynamics of the hydraulic system is neglected. Ec^Q- Reference speed ratio PID(s) PID controller Kv Initial wDRl condition of SDR/SDN Fig. 4. a) PI controller for the speed ratio and b) lookup tables The error value between the desired set-point value Tref and measured process variable of speed ratio is: e(t) = rref -T(t). (17) The control variable, which is the voltage of the servo valve: 1 1 u(t) = Kpe(t) + - J e(t)dt (18) where Kp and T are the proportional gain and integral time gain, respectively. i 0 The lookup table of the clamping force ratio {(SDR/ SDN)eq = shown in Fig. 4 is obtained by the theoretical model, following the flowchart given in Fig. 3 for different speed ratios and traction coefficients in the range of 0.5 to 2 and 0.05 to 0.1, respectively. Another lookup table (m) is obtained from the relationship m(r) = kg(r) where k and g given in Eq. (14) and (15). The pitch radius of the chain {R(t)} is also given as a lookup table, which is calculated from Eqs. (11), to (13) and (24). The amplitude of the deformed groove angle {A(t, SDR)} is obtained from finite element FE methods [15]. The term Kv shown in Fig. 4 is the gain coefficient of the valve and hydraulic system, which equals 13180. The parameters of the PI controller are calculated with numerical experiments to ensure the same time response of the system. A gain scheduling approach is proposed, by which the gains of a PI controller are calculated through numerical experiments. These results are given in Section 4. The parameters of PI controller is selected in the experiments as Kp = 1 and T = 0.1 for avoiding aggressive speed ratio changing that may damage the test rig equipment. 3 CHAIN CVT TEST RIG AND EXPERIMENTAL SETUPS The test rig and its setups are described in this part. The measurements are carried out on a chain CVT mounted on the power loop test rig in Polytechnic of Bari. The chain CVT test bench is illustrated in Fig. 5, where EM, EB, HE, LE, and TS represent electrical motor, electromagnetic brake, heat exchanger, linear encoder and torque sensor, respectively. As shown in Fig. 5, a three-phase asynchronous four-pole, 30 kW Siemens motor is located to ensure the input power of CVT. The angular velocity of the motor is controlled by a Berges inverter. To provide the resistant torque, an electromagnetic brake connected to the output pulley of CVT by a belt transmission is used. The pressure, torque and displacement sensors are located on the input and output pulley shafts. The characteristic data of the chain used for experiments are as follows: chain length L = 649.49 mm, chain width b = 24 mm, the distance of the input and output pulleys' centres d= 155 mm, the sheave angle f}0 = 11 deg. The areas of the input and secondary pulleys are ADR = 0.01979 m2, Adn = 0.009719 m2 respectively. The clamping force of drive side is calculated as follows: SDR PDRADR + mDRfc,DR ' where the coefficient fcDR is: fc,DR o )(rDRl rDR 2)' (19) (20) where poi/ = 870 kg/m3 at 20 0C and rDR1 = 0.0825 m and rDR2 = 0.0225 m. The driven clamping force is calculated as follows: SDN = PDNADN + ®DNfc,DN + Fspring ' (21) where Fspring is the spring force on the driven pulley and the coefficient fc DN is: fc ,DN =n(p )(2rDN r r 2 2 _ 4 rDN irDN 3 rDN 2 ■r4DN 3). (22) where the quantities rDN1 = 0.06 m, rDN2 = 0.0225 m and rDN3 = 0.0321 m. The spring force is calculated using the following equation: F = FS0 + 2ks (rDNlow - rDN)tan(^), (23) Fig. 5. Chain CVT test bench; a) general view, and b) schematic layout where FS0 = 536 N, the stiffness of the spring ks = 20 N/m, the pitch radius of the driven pulley at minimum speed ratio rDNjow = 0.0741 m. The axial position of the input pulley is measured by a linear encoder as a function of time in the experiments, as shown in Fig. 6. After that, the speed ratio of variator t is determined. In order to examine the geometric speed ratio of variator t, the pitch radius Rdr of the drive pulley is calculated with to following equation: L = 2d cos 0 + R1a1 + R2a2. (24) The following equation is used to calculate the RDN: RDN RDN ,max 2 tan ' (25) where the RDN^ax = 0.0741 m. Fig. 6. Representation of the measurement of the axial position of a moveable input pulley by a linear encoder A control program is developed to perform the measurements and save all data in real time in a Labview environment. The measured data are the speed ratio, the pressures of the input and output pulleys, the resistant torque and the angular velocity. In the experiments, firstly, the angular velocity of the driven pulley is fixed, and then a torque load is implemented. After that, the pressure of the drive side is controlled and so the desired speed ratio is obtained due to the PI controller. In order to investigate the transient dynamics of CVT, a step of the set point of the speed ratio is implemented. 4 RESULTS AND DISCUSSION In this section, the theoretical and experimental results and their comparison are demonstrated. First, the steady state results of the clamping forces ratios are given in Fig. 7 for the speed ratio equal to 1. It can be seen from this figure that the theoretical model is valid in steady state conditions, since it has a very good match with experimental results. Similar good matches are obtained for different speed ratios, which are not presented here. These steady state data are recorded to be used as a lookup table during the simulation. The developed model for the speed ratio control is verified using several experiments for the different initial conditions of the speed ratio, as indicated in Fig. 8. 1.5 1.4 Cr Q 1.3 CO ce Q 1.2 CO 1.1 -Experiment -Theory 0.02 0.04 0.06 0.08 tr,DN Fig. 7. Theoretical and experimental results of the clamping force ratio as a function of traction coefficient at SDN = 10 kN and t=1 For the further comparison, the maximum and minimum speed ratios are implemented as an initial condition of the experiments and simulations shown in Fig. 9. It is observed from all these figures that the theoretical and experimental results are in good agreement. Thus, the developed model is verified, and it can be used in different controllers in automotive applications. 100 150 Time [s] Fig. 8. Theoretical and experimental results of the transient response of the speed ratio for different initial condition of speed ratio at Kp = 1, T = 0.1, rndr=500 rpm, TdR = 50 Nm and Sdn=20 kN At the jumping point, there are some small errors due to the neglected factors, such as hydraulic system dynamics and the friction in the bearings of the test bench; however, these are negligible when the overall consistency of the numerical and experimental results are considered. Since the CVT has a non-linear dynamic characteristic, the shifting response will be different for the different conditions especially for the different input speeds. In other words, the input angular speed is the most dominant term that affects the shifting speed linearly. Thus, the PI coefficients handled in Figs. 8 and 9 need to be tuned for the different angular velocities. 20 Time [s] Fig. 9. Theoretical and experimental results of the transient response of the speed ratio for testing maximum and minimum range of speed ratio at Kp=1, T = 0.1, mDR=500 rpm, TDR=50 Nm and SDN=20 kN 10 Time [s] Fig. 10. Time response of the speed ratio for different angular velocity at Kp=1, T = 0.3, TDR=50 Nm and SDN=20 kN Fig. 10 indicates the time responses of the speed ratios for the different angular velocities of the input pulley with the same PI coefficients. This figure is evidence that PI gains of the speed ratio controller of the chain CVT have to be set for the different angular velocities. To ensure the regular shifting response such as the same slope of the time response, the gains are scheduled by numerical experiments for different input speeds. First, the T is kept constant and Kp is changed, and the settling time of speed ratio is obtained in case of no over-jumping conditions. The same procedure is followed for different integral time gains. The obtained data are recorded in a matrix considering the settling time of the speed ratio in case of the absence of over-jumping conditions. Fig. 11 indicates the settling times of the speed ratios for three different input angular velocities. In Fig. 11, the minimum points of the surfaces show the best PI coefficients in terms of the settling time with no over jump. In this manner, a PI gains schedule is obtained as given in Table 1. That the difference between the desired and instant speed ratios be less than 0.001 is used as the criterion to determine the settling time seen in Fig. 11. Fig. 11. Settling time of the speed ratio for the different input angular velocities Table 1. PI Gains for different initial angular speed of the CVT Input angular velocity (mDR) Kp T 500 rpm 4.3 2.2 750 rpm 3.25 2.35 1000 rpm 2.5 2.35 1250 rpm 2.075 2.35 1500 rpm 1.75 2.35 1750 rpm 1.575 2.35 2000 rpm 1.35 2.35 The time responses of the speed ratio are simulated again for the visibility using the data given in Table 1, and the results are given in Fig. 12. It can be seen from this figure that all time responses have the same slope and there is no overshooting. It should be noted that the time constant of the speed ratio is selected in the beginning and the gains are tuned depending on the demand. -%R = 2000 rpm ■ / -%R = 1750 rpm / — (aDR = 1500 rpm / %R = 1250 rpm / — %R = 1000 rpm / ^DR ~ 750 rpm / -%R = 500 rpm 10 0.5 1 1.5 Time [s] Fig. 12. The transient response of the speed ratio for the different input speed with the PI coefficients given in Table 1 (Tdr=50 Nm and SDN=20 kN) The theoretical and experimental results show that the control-oriented modelling proposed in this paper is highly accurate, thus; this simplified formula can be used for different control design as a tool with specified lookup tables. 5 CONCLUSIONS In this study, the shifting dynamics of the chain CVT is investigated to verify the theoretical model in a control application. A speed ratio controller is designed to provide a relative simple and easy-to-implement formulation as a look-up table and various experiments are performed to validate this model. The theoretical and experimental results indicate that the control-oriented modelling of transient dynamics of the variator is highly accurate. The verified model is used to obtain appropriate PI gains for having regular shifting speeds of the CVT under the different angular velocities via numerical experiments. Finally, a gain schedule is proposed for possible input speeds. These results are vital for engineers to design different control algorithms in automotive applications. 6 ACKNOWLEDGEMENTS Ahmet YILDIZ would like to thank the TUBITAK for the scholarship granted. Special thanks are due to Prof. Dr. Giuseppe Carbone for suggestions and allowance the first author to perform the experiments in the test bench in Polytechnic of Bari, Italy. 7 REFERENCES [1] Brace, C., Deacon, M., Vaughan, N.D., Horrocks R.W., Burrows C.R. (1999). The compromise in reducing exhaust emissions and fuel consumption from a diesel CVT Powertrain over typical usage cycles. Proceedings of CVT Congress, Eindhoven, p. 27-33. [2] Bottiglione, F., De Pinto, S., Mantriota, G., Sorniotti, A. (2014). Energy consumption of a battery electric vehicle with infinitely variable transmission. Energies, vol. 7, no. 12, p. 8317-8337, D0l:10.3390/en7128317. [3] Bottiglione, F., Contursi, T., Gentile, A., Mantriota, G. (2014). The fuel economy of hybrid buses: the role of ancillaries in real urban driving. Energies, vol. 7, no. 7, p. 4202-4220, D0I:10.3390/en7074202. [4] Mangialardi L., Mantriota G. (1996). Dynamic behavior of wind power systems equipped with automatically regulated continuously variable transmission. Renewable Energy, vol. 7, no 2, p. 185-203, D0I:10.1016/0960-1481(95)00125-5. [5] Yildiz, A., Kopmaz, O. (2015). Dynamic analysis of a mechanical press equipped with a half-toroidal continuously variable transmission. International Journal of Materials and Product Technology, vol. 50, no. 1. p. 22-36, D0I:10.1504/ IJMPT.2015.066864. [6] Yildiz, A., Kopmaz, O., Telli, C.S. (2015). Dynamic modeling and analysis of a four-bar mechanism coupled with a CVT for obtaining variable input speeds. Journal of Mechanical Science and Technology, vol. 29, no. 3, p. 1001-1006, D0I:10.1007/s12206-015-0214-y. [7] Carbone, G., Mangialardi, L., Mantriota, G. (2001). Theoretical model of metal V-Belt drives during rapid ratio changing. Journal of Mechanical Design, vol. 123, no. 1, p. 111-117, D0I:10.1115/1.1345521. [8] Yildiz, A., Piccininni, A., Bottiglione, F., Carbone, G. (2016). Modeling chain continuously variable transmission for direct implementation in transmission control. Mechanism and Machine Theory, vol. 105, p. 428-440, D0I:10.1016/j. mechmachtheory.2016.07.015. [9] Carbone, G., Mangialardi, L., Mantriota, G. (2005). The influence of pulley deformations on the shifting mechanisms of MVB-CVT. ASME Journal of Mechanical Design, vol. 127, no. 1, p. 103-113, D0I:10.1115/1.1825443. [10] Carbone, G., Mangialardi, L., Bonsen, B., Tursi, C., Veenhuizen, P.A. (2007). CVT dynamics: theory and experiments. Mechanism and Machine Theory, vol. 42, no. 4, p. 409-428, D0I:10.1016/j.mechmachtheory.2006.04.012. [11] Bonsen, B., Klaassen, T.W.G.L., Pulles, R.J., Simons, S.W.H., Steinbuch, M., Veenhuizen, P.A. (2005). Performance optimization of the push-belt CVT by variator slip control. International Journal of Vehicle Design, vol. 39, no. 3, p. 232256, D0I:10.1504/IJVD.2005.008473. [12] Besangon-Voda, A. (1998). Iterative auto-calibration of digital controllers: methodology and applications. Control Engineering Practice, vol. 6, no. 3, p. 345-358, D0I:10.1016/ S0967-0661(98)00005-7. [13] Jin, B.Q., Liu, Q. (2014). IMC-PID design based on model matching approach and closed-loop shaping. ISA Transactions, vol. 53, no. 2, p. 462-473, D0I:10.1016/j.isatra.2013.11.005. [14] Precup, R.-E., Preitl, S. (2006). PI and PID controllers tuning for integral-type servo systems to ensure robust stability and controller robustness. Electrical Engineering, vol. 88, no. 2, p. 149-156, D0I:10.1007/s00202-004-0269-8. [15] Carbone, G., De Novellis, L., Commissaris, G., Steinbuch, M. (2010). An enhanced CMM model for the prediction of steady state performance in CVT chain drives. ASME Journal of Mechanical Design, vol. 132, no. 2, p. 1-8, D0l:10.1115/1.4000833. [16] Kong, L., Parker, R.G. (2008). Steady mechanics of layered, multi-band belt drives used in continuously variable transmissions (CVT). Mechanism and Machine Theory, vol. 43, no. 1, p. 171-185, D0I:10.1016/j. mechmachtheory.2007.02.003. [17] Tenberge, P. (2004). Efficiency of chain-CVTs at constant and variable ratio: A new mathematical model for a very fast calculation of chain forces, clamping forces, clamping ratio, slip and efficiency. International Continuously Variable and Hybrid Transmission Congress, paper no. 04CVT-35, UC Davis. [18] Ide, T., Uchiyama, H., Kataoka, R. (1999). Experimental investigation on shift speed characteristics of a metal V-belt CVT. JSAE paper no.96363. [19] Srivastava, N., Haque, I. (2009). A review on belt and chain continuously variable transmissions (CVT): dynamics and control. Mechanism and Machine Theory, vol. 44, no. 1, p. 1941, D0I:10.1016/j.mechmachtheory.2008.06.007. Strojniški vestnik - Journal of Mechanical Engineering 63(2017)6, 383-393 © 2017 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2017.4449 Original Scientific Paper Mixed Convection Heat Transfer in a Square Lid-Driven Cavity Filled with Al2O3-Water Nanofluid Yazan Taamneh1* - Kahled Bataineh2 1 Jordan University of Science and Technology, Department of Aeronautical Engineering, Jordan 2 Jordan University of Science and Technology, Department of Mechanical Engineering, Jordan This work is focused on the steady laminar mixed convection flow in a lid-driven square cavity filled with Al2O3 water-nanofluid using computational fluid dynamics. The top lid of the cavity was kept at a higher temperature in comparison with the bottom wall and moving at a constant speed while the left and right walls were kept insulated. Simulations were performed using Fluent ver. 6.3 where the water based nanofluid was considered as a single phase. A parametric study was conducted, and the effects of the Richardson number (0.1 < Ri < 10), the volume fraction of the nanoparticle (0 < < < 0.04) on the fluid flow, and heat transfer inside the cavity were investigated. It was found that when (1 < Ri < 10) the average Nusselt number increases accordingly for the whole range of solid volume fraction of the nanofluid. The simulation results showed that the presence of nanoparticles in the base fluid increases the heat transfer rate. The variations of isotherm and streamline patterns inside the cavity with different volume fractions of nanoparticle and Richardson number were discussed and demonstrated. Keywords: mixed convection, nanofluid, square cavity, CFD, lid-driven Highlights • The research focused on steady laminar mixed convection flow in a lid-driven cavity. • The cavity was filled with Al2O3 water nanofluid. • The effect of the volume fraction of the nanoparticles was investigated. • The heat transfer rate was enhanced by increasing the solid volume fraction. 0 INTRODUCTION The combination of the mixed convection flow and the heat transfer analysis in a lid-driven cavity has recently received significant attention by many researchers due to its significance in many practical engineering and industrial applications. The flow and heat transfer due to mixed (free and forced) convection are a combination of the flow induced by the shear forces resulting from the motion of the upper lid and the bouncy forces due to the temperature gradient. The practical flow situation is considered to be a complex phenomenon because of the interaction and coupling of these effects. In addition, there are various combinations of imposed temperature gradients and many cavity configurations of a lid-driven cavity encountered in several engineering applications. Cha and Jaluria [1] studied in detail the impact of bounciness, velocities, and temperature fields on the heat transfer enhancement inside storage regions. Ghia et al. [2] investigated the driven flow in a square cavity at a high-Re fine-mesh flow. The solid particles are used as an additive suspended in the base fluid to enhance the heat transfer in a variety of applications, such as micro-electromechanical systems (MEMS), oil extraction and automobiles [3], electronic cooling, thermo-hydraulics of nuclear reactors, and solar ponds [4]. Kumar et al. [5] developed bactericidal coatings by using synthesized metal-nanoparticle (MNP)-embedded paint. They showed that the surfaces coated with silver-nanoparticle paint have excellent antimicrobial properties. Choi [6] proposed a new class of engineering fluids with superior nanofluid thermal conductivity by adding nanoparticles. Vajjha and Das [7] studied how nanofluid thermophysical properties varied with temperature and concentration. They analysed their effects on the heat transfer coefficient, the pumping power, and friction factor. They asserted that enhancing the heat transfer rate by adding solid materials reduces operation costs. Numerous experimental and numerical investigations have been conducted on the heat transfer enhancement using nanofluid. Kaheld and Vafai [8] investigated the heat transfer enhancement inside channels by commanding thermal dispersion influences. They found that the distribution of dispersive elements increases the heat transfer. Torrance et al. [9] examined early the fluid motion inside different aspect ratio cavities generated by a moving wall under natural and mixed convection. They revealed that when Grashof number increases especially in cavities with high aspect ratio, the buoyancy effect would be significant. The buoyancy-driven heat transfer enhancement of nanofluid in a two-dimensional enclosure has been analysed by Khanafer et al. [10]. Furthermore, various studies of *Corr. Author's Address: Jordan Univesity of Science and Technology, P. O. Box 3030, 22110 Irbid, Jordan, ymtaamneh@just.edu.jo 383 the thermophysical properties of nanofluid and energy transport were carried out. Jung et al. [11], Kang et al. [12], and Lee et al. [13] measured the effective thermal conductivity for different types of nanofluid using the transient hot wire method. Wang and Mujumdar [14] outlined the recent research on fluid flow and heat transfer characteristics of nanofluid in forced and free convection flows. Pak and Cho [15] showed experimentally that the convection heat transfer coefficient for Al2O3-water and TiO2-water nanofluid was enhanced by increasing the volume fraction as well as the average velocity. Xuan and Li [16] and [17] conducted an experimental study to investigate the effect of copper nanoparticles on the thermal conductivity of the nanofluid. Their results demonstrated that the suspended nanoparticles clearly enhanced the thermal conductivity of the base liquid. Furthermore, increasing the volume fraction of the nanoparticles leads to an increase of the thermal conductivity of the nanofluid. Abu-Nada and Oztop [18] carried out a numerical investigation of a mixed convection in an inclined square cavity of an aluminium-based nanofluid. They concluded that convection heat transfer was significantly enhanced by the existence of different types of nanoparticles. Tiwari and Das [19] confirmed that increasing the volume fraction of nanoparticles significantly enhances the convection heat transfer coefficient. They inspected the flow and heat transfer in a square enclosure utilizing the finite volume method by considering that the top and the bottom were insulated and the side walls remain at different constant temperatures. It is evident from the literature that the thermal conductivity of nanofluid highly depends on the volume fraction, the thermal conductivity (of both the base fluid and the nanoparticle material), the surface area, and the shape of the nanoparticles suspended in the base liquid. There are few studies in the literature that investigated the effect of a nanofluid in forced convection. In contrast, a large number of studies were devoted to the effect of nanoparticle in the natural convection. The free convection in a rectangular cavity was numerically studied by Jang and Choi [19] and Jou and Tzeng [21]. They presented their numerical results of the heat transfer enhancement of nanofluid in a two-dimensional cavity. The influence of gravity on sedimentation and the agglomeration of nanoparticle on a natural convection heat transfer was investigated by Jafari et al. [22]. They proposed a single-phase approach and a mixture model. Ismail et al. [23] studied the buoyancy forces as the driving heat transfer of nanofluid, using FLUENT. The effect of particle sizes and shapes on the heat transfer characteristics and the pressure losses were studied by Merilainen [24]. They determined that the average convective heat transfer coefficients of nanofluid were improved up to 40 % compared to clear fluids. The mixed convection flow and the heat transfer analyses in a lid-driven inclined enclosure filled with a nanofluid were conducted by Iwatsu et al. [25]. The flow and the heat transfer in a square cavity with insulated top and bottom walls, and differentially-heated moving sidewalls were also investigated using the finite volume approach by Mansour and Ahmed [26]. They investigated the effects of the Richardson number and the volume fraction of the nanoparticles on the heat transfer. It was observed that when the Richardson number equals unity, the average Nusselt number increased substantially with the increase in the volume fraction of the nanoparticles. The mixed convection flow in a lid-driven enclosure filled with a fluid-saturated porous medium was investigated by Abu-Nada and Chamakha [27]. They reported the effects of the Darcy and Richardson numbers on the flow and the heat transfer characteristics. A numerical study of laminar mixed convection in shallow driven cavities with a hot moving lid on top and cooled from the bottom has been carried out by Sharif [28]. It was observed that the average convection heat transfer increased slightly with the cavity inclination angle for the forced convection while it increased more significantly under natural convection. Oztop et al. [29] investigated the heat transfer and the fluid flow due to buoyancy forces in a partially heated enclosure using different nanofluids. Calculations were performed for different Rayleigh numbers, heights of heater, locations of heater, aspect ratios, and volume fractions of nanoparticles. They found that the heat transfer increases with increases of the Rayleigh number, the height of heater and the volume fraction of nanoparticles. Furthermore, the heat transfer enhancement was found to be dependent on the type of nanofluid and more pronounced at lower aspect ratios. Polidori et al. [30] utilized the integral formalism approach to investigate the natural convection heat transfer of a Newtonian nanofluid in a laminar external boundary-layer. They dealt with y-Al2O3/ water nanofluid whose Newtonian behaviour was experimentally confirmed for particle volume fractions less than 0.04. They concluded that generalized conclusions about the heat transfer enhancement with the use of a nanofluid needs to be carefully drawn. In addition, they found that the natural convection heat transfer is not solely characterized by the nanofluid's effective thermal conductivity. Moreover, they concluded that the sensitivity to the used viscosity model cannot be ignored as it plays a key role in the heat transfer behaviour. Putra et al. [31] executed an experimental work on the natural convection of Al2O3 and CuO-water nanofluid inside a horizontal cylinder heated and cooled from both ends, respectively. They found that the presence of nanoparticles spoiled the natural convective heat transfer systematically. They also observed a systematic degradation of the natural convective heat transfer with increased concentrations of particles. Haddad et al. [32] and [33] studied numerically the natural convection heat transfer and fluid flow of CuO-water nanofluid in an open cavity heated from the bottom. They considered in their calculation the variation in the viscosity and the thermal conductivity. They found that the heat transfer decreases with increasing of solid volume fraction as a result of increasing the viscosity of the nanofluid. Although many studies are found in the literature investigating the heat transfer characteristics of nanofluids, there is a lack of studies devoted to combined mixed convection flow and heat transfer. Moreover, due to nanofluid's superior thermophysical properties, utilizing nanofluids in industrial applications has recently witnessed rapid growth. However, there is still no accurate understanding of the effect of using nanoparticles in the combined convection heat transfer. Therefore, the main objective of this study is to investigate numerically the combined convection flow and heat transfer in a square lid-driven cavity utilizing nanofluid. General correlations for the effective thermal conductivity, the viscosity, and the thermal expansion coefficient of a nanofluid were adopted in this study in terms of the volume fraction, the particle diameter, the temperature, and the base fluid physical properties. These theoretical models were employed utilizing a user-defined function (UDF) in Fluent. move at a constant speed Ux. The physical model considered in this study is shown in Fig. 1. The shape and the size of particles are assumed to be uniform with diameter equal to 100 nm. The fluid properties of the nanofluid vary when nanoparticles are suspended. In this study, the nanofluid is treated as a single phase. The thermophysical properties of the nanoparticles and the fluid phase at T = 300 K are presented in Table 1. Mixed convection flow and heat transfer in a square lid driven cavity are simulated by using the mass, momentum, and energy conservation equations. The governing non-linear partial differential equations can be written as follows: u dU+v U dX dY dU dV n -+ — = 0, dX dY dP 1 Heff dX ReVfPf U dV+V dV=-
Tc. The top wall of the cavity is allowed to
where L is the reference length, Ux is the reference velocity, and v is the kinematic viscosity. The Reynold's number Re is the ratio of inertial to viscous forces, which influences the fluid flow features within the cavity. The ratio of Gr/Re2 is the convection parameter and is called the Richardson number Ri. It measures the relative strength of the natural convection and the forced convection for the present problem.
+
Table 1. Thermophysical properties of the base fluid and aluminium [26]
Physical properties Water AI2O3
Cp [J/(kg K)] 4179 765
k [W/(m2-K)] 0.613 25
p [kg/m3] 997.1 3970
n [Pa-s] 0.000891 -
P [1/K] 0.00021 0.0000017
1.1 Boundary Conditions
A schematic of the configuration analysed in this study is shown in Fig. 1. The appropriate dimensionless boundary conditions for the present study are as follows:
da
— = 0, U = V = 0, X = 0, dY
da
— = 0, U = V = 0, X = 1, dY
0 = 0, U = V = 0, Y = 0, 0 = 1, U = 1, V = 0, Y = 1.
(7)
L= 1 m"
Fig. 1. Schematic of the square lid-driven cavity
coefficient of nanofluid were developed in terms of the solid volume fraction, the particle diameter, the temperature, and the base fluid physical properties, Kahnafer et al. [10]. These correlations have been implemented in Fluent using a user defined function.
1.2.1 Density
The density of nanofluid is based on the physical principle of the mixture rule [10].
= (1 )Pf + V„P„,
(8)
where f and p refer to the fluid and nanoparticle respectively, q>p is the solid volume fraction of the nanoparticles, and p is the density.
1.2.2 Viscosity
Different models of viscosity have been used by researchers to model the effective viscosity of nanofluid as a function of solid volume fraction. For low solid volume fraction, Einstein's model can be used to predict the viscosity of the nanofluid.
It is worth mentioning that Einstein's model underestimates the nanofluid viscosity and does not consider the effect of temperature variations [34]. In this study, general correlation formulas to determine the effective viscosity of Al2O3-water considering the temperature effect proposed by Nguyen et al. [35] are used.
Meff = (1.125 - 0.0007 x T , Vp = 1%, 20 < T [ oC]< 70, (9)
= (2.1275 - 0.0215 x T + 0.0002 x T 2)^f,
yp = 4%, 20 < T[ oC]< 70. (10)
1.2.3 Thermal Conductivity
1.2 Thermo-Physical Properties of A|2O3-Water Nanofluid
The nanofluid mixture is considered as a single phase; thus, the nanoparticles and the base fluid are in a thermal equilibrium with each other and the relative velocity is negligible or equal to zero. Therefore, the effective thermophysical properties that depend mainly on the temperature, volume concentration of the nanoparticles, and the properties of the base fluid and the suspended particles are adopted in this study. General correlations for the effective thermal conductivity, the viscosity, and thermal expansion
Numerous studies were conducted in the literature to model the thermal conductivity of nanofluid. A general thermal conductivity correlation for Al2O3-water nanofluid was proposed by Kahnafer et al. [10] for various temperatures, nanoparticles diameter, and solid volume fraction. The proposed correlation is given as follows:
^ = 1.0 +1.0112® +
kf ®
+Z4375® (
d (nm)
0.0248^ (■
p 0.613
). (11)
1.2.4 Thermal Expansion Coefficient
The effect of temperature and solid volume fraction on the thermal expansion coefficient of Al2O3-water nanofluid were investigated by Ho et al. [36]. They developed a correlation for the thermal expansion coefficient of Al2O3 water nanofluid as a function of the temperature and the solid volume fraction of nanoparticles. This correlation is given as:
4 7211
peff = (-0.-79^ + 9.3158 x10-3T - -T^) x10-3
0 <çp < 0.04, 10 < T[°C]< -0. (12)
In this study, the values of the thermal expansion coefficient are evaluated utilizing Eq. (12) and based on the average temperature between cold and hot walls.
1.2.5 Heat Capacity
The specific heat of nanofluid is formulated by assuming a thermal equilibrium between the nanoparticles and the base fluid phase as follows [10]:
_(1 -Vp )Pfcf +Vp Ppcp
peff
(13)
where pp is the density of the nanoparticle, pf is the density of the base fluid, pf is the density of the nanofluid, and cp and Cf are the heat capacity of the nanoparticle and the base fluid, respectively.
2 NUMERICAL METHODS
The finite volume method is employed using Fluent 6.3 commercial software to solve the governing equations subject to specified boundary conditions. Gambit commercial software is used for the generation of grids. A high number of cells are constructed near the surface of cavity walls to compensate for the high velocity gradient in the boundary layer region of the viscous flow. The coupling between the pressure and velocity fields is achieved using Simple scheme. A second order upwind scheme is used for the mixed convection. The coupled equations are solved sequentially. Time independent solver was used for all the simulation. Laminar model was used to simulate the mixed convection flow. A second order discretization scheme was used for all simulations. All the thermophysical properties material described above are implemented in the Fluent software using user defined function (UDF). The simulation is terminated when the residuals for continuity and
momentum equations attain 10-6 and the residual for the energy equation attain 10-8. The local Nusselt number based on the length of the square cavity is expressed as:
Nu = -k-f Ê°\ kf dY
(14)
The average Nusselt number of the hot wall can be obtained by integrating the local Nusselt number along the wall as:
_ i
Nu = J Nu(X )dX. (15)
3 RESULTS AND DISCUSSION
In order to ensure that computational results are grid-independent, different grid sizes were tested (see Table 2).
Table 2. Grid independence results for pure water (p = 0)
Number of grids 50 x 50 80 x 80 100x100
Average Nusselt number at the hot wall for Ri = 0.1 7.0124 7.294 7.345
Average Nusselt number at the hot wall for Ri = 1.0 2.2897 2.345 2.346
Average Nusselt number at the hot wall for Ri = 10 1.653 1.671 1.677
The grid independence tests are performed for pure water (p=0). Table 2 illustrates the results of independence studies for average Nusselt number Nu at the hot wall for different Richardson numbers, Ri. From Table 2, it is noticeable that grid sizes of 80 x 80 and 100 * 100 give almost the same results for the average Nusselt number. For further validation, the present numerical solutions are validated against experimental as well as numerical results.
The numerical results are compared in terms of Nusselt number and ^-velocity at the mid-section of the cavity and found to be in good agreement with previous works (see Figs. 2 and 3). It is worth mentioning that when Ri ^ », the heat transfer by natural convection is dominant, while when Ri ^ 0 , the heat transfer by forced convection is dominant. In contrast, when Ri ^ 1, the heat transfer by the mixed convection is dominant. The Richardson number in this study varies from 0.1 to 10. Numerical simulations are conducted to demonstrate the effects of Richardson number Ri and solid volume fraction ^ on mixed convection heat transfer in a square cavity.
0
The value of the solid volume fraction varies from 0.0 to 0.04. The local Nusselt number distributions along the hot moving wall for various Richardson numbers and solid volume fractions are illustrated in Figs. 4 and 5. The local Nusselt number is high near the left wall and decreases towards the right wall. The effect of the solid volume fraction and the Richardson number is clearly discernible in these figures.
Fig. 2. Comparison of the local Nusselt number variation for hot wall for pure water and Ri = 1
Fig. 3. Comparison of U-velocity at mid-section of the cavity (X= 0.5) for pure water and Ri = 0.1
In general, the local Nusselt number started to decrease sharply at the beginning (i.e. at the left wall) and then continues to decrease rapidly at the right wall. It is clear from these figures that for a laminar flow with a constant solid volume fraction, decreasing the Richardson number leads to a significant increase in the local Nusselt number. The reason for this enhancement is that the forced convection is becoming dominant when the Richardson number is low. It can be seen in Figs. 5a and b that the local Nusselt number increases as the solid volume fraction y increases for
both Ri =10 and Ri =1. However, for Ri = 0.1 (see Fig. 5c) the Nusselt number does not change significantly with the solid volume fraction.
The variations of the average Nusselt number with Richardson number along the hot wall for
Fig. 4. Local Nusselt number for different solid volume fraction; a) y = 0.0, b) y = 0.01 and c) y = 0.04
Fig. 5. Local Nusselt number for different Richardson number; a) Ri = 10, b) Ri = 1, and c) Ri = 0.1
different solid volume fractions are depicted in Fig. 6. The figure shows that the average Nusselt number increases with the solid volume fraction and decrease with Richardson number.
Fig. 7 demonstrates the influence of the Richardson number on isotherms for the solid
Nu 5
-ip = 0.0
..............ip =0.01
---ip = 0.004
.......... '
o.i l id
Ri
Fig. 6. Average Nusselt number verses Richardson number for different solid volume fraction
volume fraction y = 0.04. For Ri = 0.1, shear forces are dominant in the cavity and, thus, the isotherms are mostly congregated to the area near the bottom surface of the square cavity. The high temperature gradient near the bottom wall is due to these high shear forces developed due to movement of the cavity lid. Moreover, a thin boundary layer is developed on the top of the hot wall. Furthermore, a strong circulation region is observed in the centre of rotating vortex inside the cavity. Increasing Ri up to 1, the shear forces and bouncy forces are the same magnitude, thus, coarse isotherms existed near the bottom wall. As a result, a moderate temperature gradient in the vertical direction can be seen. Further increases in Ri lead to wider expansion of isotherms covering the whole cavity. This leads to a further decrease of the heat transfer intensity near the bottom wall due to a stabilized free convection effect. It should be noted that for this case, the inertia force is weak. Therefore, the bouncy force defines the formed fluid flow and heat transfer inside the cavity.
Fig. 8 demonstrates the effect of a solid volume fraction on isotherms under the natural convection condition (Ri = 10). For the case of pure water (y = 0), most of the entire cavity is in a thermally stratified state. In other words, the isotherms are almost parallel in the horizontal direction, except for the isotherm area in the right upper corner, where the isotherms show significant changes. Increasing y up to 0.01, the isotherm area in the right upper corner becomes bigger and the isotherms show significant changes. When the solid volume fraction is increased up to y = 0.04, the isotherms are clustered near the cold
Fig. 7. Isotherm for different Richardson number, p = 0.04; Fig. 8. Isotherm for different nanoparticle solid volume fraction,
a) Ri = 0.1, b) Ri = 1, and c) Ri = 10 Ri = 10; a) p = 0, b) p = 0.01, and c) p = 0.04
bottom wall. These changes are due to a high temperature gradient near the bottom wall. These changes have led to the increased thermal conductivity of nanofluid.
Fig. 9 demonstrates the influence of solid volume fraction on streamlines for the free convection case (Ri = 10). For Ri = 10, the bouncy force is dominant in the cavity and streamlines are mostly concentrated
a)
b)
c)
Fig. 9. Stream function for different nanoparticle solid volume fraction for Ri = 10 (kg/s); a) q = 0, b) q = 0.01, and c) q = 0.04
towards the top moving wall. Moreover, a denser distribution of streamlines near the top moving wall is observed for pure water (q = 0). Further increase of ^ leads to expansion of streamlines covering the whole
cavity due to the free convection enhancement. When the solid volume fraction is increased to y = 0.04, flow circulation is slightly increased. This can be explained as follows: increasing the solid volume fraction leads to increasing the bouncy force, which ultimately enhances the thermal conductivity.
5 CONCLUSIONS
The utilization of nanofluid in industrial applications has recently witnessed rapid growth due to their superior thermophysical properties. There is a need for a deeper and an accurate understanding of the effect of using nanoparticles in the combined convection heat transfer. Thus, this study focuses on the numerical investigation of steady, laminar mixed convection flow, and heat transfer in a lid-driven square cavity filled with water-Al2O3 nanofluid. The CFD model that includes the temperature variation effects on the thermophysical properties is validated against experimental results and found to be in good agreement. It was found that the local Nusselt number at Ri = 1 is augmented considerably by increasing the solid volume fraction of nanoparticles. At a relatively lower Richardson number (Ri = 0.1), the local Nusselt number does not change significantly with the presence of nanoparticles. It was concluded that when the solid volume fraction and Richardson number are varied, different typical patterns of isotherms are obtained. Furthermore, the streamline pattern in the square lid-driven cavity was almost completely described by primary recirculation, and there were no significant changes when solid volume fraction changes.
6 REFERENCES
[1] Cha, C.K, Jaluria, Y. (1984). Recirculation mixed convection flow for energy extraction. International Journal of Heat and Mass Transfer, vol. 27, no. 10, p. 1801-1810, DOI:10.1016/0017-9310(84)90162-5.
[2] Ghia, U., Ghia, K.N., Shin, C.T. (1982). High-Re solutions for incompressible flows using the Navier-Stokes equations and a multigrid method. Journal of Computational Physics, vol. 48, no. 3, p. 387-411, DOI:10.1016/0021-9991(82)90058-4.
[3] Senthilraja, S., Karthikeyan, M., Gangadevi, R., (2010). Nanofluid applications in future automobiles: comprehensive review of existing data. Nano-Micro Letters, vol. 2, no. 4, p. 306-310, DOI:10.1007/BF03353859.
[4] Imberger, J., Hamblin, P.F. (1982). Dynamics of lakes, reservoirs, and cooling ponds. Annual Review of Fluid Mechanics, vol. 14, p. 153-187, DOI:10.1146/annurev. fl.14.010182.001101.
[5] Kumar, A., Vemula, P.K., Ajayan P.M., John, G. (2008) Silver-nanoparticle-embedded antimicrobial paints based on vegetable oil. Nature Materials, vol. 7, p. 236-241, DOI:10.1038/nmat2099.
[6] Choi, S.U.S. (1995). Development and applications of Non-Newtonian flows. Singer, D.A., Wang, H.P. (eds.) Development and Application of non-Newtonian Flows, p. 99-105, ASME, New York.
[7] Vajjha, R.S., Das, D.K. (2012). A review and analysis on influence of temperature and concentration of nanofluid on thermophysical properties, heat transfer and pumping power. International Journal of Heat and Mass Transfer, vol. 55, no. 15-16, p. 4063-4078, DOI:10.1016/j. ijheatmasstransfer.2012.03.048.
[8] Khaled, A.-R.A., Vafai, K. (2005). Heat transfer enhancement through control of thermal dispersion effects. International Journal of Heat and Mass Transfer, vol. 48, no. 11, p. 21722185, DOI:10.1016/j.ijheatmasstransfer.2004.12.035.
[9] Torrance, K., Davis, R., Eike, K., Gill, P., Gutman, D., Hsui, A., Lyons, S., Zien, H. (1972). Cavity flows driven by buoyancy and shear. Journal of Fluid Mechanics, vol. 51, no. 2, p. 221-231, DOI:10.1017/S0022112072001181.
[10] Khanafer, K., Vafai, K., Lightstone, M. (2003). Buoyancy-driven heat transfer enhancement in a two-dimensional enclosure utilizing nanofluids. International Journal of Heat and Mass Transfer, vol. 46, no. 19, p. 3639-3653, DOI:10.1016/S0017-9310(03)00156-X.
[11] Jung, J.-Y., Cho, C., Lee, W.H., Kang, Y.T. (2011). Thermal conductivity measurement and characterization of binary nanofluids. International Journal of Heat and Mass Transfer, vol. 54, no. 9-10, p. 1728-1733, DOI:10.1016/j. ijheatmasstransfer.2011.01.021.
[12] Kang, H.U., Kim, S.H., Oh, J.M. (2006). Estimation of thermal conductivity of nanofluid using experimental effective particle volume. Experimental Heat Transfer, vol. 19, no. 3, p. 181191, DOI:10.1080/08916150600619281.
[13] Lee, S., Choi, S.U.-S., Eastman, J .A. (1999). Measuring thermal conductivity of fluids containing oxide nanoparticles. Journal of Heat Transfer, vol. 121, no. 2, p. 280-289, DOI:10.1115/1.2825978.
[14] Wang, X.-Q., Mujumdar, A.S. (2007). Heat transfer characteristics of nanofluids: a review. International Journal of Thermal Sciences, vol. 46, no. 1, p. 1-19, DOI:10.1016/j. ijthermalsci.2006.06.010.
[15] Pak, B.C., Cho, Y.I. (1998). Hydrodynamic and heat transfer study of dispersed fluids with submicron metallic oxide particles. Experimental Heat Transfer, vol. 11, no. 2, p. 151170, DOI:10.1080/08916159808946559.
[16] Xuan, Y., Li, Q. (2000). Heat transfer enhancement of nanofluids. International Journal of Heat and Fluid Flow, vol. 21, no. 1, p. 58-64, DOI:10.1016/S0142-727X(99)00067-3.
[17] Xuan, Y., Li, Q. (2003). Investigation on convective heat transfer and flow features of nanofluids. Journal of Heat Transfer, vol. 125, no. 1, p. 151-155, DOI:10.1115/1.1532008.
[18] Abu-Nada, E., Oztop, H.F. (2009). Effect of inclination angle on natural convection in enclosures filled with Cu-water nanofluid. International Journal of Heat and Fluid Flow, vol. 30, no. 4, p. 669-678, DOI:10.1016/j.ijheatfluidflow.2009.02.001.
[19] Tiwari, R.K., Das, M.K. (2007). Heat transfer augmentation in a two-sided lid-driven differentially heated square cavity utilizing nanofluids. International Journal of Heat and Mass Transfer, vol. 50, no. 9-10, p. 2002-2018, DOI:10.1016/j. ijheatmasstransfer.2006.09.034.
[20] Jang, S.P., Choi, S.U.S. (2004). Free convection in a rectangular cavity (Benard convection) with nanofluids. Proceedings of International Mechanical Engineering Congress and Exposition, Anaheim, p. 147-153, DOI:10.1115/ IMECE2004-61054.
[21] Jou, R.-Y., Tzeng, S.-C. (2006). Numerical research of natural convective heat transfer enhancement filled with nanofluid in rectangular enclosures. International Communication of Heat and Mass Transfer, vol. 33, no. 6, p. 727-736, DOI:10.1016/j. icheatmasstransfer.2006.02.016.
[22] Jafari, A., Mousavi, S.M., Tynjala, T., Sarkomaa, P. (2009). CFD simulation of gravitational sedimentation and clustering effects on heat transfer of a nano-ferrofluid. PIERS Proceedings, p. 1352-1356.
[23] Ismail, A.F., Rashmi, W., Khalid, M. (2008). Numerical study on buoyancy driven heat transfer utilizing nanofluid in a rectangular enclosure. Proceedings of the UK-Malaysia Engineering Conference, p. 118-123.
[24] Merilainen, A., Seppala, A., Saari, K., Seitsonen, J., Ruokolainen, J., Puisto, S., Rostedt, N., Ala-Nissila, T. (2013). Influence of particle size and shape on turbulent heat transfer characteristics and pressure losses in water-based nanofluids. International Journal of Heat and Mass Transfer, vol. 61, p. 439-448, DOI:10.1016/j.ijheatmasstransfer.2013.02.032.
[25] Iwatsu, R., Hyun, J.M., Kuwahara, K. (1993). Mixed convection in a driven cavity with a stable vertical temperature gradien. International Journal of Heat and Mass Transfer, vol. 36, no. 6, p. 1601-1608, DOI:10.1016/S0017-9310(05)80069-9.
[26] Mansour, M.A., Ahmed, S.E. (2013). Mixed convection in double lid-driven enclosures filled with Al2O3-water nanofluid. Journal of Thermophysics and Heat Transfer, vol. 27, no. 4, p. 707-718, DOI:10.2514/1.T4102.
[27] Abu-Nada, E., Chamkha, A.J. (2010). Mixed convection flow in a lid-driven inclined square enclosure filled with a nanofluid. European Journal of Mechanics - B/Fluids, vol. 29, no. 6, p. 472-482, DOI:10.1016/j.euromechflu.2010.06.008.
[28] Sharif, M.A.R. (2007). Laminar mixed convection in shallow inclined driven cavities with hot moving lid on top and cooled from bottom. Applied Thermal Engineering, vol. 27, no. 5-6, p. 1036-1042, DOI:10.1016/j.applthermaleng.2006.07.035.
[29] Oztop, H.F., Estellé, P., Yan, W-M., Al-Salem, K., Orfi, J., Mahian, O. (2015). A brief review of natural convection in enclosures under localized heating with and without nanofluids. International Communication in Heat and Mass Transfer, vol. 60, p. 37-44, DOI:10.1016/j.icheatmasstransfer.2014.11.001.
[30] Polidori, G., Fohanno, S., Nguyen, C.T. (2007). A note on heat transfer modeling of Newtonian nanofluid in laminar free convection. International Journal of Thermal Sciences, vol. 46, no. 8, p. 739-744, DOI:10.1016/j.ijthermalsci.2006.11.009.
[31] Putra, N., Roetzel, W., Das, S.K. (2003). Natural convection of nano-fluids. Heat and Mass Transfer, vol. 39, no. 8, p.775-784, DOI:10.1007/s00231-002-0382-z.
[32] Haddad, Z., Öztop, H.F., Abu-Nada, E., Mataoui, A. (2012). A review on natural convective heat transfer of nanofluids. Renewable and Sustainable Energy Reviews, vol. 16, no. 7, p. 5363-5378, DOI:10.1016/j.rser.2012.04.003.
[33] Haddad, Z., Abu-Nada, E., Öztop, H.F., Mataoui, A. (2012). Natural convection in nanofluids: Are the thermophoresis and Brownian motion effects significant in nanofluid heat transfer enhancement?, International Journal of Thermal Sciences, vol. 57, p. 152-162, DOI:10.1016/j.ijthermalsci.2012.01.016.
[34] Einstein, A. (1956). Investigation on the Theory of Brownian Movement. New York, Dover.
[35] Nguyen, C.T., Desgranges, F., Roy, G., Galanis, N., Mare, T., Boucher, S., Mintsa, H.A. (2007). Temperature and particle-size dependent viscosity data for water-based nanofluids - Hysteresis phenomenon. International Journal of Heat and Fluid Flow, vol. 28, no. 6, p. 1492-1506, DOI:10.1016/j. ijheatfluidflow.2007.02.004.
[36] Ho, C.J., Liu, W.K., Chang, Y.S., Lin, C.C. (2010). Natural convection heat transfer of alumina-water nanofluid in vertical square enclosure: An experimental study. International Journal of Thermal Sciences, vol. 49, no. 8, p. 1345-1353, DOI:10.1016/j.ijthermalsci.2010.02.013.
Strojniški vestnik - Journal of Mechanical Engineering 63(2017)6, 394-404 © 2017 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2017.4332 Original Scientific Paper
A Practical Method to Detect a Transverse Cracked Rotor
Using Transient Response
Xiaofeng Wang - Jun Liu* - Weimin Ge
Tianjin University of Technology, Tianjin Key Laboratory of the Design and Intelligent Control of the Advanced Mechanical System, China
To detect a transverse crack caused by fatigue or creep, most of the research has thus far paid attention only to resonances of steady-state oscillations created by the crack and proposed diagnosis systems utilizing these vibration phenomena. However, from a practical view point, these diagnosis systems have the following flaws: (1) the probability that a resonance occurs due to a crack in the rated rotational speed range is a lower position; (2) It is very dangerous to observe vibration characteristics in resonance ranges. In order to solve these problems, this paper uses a practical detection method utilizing the characteristic changes in a transient oscillation during the start-up, the shutdown, or the variable running speeds of rotating machinery. This method has great advantages, because it can check the occurrence signals of a crack in a wide speed range using a single sweep and avoid the operation in dangerous resonance ranges. Non-stationary characteristics during passages through the main resonance and various kinds of resonances are studied numerically and experimentally. Keywords: transient response, nonlinear rotor, cracked rotor, crack identification, experiments
Highlights
• A practical method for the diagnosis of the cracked rotor has been investigated theoretically and verified systematically in many experiments.
• The characteristics of the transient response have been shown on the cracked rotor, and there is not only the harmonic resonance, but also the sub-harmonic resonance and the super harmonic resonance are analysed under the transient state.
• A sensitive and accurate experimental setup has been developed with a high-quality data acquisition system to obtain precise measurement data.
• The experimental results have been verified well above the analyses from theoretical simulations.
0 INTRODUCTION
A transverse crack occurs due to fatigue, creep or both in the rotating machinery during operation, To prevent a serious accident caused by a cracked rotor, it is very important to discover the crack at the early stage of the crack propagation. Therefore, various kinds of diagnosis system have been developed to detect a crack in the rotating machinery. Verney and Green presented an on-line crack diagnosis regimen hinging on the accuracy of the crack model which should account for the crack's depth and location [1]. Silani et al. introduced a new finite element (FE) approach to detect small cracks and calculated the flexibility matrix of crack elements with modified integration limits [2]. Li and Chu developed an HHT signal processing technique on the AE feature extraction of natural fatigue cracks in rotating shafts [3]. Guang and Chen introduced a FE model for the crack identification of a static rotor with an open crack [4]. Ma et al. also used a FE model to calculate time-varying mesh stiffness for the effects of profile shift and tooth crack in a gear rotor system [5]. Darpe presented a novel method of the transient torsion excitation to detect fatigue transverse cracks in rotating shafts [6]. Xie et al. studied the motion
stability of the flexible rotor-bearing system under the unsteady oil-film force and other faults by calculating the maximum Lyapunov exponent of the system [7]. Zhu et al. [8] and Ishida and Inoue [9] theoretically and experimentally analysed the dynamic characteristics of a cracked rotor with an active magnetic bearing. Ferjaoui et al. investigated the effect of the presence of a transverse crack in a rotor supported by two hydrodynamic journal bearings [10]. In general, the detection methods of a crack are classified into two groups. One, such as [4], is the static examination where a rotating machine is dissolved, and the parts are examined independently. The other is the dynamic examination in which the changes of the vibration characteristics could be observed during the operation.
Aiming at a practical generator rotor, Chu and Wang reported that the magnitude of the harmonic resonance at the main critical speed increased and a super-harmonic resonance at the secondary critical speed occurred due to a crack [11]. Concerning other kinds of resonances, there is no report on practical rotors. However, Ishida et al. [12] and Ishida and Hirokawa [13] observed a sub-harmonic oscillation of order 1/2, a super sub-harmonic oscillation of order 3/2 and a summing-and-differential harmonic oscillation owning to a crack in the experimental
394
*Corr. Author's Address: Tianjin University of Technology, No. 391 Binshuixi Road, Tianjin 300384, China, liujunjp@tjut.edu.cn
setup. Guo et al. examined the identification of the early crack propagation with the empirical mode deposition (EMD) method [14]. Gomez et al. analysed the vibration signals based on energy utilizing the wavelet theory. The results demonstrated the good reliability of crack diagnosis with the 3* energy [15]. Ishida et al. investigated transient responses at the 1/2 order sub-harmonic oscillation [16]. Li and Zhang used the Hilbert-Huang Transform (HHT) to identify a crack in a rotor-bearing system under transient oscillations [17] and [18]. Wang etc. proposed the application of order tracking to investigate a crack when the rotor system has a varying speed [19]. Therefore, there is a possibility that those oscillations also occur in practical machines under stationary responses or transient responses.
Although the static method is more reliable, it requires much time and money. Therefore, the dynamic examination is preferable for the early detection of a crack. Most of the dynamic monitoring systems focus on the changes of vibration characteristics in the steady-state oscillation. They have the following defects from the perspective of practice: (1) Many symptoms due to a crack do not occur in the rotational speed range; therefore, it is impossible to detect them during the normal operation; (2) When some changes occur due to a crack, it is dangerous to investigate the characteristics of the resonance range because there is a possibility that the crack develops rapidly during the investigation.
This study focuses on a practical detection method using non-stationary vibrations to overcome those defects. When a rotating machine starts up or shuts down, the rotor sweeps all over the rotational speed range below the rated rotational speed. If these non-stationary data are used to detect a crack, defects (1) and (2) could be avoided. In addition, there are some studies to illustrate systematically the transient response with simulations.
In this paper, a typical open-close model is used to investigate the characteristics of non-stationary oscillations of a cracked rotor during the passages through the main critical speed and various kinds of subcritical speeds. In particular, the study focused attention on the influences of an angular acceleration and the magnitude and phase of an unbalance on the maximum amplitude. The effectiveness of this method is verified systematically through simulations and many experiments.
To solve the above problems, the next section proposes the theoretical modelling and the motion equations. The resonances of steady-state oscillations with a crack are investigated in Section 2. The
method to detect a crack using the change of the characteristics of non-stationary oscillations passing through the harmonic resonance is explained in Sections 3 and 4. Non-stationary oscillation during passages through a forward super-harmonic and a forward sub-harmonic resonance are interpreted in Section 5. The characteristics of non-stationary oscillations during passages through a variety of resonances are summarized and shown in Section 6. The experimental setup is explained in Section 7, and the experiment results are presented in Section 8. Finally, the concluding discussion is given in Section 9.
1 THEORETICAL MODELLING AND MOTION EQUATIONS
1.1 Theoretical Modelling and Spring Characteristics
The theoretical model and the coordinate systems are shown in Fig. 1. In the experimental setup mentioned in Section 7, the rotor system where the deflection and the inclination couple each other is a four-degree-of-freedom (4DOF). The disk is not located at the shaft centre. On the contrary, if the disk is located at the shaft centre, it can be divided into two separate 2DOF systems, that is, a deflection model (the Jeffcott rotor) and an inclination model. In the latter system, the natural frequencies change due to the gyroscopic moment similar to the 4DOF model [20]. Therefore, we use the inclination model in the theoretical analysis. The origin of the static Cartesian coordinate system O-xyz is at the midpoint of the bearing centreline (the connecting line of the right and left bearings). The z-axis coincides with the breath ring centreline. The inclination angle of the elastic shaft at the disk mounting position can be expressed by 6 and its projection angles of 6 to the xz- and yz-planes can be expressed by 6x and 6y , respectively. It is supposed that a crack appeared on the half of the shaft's length. The rotating coordinate system O-x'y'z' is also considered where the x'-axis coincides with the crack boundary. The projection angles of 6 to x'z- and y'z-planes can be represented by 9'x and Q'y, respectively. When the crack is opened and 0'y > 0, the stiffness of the shaft becomes small. When the crack is closed and 0'y < 0, the shaft stiffness returns to the same value as the shaft with no crack. Therefore, the restoring moment has the spring characteristics with a piecewise linearity and its components M'x and M'y in the x'z'-and y'z'-planes, respectively, are shown in Fig. 2. They are represented as follows.
-K=m
-M'y = (ô'2-A5'2 )6'y -M'y = (82 + A82)6'y
(6'y > 0) (6'y < 0)
(1)
where, and S' are the spring constants and AS2 is the directional difference.
Crack
^OPen)^
X X
Fig. 1. Model of a cracked rotor and coordinate systems
Fig. 2. Spring characteristics of a cracked rotor
1.2 Motion Equations
The ratio of the polar inertia moment to the diametric inertia moment can be represented by ip , the rotational speed and the damping coefficient can be represented by m and c, respectively. The dynamic unbalance's magnitude and its phase angle can be expressed by t and a, respectively. The rotational angle of the x'-axis is represented by y. Corresponding to the gravitational force, the constant moment M0 that works in the ^-direction is considered. The motion equations governing non-stationary oscillations in a symmetrical 2DOF inclination model with no crack is given by Ishida et al. [12] and Ishida and Yamamoto [20]. First, we transfer Eq. (1) into the expression of the stationary coordinate system. By replacing the part representing the restoring force in symmetrical system by this expression, we can obtain the non-dimensional motion equations for a cracked rotor as follows.
Ï + ipwOy + ipvOy + céx + (1T A 2)0X + (Aj ± A2)(éx coslrnt + éy sinlrnt) = = (1 - ip )t {2 cos (y/ +a) + y/ sin(i/ +a)}
éy - iP¥éx - ipy>ix + céy + (1T Ai)éy +
+(A1 ± A 2)(éx sin 2œt - éy coslrnt ) = = (1 - ip )t {y/2 sm(y/- + a ) - y/ cos(y/ + a) + M0
where
8' = (8!+ 52)/2,
A2 =A8'2/28'. As for the symbol
A1 = (Ô'-Ô2)/2Ô\ in the motion
equations, all upper signs are shown for 9' > 0 and all lower signs are also shown for 9' < 0.
The above motion equations have the dynamic characteristics as follows: (a) time-varying coefficients similar to the asymmetrical rotor, (b) rotating piecewise nonlinearity and (c) unbalance excitation.
2 STEADY-STATE OSCILLATIONS
The resonances of steady-state oscillations are investigated before studying non-stationary oscillations. When the rotational speed y = œ is a constant, the angular position of the x'-axis can be expressed by:
y = œt +
(3)
where y0 is the initial angle.
With this condition, Eq. (2) is integrated numerically by the Adams method. Let pf and pb are the natural frequencies of a forward and a backward whirling motions, respectively. In view of a vertical rotor, only a harmonic resonance [pf= m] appears in the vicinity of the main critical speed. In the following, the notation [20] is used to show the relationship between the natural frequency and the rotational speed when the resonance occurs. When an unbalance and the crack are on the same side, there exists an unstable range at the main critical speed. Otherwise, the unstable range will disappear [20].
Fig. 3 shows the case of a horizontal rotor. The symbol 0 represents the amplitude obtained numerically. These values are connected smoothly by a full line. Various kinds of resonances occur in a wide rotational speed range due to a crack in view of a horizontal rotor. In addition to the harmonic resonance [pf= m], the backward harmonic resonance [pb=-m], the super-harmonic resonances [pf= 2m] and [pf= 3m], the sub-harmonic resonance [pf= (1/2)m], the super-sub-harmonic resonance [pf= (3/2)m] and
the combination resonance [m = pf -pb] occur. This is because the equilibrium position of the rotor shifts due to gravity and, as a result, the rotor system has both characteristics of more complex nonlinearity and more complex parametric excitation. Within all the above resonances, only harmonic resonances of these [pf = m ] can vary these characteristics significantly, depending on the angular position of the unbalance.
3 NON-STATIONARY OSCILLATIONS (A VERTICAL ROTOR)
This section explains the non-stationary characteristics of a cracked rotor. The governing motion equations of it are given by putting M0 = 0 in Eq. (2). The acceleration of a rotor is a constant X, regardless of acceleration or deceleration. The angular position of the x'-axis can be obtained by:
y = (1/2)Xt1 +wt + y0. (4)
The response curves are shown by full lines in Fig. 4, and the curves are obtained numerically via the Adams method. The results for three kinds of angular acceleration X are also shown in Fig. 4. For
comparison, the amplitudes of steady-state oscillations (X = 0) are shown by the symbol
When the unbalance and the crack are on the same side, the result is shown in Fig. 4a. In this case, an unstable range exists. If the rotational speed of a rotor sweeps this unstable range, the large amplitude appears for the small angular acceleration X. If they are on the opposite side, another result is shown in Fig. 4b. Since there no exists an unstable range, the amplitude is comparatively small regardless of any value of the angular acceleration X. Fig. 5 shows that the maximum amplitude rmax changes with the X. If the unbalance and the crack are on the same side, the cracked rotor cannot pass the main critical speed range due to having the very large amplitude when the angular acceleration X is less than a certain critical value.
Different from the case of steady-state oscillations, the whirling speed and the rotational speed are different from each other, and the crack opens and closes repeatedly during the passage through the main critical speed. Therefore, it is imagined that the difference in the angular position of the unbalance does not influence the maximum amplitude in the case
1 2
I
<
/,=0.1 c = 0.02 [+2 ®] l i [ P P. CO]
Aj = 0,05 J [+(3/2)«] J > \ 1 (i :>'!
■ A2 = 0.05
M0 = -1.0 j
r = 0.1 i
a = 270° ft I 1 a[ t J
[+3®]$ o—i
0.5
1.0
1.5 2.0 2.5 Rotational speed \j/
Fig. 3. Amplitude variation curve of a horizontal cracked rotor
3.0
a)
b)
Fig. 4. Amplitude variation curves In the main critical speed with a vertical rotor; a) unbalance within the crack side, and b) unbalance without the crack side
of non-stationary oscillations. But Fig. 5 shows that the maximum amplitude rmax changes remarkably due to the direction of the unbalance. This means that the repetition of the opening and closing states do not occur frequently in the resonant range.
2 3 4
Angular acceleration /1 (/10")
Fig. 5. Maximum amplitudes with acceleration X
These results can be interpreted from the viewpoint of vibration diagnosis as follows: Since the maximum amplitude increases remarkably as shown in Fig. 4a, the appearance of a crack can be detected from the incremental amplitude. In contrast, if they are on the opposite side, it is very difficult to find it.
4 NON-STATIONARY OSCILLATIONS (A HORIZONTAL ROTOR)
In this section, non-stationary oscillations during the passage through the main critical speed are investigated when the rotor system is supported horizontally.
Time histories obtained by numerical integration are more complicated than that of a vertical rotor. A time history is shown in Fig. 6a for the case that the unbalance and the crack are on the same side. In addition, spectrums obtained by the complex-FFT method [20] are shown in Fig. 6b, in which the positive abscissa represents the forward whirling motion and the negative abscissa represents the backward whirling motion. Because of the rotor passing through the unstable range, the amplitude changes remarkably. Different from a harmonic component case of a vertical rotor, many frequency components exist in the spectrum diagram. In addition to the harmonic component between 1 and 1.5, a constant component, a backward harmonic component and a forward super-harmonic component with a frequency of two times the rotational speed coexist. Therefore, it is impossible to obtain the amplitude of the harmonic component by
calculating ,J(6x )2 + (6y )2. Instead, the complex-FFT method is used to process these data.
Based on the steady-state oscillations at the main critical speed of a cracked horizontal rotor [9] in our previous study, we obtained the following results: (a) in the case of the comparatively large unbalance, there exists an unstable range if the unbalance and the crack are on the same side and it disappears if they are on the opposite side; (b) in the case of the comparatively small unbalance, an unstable range appears regardless of the angular position of the unbalance. Therefore, we discuss non-stationary oscillations in the following two cases.
4.1 Case with a Large Unbalance
Fig. 7 shows amplitude variation curves with a comparatively large unbalance. The value of the unbalance is the same as that of Fig. 4. If the unbalance and the crack are on the same side, the maximum amplitudes are very large as shown in Fig. 7a. In contrast, if they are on the opposite side, the maximum amplitudes are small as shown in Fig. 7b.
a)
£)J Frequency
Fig. 6. Non-stationary oscillation of a horizontal rotor; a) time history with non-stationary and b) spectrum with non-stationary
Fig. 7 is almost the same for a vertical rotor shown in Fig. 4. There is no qualitative difference between a vertical rotor and a horizontal rotor as the large unbalance. Therefore, the appearance of a crack
can be similarly detected only if the unbalance and the crack are on the same side.
a)
b) Rotational:
Fig. 7. Responses [+«] of a horizontal rotor with a large; a) unbalance within the crack side, and b) unbalance without the crack side
2.5
2.0
1.5
1.0
0.5
0.0
ip =0.25
c = 0.025 o
r = 0.005 2 = 0.0001
A, =0.02
A, = 0.02
a=270° | \ A =0.0005
M„=-1.0 /
\ y\ 2 = 0.0010
1
P \W_Ay \ -
0.9
1.0 1.1
12
1.3 1.4 1.5
Rotational speed i// Fig. 8. Responses [+«] of a horizontal rotor with a small unbalance
4.2 Case with a Small Unbalance
Fig. 8 shows amplitude variation curves with a relatively small unbalance when the unbalance and the crack are on the opposite side. As seen from
the motion equations, the cracked rotor has both characteristics of a forced vibration system owing to the unbalance and a parametrical excited system owing to the breathing mechanism of the crack and the static deflection. If the rotor is well balanced, the effect of the former diminishes then the latter appears predominantly due to the breathing of the crack and the static deflection. As a result, in the steady-state oscillation, an unstable range always exists, no matter what the angular position of the unbalance. Therefore, the maximum amplitude of the non-stationary oscillation is large and it does not vary according to the angular position of the unbalance. Fig. 8 corresponds to Fig. 7b. These figures indicate that the amplitude becomes large when the unbalance is not large. The results consider that when the rotor system is well balanced, the occurrence of a crack can be monitored by the incremental amplitude to pass through the main critical speed with any directional unbalance.
5 FORWARD SUPER-HARMONIC AND FORWARD SUB-HARMONIC RESONANCES
5.1 Forward Super-harmonic Resonance of Order 2
In the neighbourhood of the rotational speed y = 0.52 in Fig. 3, twice the rotational speed is almost equal to a forward natural frequency and the super-harmonic resonance [pf=2oi\ occurs. This resonance does not occur in a linear symmetrical rotor with a circular cross section if there is not a crack.
When the rotor passes through this resonance, the amplitude change is shown in Fig. 9. Because this super-harmonic resonance occurs more rarely than the harmonic resonance, comparatively large values of parameters Aj and A2 are used in the calculation. The maximum amplitude is almost independent to the angular position of the unbalance and does not decrease rapidly when the angular acceleration X increases. These characteristics show that this kind of resonance can be used as a correct signal for the occurrence of a crack.
5.2 Forward Sub-Harmonic Resonance of Order 1/2
In the neighbourhood of the rotational speed y = 2.23 in Fig. 3, the sub-harmonic resonance of order 1/2 [pf = 1/2w] occurs. This resonance does not occur in a linear symmetrical rotor if there is no crack. The amplitude variation curves are shown in Fig. 10a and the relationships between rmax and X are shown in Fig. 10b. Different from the cases of the harmonic resonance and the super-harmonic resonance, the
Rotational speed /
Fig. 9. Amplitude curve of super-harmonic resonance [+2a]
resonance curves of this steady-state oscillation
bifurcate from a trivial solution with the zero amplitude. In such a case, it is known that the
amplitude variation curve can change remarkably
according to the initial angular position of the
unbalance [16]. This means that the maximum
amplitude takes various values with a given
acceleration X. As a result, the relationship between rmax and X is represented by the shaded zone in Fig. 10b. Since this initial angular position cannot be controlled generally, the maximum amplitude of a non-stationary oscillation changes randomly for every operation. In Fig. 10b, the maximum amplitude decreases rapidly with the angular acceleration X increasing. Based on the above characteristics, it is considered that this resonance is not suitable to be used as a signal of the occurrence of a crack in the non-stationary oscillations. Because this resonance does not always occur, the method is still not suitable for detecting a crack. In order to find the occurrence of a crack by observing this resonance, it is necessary to conduct several running tests with the small angular acceleration X.
6 COMPARISON OF CHARACTERISTICS OF NON-STATIONARY OSCILLATIONS THROUGH VARIOUS RESONANCES
As mentioned above, the characteristics of non-stationary oscillations during the passage through resonances is different depending on the kind of
<
2.5 2.0 1.5 1.0 0.5 0.0
ip = 0.25
c= 0.01
j-= 0.1
A! =0.1
A2=0.1
A= 0.0002 °
ce= 270°
M0=-1.0
a)
2.6 2.7 2.8 2.9 3.0 3.1 3.2 Rotational speed y/
Fig. 10. Amplitude curve of sub-harmonic resonance [+1/2rn]; a) resonance curve, and b) maximum amplitude rmax for X
a)
b)
Fig. 11. The relationships between rmax and X; a) main critical speed with a small unbalance, and b) different cases of various kinds of critical speed
resonance. Here, the maximum amplitude rmax depending on the angular acceleration X is summarized and compared.
Fig. 11a shows the case of the main critical speed of the horizontal rotor with a small unbalance. Since an unstable range exists at any directional unbalance, the maximum amplitude always becomes very large if the rotor passes through the main critical speed with a small angular acceleration. This means that, as mentioned above, the harmonic resonance can be used to detect a crack when the rotor is well balanced. When the unbalance is comparatively large (this case is not shown by the figure), the maximum amplitude during the passage through the main critical speed range does not increase due to the crack as the unbalance is located in the opposite side of a crack. Therefore, in this case, we must investigate the change of the maximum amplitude by changing the unbalance's direction.
Fig. 11b shows the relationships between rmax and X to pass through various kinds of resonances. In this figure, each plot representing the maximum amplitude is obtained by drawing a response curve like Fig. 4a. And these plots are connected smoothly by full lines. In order to investigate the quantitative difference among these resonances, the same parameter values are used. In the case of various harmonic resonances where the maximum amplitude depends on the initial angular position, a result for the same initial angular position is shown in this figure. From Fig. 11, the following results are confirmed: (a) if the angular acceleration X is small, every kinds of resonance can be used as the signal of the occurrence of a crack; (b) if the angular acceleration X is comparatively large, we can use only the super-harmonic resonance [+2«] as a signal of the occurrence of a crack.
Fig. 12. Experimental system of a vertical rotor
7 EXPERIMENTAL SETUP
For safety, it was impossible to perform experiments on non-stationary oscillation using a horizontal rotor. Therefore, a vertical rotor system is used in experiments.
The experimental setup is shown in Fig. 12. A disk is mounted on the shaft at the positions of 200 mm or 150 mm below the middle of the elastic shaft that is simply supported using two double-row ball bearings at both shaft ends. The shaft is 700 mm in length and 12 mm in diameter. The diameter and the thickness of the disk are 481 mm and 5.5 mm, respectively. In order to make the same characteristics as in a transverse crack, a notch is made at the position 335 mm from the upper bearing, and it is filled with a part of the same dimensions as the notch as shown in Fig. 12. The notch has the width rn = 20 mm and the depth h = 5 mm. The phase angle a in Fig. 12 represents the angle between the crack and the unbalance. The deflections in the x- and j-directions are measured by position sensors.
8 EXPERIMENTAL RESULTS
8.1 Amplitude Variation during the Passage through the Resonance
The experimental setup where a disk is mounted on the positions of 200 mm below the middle of the elastic shaft is used.
8.1.1 Unbalance within the Crack Side
The unbalance is set a = 110° within the crack side. The experiments are made for eight different angular accelerations X = (75, 60, 50, 42.9, 37.5, 33, 30 and 27.3) rpm/s. The experimental results are represented by solid lines in Fig. 13. For comparison, the oscillation curve of the steady-state response is also shown by the symbol The unstable range drawn by the hatching appears between rn = 1015 rpm and rn = 1040 rpm when the rotor runs with a constant speed. When the angular acceleration X is 27.3 rpm/s, the maximum amplitude is large. When the angular acceleration X equals 75 rpm/s, the rotor passed the unstable range with small amplitude. The experimental results shown in Fig. 13 agree qualitatively with the numerical results shown in Fig. 4a.
8.1.2 Unbalance without the Crack Side
The unbalance is set a = 259° without the crack side. The experiments are made for six different angular accelerations k = (75, 60, 50, 42.9, 37.5 and 33) rpm/s. The experimental results are shown by solid lines in Fig. 14. There is no unstable range under the steady-state oscillation. Therefore, when the rotor passes through the main critical speed, the maximum amplitude is always small. The experimental results shown in Fig. 14 agree qualitatively with the numerical results shown in Fig. 4b.
900 1000 1100 1200 1300
Rotational speed