^ Since 1955 _ Strojniški vestnik Journal of Mechanical Engineering no. 9 year 2015 volume 61 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: Grafex, d.o.o., printed in 380 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 Vice-President of Publishing Council Jože Balič University of Maribor, Faculty of Mechanical Engineering, Slovenia Strojniški vestnik Journal of Mechanical ^jjfr Engineering Cover: The tool-wear of cutting tools has a very strong impact on the product quality as well as efficiency of the machining processes. Despite the nowadays high automation level in machining industry, tool-wear is usually measured off the machine tool. This paper presents an innovative, robust and reliable direct measuring procedure for measuring spatial cutting tool-wear in-line, with the usage of laser profile sensor. Courtesy: Laboratory for Cutting, Faculty of Mechanical Engineering, University of Ljubljana, Slovenia ISSN 0039-2480 International Editorial Board Kamil Arslan, Karabuk University, Turkey 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 Janez Kušar, UL, Faculty of Mechanical Engineering, Slovenia Gorazd Lojen, UM, Faculty of Mechanical Engineering, Slovenia Thomas Lübben, 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 Popović, 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 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. © 2015 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 61, (2015), number 9 Ljubljana, September 2015 ISSN 0039-2480 Published monthly Papers Luka Čerče, Franci Pušavec, Janez Kopač: A New Approach to Spatial Tool Wear Analysis and Monitoring 489 Matej Müller, Gorazd Novak, Franc Steinman, Gašper Rak, Tom Bajcar: Influence of the Operating and Geometric Characteristics of a Bottom-hinged Flap Gate on the Discharge Coefficient of a Side Weir 498 Mohsen Moslemi, Mohammadreza Khoshravan: Cohesive Zone Parameters Selection for Mode-I Prediction of Interfacial Delamination 507 Andrzej Milecki, Dominik Rybarczyk: Modelling of an Electrohydraulic Proportional Valve with a Synchronous Motor 517 Yu Zhang, Hongzhi Yan, Tao Zeng: Computerised Design and Simulation of Meshing and Contact of Formate Hypoid Gears Generated with a Duplex Helical Method 523 Jernej Laloš, Tomaž Požar, Janez Možina: High-Frequency Calibration of Piezoelectric Displacement Sensors Using Elastic Waves Induced by Light Pressure 533 Primož Potočnik, Ervin Strmčnik, Edvard Govekar: Linear and Neural Network-based Models for Short-Term Heat Load Forecasting 543 Strojniški vestnik - Journal of Mechanical Engineering 61(2015)9, 489-497 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2512 Original Scientific Paper A New Approach to Spatial Tool Wear Analysis and Monitoring Luka Čerče* - Franci Pušavec - Janez Kopač University of Ljubljana, Faculty of Mechanical Engineering, Slovenia The tool wear of cutting tools has a very strong impact on the product quality as well as on the efficiency of the machining processes. Despite the current high automation level in the machining industry a few key issues prevent complete automation of the entire turning process. One of these issues is tool wear, which is usually measured off the machine tool. Therefore, its in-line characterization is crucial. This paper presents an innovative, robust and reliable direct measurment procedure for measuring spatial cutting tool wear in-line, using a laser profile sensor. This technique allows for the determination of 3D wear profiles, which is an advantage over the currently used 2D subjective techniques (microscopes, etc.). The use of the proposed measurement system removes the need for manual inspection and minimizes the time used for wear measurement. In this paper, the system is experimentally tested on a case study, with further in-depth analyses of spatial cutting tool wear performed. In addition to tool wear measurements, tool wear modelling and tool life characterization are also performed. Based on this, a new tool life criterion is proposed, which includes the spatial characteristics of the measured tool wear. The results of this work show that novel tool wear and tool life diagnostics yield an objective and robust methodology allowing tool wear progression to be tracked, without interruptions in the machining process or in the performance of the machining process. This work shows that such an automation of tool wear diagnostics, on a machine tool, can positively influence the productivity and quality of the machining process. Keywords: machining process, tool wear measurement, spatial tool wear, in-line monitoring, volumetric estimator, tool life prediction Highlights • Newly developed measuring system to determine spatial cutting tool wear in-line. • Spatial cutting tool wear measurements were performed on a case study. • Tool wear modelling and tool life characterization has been performed. • A new tool life criterion is proposed. 0 INTRODUCTION AND STATE OF THE ART The tool wear of cutting tools has a very strong impact on product quality as well as on the efficiency of machining processes overall. Despite the current high automation level in the machining industry, a few key issues prevent complete automation of the entire turning process. One of these issues is tool wear, which is usually measured off the machine tool and is still done by hand under a toolmaker's microscope. There is no any industry wide procedure for automating the process of measuring wear. Furthermore, such conventional wear measurement requires stopping the automated turning, removing the tool, measuring the tool and putting the tool back to the holder, which is a considerable time loss relative to the tool's life. Therefore, the in-line characterization of cutting tool wear is crucial for cutting cycle times and costs, as well increasing the overall efficiency of the machining process. This paper presents an innovative, robust and reliable direct measuring procedure for measuring spatial cutting tool wear in-line using a laser profile sensor. This technique allows 3D wear profiles to be determined, which is an advantage over the currently used 2D subjective techniques (microscopes, etc.). The use of the proposed measurement system removes the need for manual inspection and minimizes the time used for wear measurement. A tool life characterization is also performed. Based on this, a new tool life criterion is proposed that includes the spatial characteristics of the measured tool wear. In the next three subsections, extensive state reviews are performed. 0.1 Tool Wear Characterization Review The damages to a cutting tool are the consequence of the stress state and thermal loads on the tool surfaces, which in turn depend on the cutting mode, i.e. turning, milling or drilling, cutting parameters and the cooling/ lubrication conditions. During machining, the cutting tool wear mechanisms and their rates are very sensitive to changes in the cutting operation and the cutting conditions. To minimize machining costs, one must not only to find the most suitable cutting tool and work material combination for a given machining operation, but also reliably predict the tool life and also the tool changing/replacement protocol. Tool wear mainly occurs at the rake and flank face. Flank wear is caused by friction between the flank face of the cutting tool and the machined workpiece surface and leads to loss of the cutting edge. Therefore, flank wear affects the dimensional accuracy and surface finish quality of the product. In practice, flank wear is generally used as the cutting tool wear criteria. When critical value of tool wear has been reached, the cutting tool is assumed to have reached its tool life. This point should be placed before but close to the moment when the cutting tool fails due to excessive stresses and thermal alterations. To avoid this, the cutting tool must be replaced before reaching its critical limit. In practice, the preferred cutting tool life criteria is the tool-flank wear threshold. However, even though the wear should progress gradually and should thus be easily monitored, tool wear is hard to predict due to uncertainties in the machining process. Therefore, to robustly define the tool-changing protocol in NC (numeric control) programs, the inline tool wear measurement is inevitable [1] and [2]. In practice, some directly measured dimensional characteristics and criteria of typical wear patterns, i.e. crater, flank wear, and depth-of-cut notch wear at the extremities, for HSS (high speed steels), carbide and ceramics tools, are standardized in ISO 3685 [3], as shown in Fig. 1. Fig. 1. Typical wear pattern according to ISO 3685 The process of cutting tool wear consists of three characteristic parts: the initial (running-in) period, the longest uniform (progressive) wear period and the accelerated wear period leading often to catastrophic failure [1]. The machining process needs to be, on the one hand, stopped at the right time to prevent undesired consequences of the tool wear, such as increases in cutting forces, vibrations, noise, and temperature in the cutting zone, as well as deviation of part dimensions and surface quality from the point of view of respective tolerance values, while on the other hand, the tool should be replaced as late as possible to increase efficiency and decrease tooling costs. In reality, several wear mechanisms occur simultaneously, although any one of them may dominate the process. They can be qualitatively identified as mechanical, thermal and adhesive. Mechanical types of wear, which include abrasion, chipping, early gross fracture and mechanical fatigue, are basically independent of temperature. Thermal loads appear with plastic deformation, thermal diffusion and oxygen corrosion as their typical forms, and increase drastically at high temperatures, thus accelerating the tool failure by causing easier tool material removal (by abrasion or attrition) [4]. At high cutting speed, temperature-activated wear mechanisms including diffusion (solution wear), chemical wear (oxidation and corrosion wear), and thermal wear (superficial plastic deformation due to thermal softening effect) occur. By contrast, at lower cutting speeds adhesive and abrasive wear are the most significant types of wear. ISO 3685 notes that the type of tool wear is difficult to determine when machining with grooved tools. Ee et al. [5] found that the tool wear pattern on a worn grooved tool insert is influenced by the three-dimensional chip-flow and by the complex chip-groove configurations. They also found that some grooved tools fail long before the major flank wear (VBB) reaches its failure criterion in many cases. It is also shown, that grooved tools demonstrate tool-failure as a result of concurrent multiple tool wear parameters. This confirms that the tool wear, or rather tool life, is significantly affected by the combined effects of cutting conditions and chip-groove configurations. This also answers the question as to why several new measurable parameters have been proposed for tool wear in machining with grooved tools (Fig. 2) [5]. From the grooved tools wear parameters, it can be seen that the wear regions in grooved tools can be very specific, such as the notch, etc. However, in general they can be divided into three major wear regions. These are: (a) the edge wear region, (b) the land and secondary face wear region, and (c) the backwall wear region (Fig. 2). Fig. 2. Measurable tool wear parameters in grooved tools 0.2 Tool Wear Measurement Review Tool wear can be measured using direct measuring techniques or estimated by indirect measuring techniques [6] and [7]. In indirect measuring techniques, tool wear is estimated using other more measurable machining process variables such as cutting force, surface finish, acoustic emission, vibrations, energy consumption, temperature, etc. [8]. A survey of the literature indicates that a variety of approaches have been applied to tool wear prediction [9] to [12]. Conversely, with direct measuring techniques, condition monitoring is carried out by analysing the change in geometry of the cutting tool. Usually these techniques can be taken only between two sequential machining runs, because the major flank of the tool is not visible during actual machining. Direct tool condition monitoring techniques can be divided into two-dimensional and three-dimensional techniques. Additionally, they can be divided into measurements made directly on the machine tool or outside it (using a microscope, etc.). Kurada and Bradley [13] have designed a system consisting of a fiber-optic light source to illuminate the tool and a CCD camera, which is used in combination with a high resolution video zoom microscope. The resolution in the captured image was approximately 3 дт per pixel. Images are captured orthogonally to the flank face. The actual wear region identification and measurement relies on the specular reflection of the wear area from the fiber-optic lights. Identification of the tool wear area is performed in 2D, through the wear area reflection. Lanzetta [14] has a different approach to the problem. Using two cameras, a grey-scale image is taken from the nose of the tool and image of the tool silhouette is taken from the top of the tool. The silhouette is taken with a back-light and the tool nose is directly illuminated. Images are also taken with a directional light from the top of the tool to detect crater wear. Flank wear is measured from the silhouette image. The silhouette of the worn tool is compared to the silhouette of the unused tool, which gives an estimate of the tool flank wear. Measurement error is reported to be less than 5 %. Niranjan Prasad and Ramamoorthy [15] measured the crater wear using stereo vision. Stereo vision is not achieved by using two cameras but with the tool moving between two images taken with the same grey-scale camera. Resolution is 12 дт x 10 дт per pixel. The depth and width of the crater wear is measured and the measurements can be made accurately. However, the images of the tool are taken outside the turning machine and no method for automatically measuring the flank wear is presented. On the other hand, Jurkovic et al. [16] used structured light from a laser diode to capture elevation information from the tool surface. A monochrome CCD camera was used to capture the images with a halogen light. The resolution of the image was approximately 3 дт per pixel. Flank wear can be measured as well as crater wear using the elevation information. The system was combined with a CNC lathe. Additionally, Devillez et al. [17] shows that white light interferometry is a useful tool for measuring cutting insert crater wear. The measurement methodology allows examination of all aspects of crater wear after the metal cutting process, although this technique does not permit in-line tool crater measurement and flank wear measurement. Dawson and Kurfess [18] present a novel method for quantifying tool wear using three-dimensional computational metrology that is also based on white light interferometry. The presented measurement technique allows measurement of flank and crater wear. Furthermore, a technique was developed to measure the deviations between the three-dimensional tool data and an 'ideal' version of an unworn cutting tool. The calculated deviations were then used to define the tool wear volume. Wang et al. [19] employ a phase shifting method using fringe patterns to measure crater wear by constructing a 3D map of the tool insert. Four fringe patterns with various phase-shifts are projected onto the rake face of the tool upon which four grey-level images are captured. With appropriate setup, the system can be used for on-line crater wear measurement in industry, but the accuracy level is limited. A combination of the fringe projection system and the white light interferometry is used by Weckenmann and Nalbantic [20]. The main disadvantages of the described methods are the inability to measure wear profiles in depth (three-dimensional geometry of crater wear - KT, etc.) and/or that they cannot be robustly used directly in the machine tool. To perform the measurement with these methods, the cutting insert should be removed from the machine tool. This causes time loss and possible problems with the accuracy of subsequent processing. This paper presents a novel method belonging to direct measuring method of determining of cutting tool wear. This method offers the possibility of three-dimensional tool wear measurement directly on the machine tool, without the need to remove cutting inserts from the tool holder (machine tool). 0.3 Tool Life Assessment Review Flank wear of cutting tools is often selected as the tool life criterion because it determines the diametric accuracy of machining, as well as its stability and reliability [21]. According to the ISO 3685:1993 [3], the assessment of flank wear is determined by its direct measurement. Flank wear is considered using an aged Taylor's tool life equation having a phenomenological nature. The Taylor Tool life estimation is based on experimental data that are gathered by measuring the time until the tool life criterion is fulfilled using several different cutting speeds. Based on this experimental set, a curve (regression model) representing the tool life can be constructed. The drawback is that it describes just the linear section of the tool life curve. Additionally, Taylor's basic equation does not include the effect of cutting feed and is limited to a certain range of speeds. It also does not give information about the tool wear at any particular time t during turning. Even though an improved Taylor model exists that includes both feed and depth of cut, but they still just cover the tool life, not the wear at particular times. An additional problem of such a model is that it is constructed based on tool wear measurements performed by the operator, etc. These assessments are subjective and insufficient. They do not include the tool geometry (the flank angle, the rake angle, the cutting edge angle, etc.) so they are not suitable for comparing cutting tools having different geometries. They do not include the cutting regime and thus do not reflect the real amount of work material removed by the tool during the time over which the measured flank wear is achieved [21]. The selection of these characteristics depends upon the particular objective of a tool wear study. Most often, dimensional accuracy dictates this selection, i.e. the need to manufacture parts within the tolerance limits assigned for tool wear. As such, the tool life defined by this criterion may be referred to as dimension tool life. Dimension tool life can be characterized by the time within which the tool works without adjustment or replacement (Tc-l); by the number of parts produced (N^); by the length of the tool path (Lc-l); by the area of the machined surface (Ac-l) and by the linear relative wear (hl-r) [21]. All these characteristics listed are particulars and thus, in general, do not allow for the optimal control of cutting operations, comparison of different cutting regimes, assessment of different tool materials, etc. For example, dimensional tool life is of little help if one needs to compare cutting tools that work at different cutting speeds and feeds and/or when the widths of their flank wear land are not the same. The problem lies in the complex cutting tool surface geometries as well as high stochastic behaviour of the tool wear from the point of view of aggression and location. Therefore, the combination of the dimensional wear rate with the relative surface wear and the specific dimensional tool life are much more general characteristics to be used in metal cutting tests conducted everywhere from the research laboratory to shop floor level. A new 3D tool wear based estimation is presented in Section 2. 1 NEWLY DEVELOPED MEASUREMENT SYSTEM The newly developed measuring system (Fig. 3) consists of a high-accuracy 2D laser displacement sensor Keyence LJ-G015 with a proper controller Keyence LJ-G5001 [22], a motorized linear translation stage Standa 8MT173-DCE2 and a LabVIEW application developed for process control. The measurement range of the 2D laser displacement sensor Keyence LJ-G015 is in Z-axis (height) ±2.6 mm and 7.0 mm in X-axis at the reference distance. The repeatability in the Z-axis is 0.2 ^m and 2.5 ^m in the X-axis. For linear positioning of the measurement head, a motorized linear translation stage is used with a minimal incremental motion of 0.1 ^m and bi-directional repeatability of 0.4 ^m [23]. With movement of the profile sensor across the cutting tool and the support of the developed software (LabVIEW application), the profile data are grabbed and prepared in a matrix form for further evaluation and analysis. A laser displacement sensor measures the distance from the measurement head to the points projected on the measured object. In this way the Z-coordinates of the point cloud are measured. The X-coordinate is defined by the specification of the laser displacement sensor [22], while the Y-coordinate represents linear stage feed direction. While the relative position/orientation of the measuring head has a significant influence on the measurements, the analysis of its orientation has been experimentally and empirically analysed in one of our previous works [24]. The accuracy was determined to reach ± 8 (tm. The measurement system was mounted on Mori Seiki SL-153 CNC lathe as shown in Fig. 3. This allows a quick measurement of tool wear on the machine tool itself. Between two machining operations, the average measuring procedure can be performed in approximately 20 s, without stopping the machining process. Fig. 3. Measurement system mounted on a Mori Seiki machine tool Thus, we can take a measurement that includes far more information about tool wear, without having to remove the cutting inserts from the holders, which would be necessary if we were to measure tool wear on a toolmaker's microscope. 2 TOOL LIFE PREDICTION FOR GROOVED CUTTING INSERTS Tool life prediction in grooved tools involves simultaneously occurring multiple wear parameters, as presented in Fig. 2 (wear on flank, wear on rake, etc.). However, there are no universally accepted standards for tool life prediction in grooved tools. The ISO 3685 based criterion, established for flat-faced tools, is usually not implementable on grooved tools. The geometry of the cutting tool significantly affects tool life, as the geometry defines the magnitude and direction of the cutting force and its components, the sliding velocity at the tool-chip interface, the distribution of the thermal energy released in machining, the temperature distribution in the cutting edge, etc. [21]. Therefore, a new tool life prediction method is presented based on the rate of cutting tool geometry change and consequent change in the stresses at the cutting edge. New cutting edge -Used cutting edge Fig. 4. Measured cutting tool wear geometry The greatest impact on prediction of the tool life for finishing operations is the radial wear hr (Fig. 4). Thus, we can say that the cutting tool is useful until the radial wear hr exceeds the prescribed tolerance of the workpiece. With the previously presented measurement system, current radial wear hr can be measured directly on the machine tool. In this way we can control, measure and predict the cutting tool life for finishing operations. In roughing operations, productivity is important (i.e. material removal rate), so the tool breakage is under consideration. For this purpose the influence of mechanical load on the rake face of the cutting insert will be simulated using the finite element method (FEM). Tamizharasan and Kumar [25] use FEM for simulating the effect of tool geometry on flank wear, surface roughness and cutting forces. In our case, the 3D model of the cutting insert is defined on the basis of measurements results (KT, VB, KE and yef) performed in a case study. The mechanical pressure distribution on the rake face is applied based on the Zorev model [26]. The pressure decreases nonlinearly from the maximum value at the cutting edge to zero at the end of chip-tool contact zone (where the chip separates from the rake face). Apart from the model of contact load distribution on the tool faces, the contact load also depends on the contact area of the tool with the chip and the workpiece. On the rake face, the contact area depends on the contact length between the chip and tool rake face, lcr, as shown in Fig. 5, and chip width. The length can be measured from the rubbing marks on the rake face, which can be seen from the 3D wear measurements. The case study will be shown in the next section. As a benefit of these FEM simulations, stress analyses show how changes in the geometry (wear) bring the tool to the stress limit/structure breakage. Based on this, the tool life can be predicted and with it the undesired tool breakage in real machining operations. -New cutting edge-Used cutting edge Fig. 5. a) mechanical pressure distribution on the cutting tool edge, and b) 3D model of cutting insert with pressure distribution on the rake face As a proof of the proposed concept, the measurement system presented in Section 1 was tested in a machining environment. The results of a case study were analysed as described above and are presented in the next section. 3 CASE STUDY AND RESULTS The presented measurement system has been tested on a case study of bearing steel machining (100Cr6 -AISI52100) performance determination. A bar shaped workpiece with diameter 40 mm and length 290 mm was machined. Machining tests were conducted on a Mori Seiki SL-153 turning machine. Commercially available Sumitomo DCGT 11T304 R-FX cutting inserts, grade ACZ310, were used with SDCJR 2020K11 tool-holder. All the machining experiments were performed dry, while the cutting parameters have been defined according to the manufacturer's recommendations: ap = 1 mm, vc=250 m/min and f= 0.12 mm/rev. To follow the progress of the wear, the sequential longitudinal workpiece turning operations were performed on 20 mm length intervals. The result of the corresponding tool wear progress is presented in Fig. 6. The measurement itself is executed in approximately 20 seconds, grabbing the 6 mm * 3.5 mm area, with a grid size of 0.01 mm * 0.005 mm (AX * Д Y). To determine tool wear progress, six repeated measurements were performed. The first picture shows the new (unworn) cutting insert, while the sixth picture show the critical wear that occur on the rake face of the cutting tool after 7.4 minutes of machining. The colour on those plots corresponds for the wear magnitude. Red color shows the presence of BUE (build up edge), while the blue color shows the actual wear. Nevertheless, the BUE regions were not dominant. Furthermore, from the result it can be seen that abrasive, adhesive and diffusion wear occurred on the cutting edge. Abrasive and diffusion wear can be observed on the rake face (region of blue colour), while adhesive wear is seen in the growing region of yellow colour. The progression of crater wear on rake face and chipping of cutting edge can also be observed (Fig. 6, darker blue zone on cutting edge). Observing the rake face crater wear depth formation, with a growing region of darker blue colour, high temperatures in cutting area are suspected. t=5.7 min ^^^P t~6.8 min Ig';- ' t~7.4 i Fig. 6. 3D deviation results (new vs. worn cutting tools) Fig. 7. Measurement of typical wear patterns, according to ISO 3685, on a cross section profile, 0.5 mm from secondary flank face after 7.4 min of cutting For analysis of the extent and shape of the wear, an example of a cross section perpendicular to the primary flank face is shown and analysed in Fig. 7. The wear corresponds to the 7.4 min of machining. The depth of flank wear is clearly visible from the comparison of cutting inserts cross-sections profiles (new vs. worn). Fig. 7 shows that with the use of the presented measurement system, the extent of the wear formation can be precisely determined. The depth of the crater for this particular case is KT = 0.115 mm and the flank wear extends to a depth VB = 0.058 mm. The 3D deviation results showed that the maximum wear occurs 0.5 mm from the secondary flank face. Therefore, these profiles have been analysed and the results of flank wear land width ( VB) and total wear volume (VWear) as a function of cutting time are presented in Fig. 8. Initial (running- Uniform (progressive) Accelerated in) wear period wear period wear period r. / / / / Ф 0.03 E D О > 0.02 го CD 0.01 ГО .О H 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 Cutting timet [min] -*-Tool wear VB [mm] -♦-Total wear volume V Fig. 8. Flank wear land width (VB) for cross section profile 0.5 mm from the secondary flank face and total wear volume (VWear) as function of cutting time t E 0.10 E, fe 0.08 Initial (running-in) wear period Uniform (progressive) wear period Accelerated wear period / / / / / 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 Cutting time t [min] Fig. 9. Maximum crater depth (KT) as a function of cutting time t for the cross section profile 0.5 mm from secondary flank face The results confirm that observing just the VB is not sufficient for tool life prediction. Higher wear than VB is observed on the rake face (although values are significantly smaller than KT). The VB will result in higher cutting forces and missing tolerances, while KT will cause cutting tool breakage. In practice, the most important consideration in selecting cutting tools and cutting conditions is tool life (T). T is defined as the time when the measured wear levels exceed a predefined threshold. A standard measure of tool life is the time to develop a flank wear length VBB of 0.3 mm or VBBmax of 0.6 mm, in case of cemented carbides. Standard ISO 3685:1993 also recommends a maximum value of KT (KTmax). It can be empirically calculated via following equation [2]: KTmax = 0.06 + 0.3 • f , (1) where f is feed. For the cutting parameters used, the KTmax is 0.096 mm. From the results presented in Fig. 9 it is apparent that KT is critical, while the critical value for VB is not reached (Fig. 8) (VB = 0.06 mm after 7.5 min of machining). However, this empirical KTmax equation does not include a lot of important parameters in machining (material, depth of cut, etc.) and is questionable. Therefore, this paper also aims to present a new approach in tool life prediction using FEM simulations. For this purpose, geometry data from cross-section profiles of the worn cutting inserts have been used to build an accurate 3D model of the cutting edge. This model has been used for stress analysis. The maximum pressure p = 0.95 GPa has been calculated from cutting forces. This has been performed with a back calculation of online measured cutting forces. The cutting tool material properties presented in Table 1 were found in the literature [27] and [28] and used in this work. Table 1. Cutting tool material properties Elastic modulus 600000 N/mm2 Poisson's ratio 0.225 Shear modulus 263000 N/mm2 Mass density 14.5 kg/m3 Tensile strength 450 N/mm2 Compressive strength 6000 N/mm2 Maximum crater depth (KT) as a function of cutting time is shown in Fig. 9. In all cases, three characteristic parts are observed: the initial (running-in) period (up until 1.5 min), the longest uniform (progressive) wear period (from 1.5 min to 5.5 min) and the accelerated wear period (from 5.5 min on). The results of the FEM simulation are presented in the Fig. 10. In Fig. 10a, the results of stress analyses of the new cutting insert are presented, while Fig. 10b presents results for the critically worn insert. The wear corresponds to the 7.4 min of machining. The maximal stress for a new insert (a) is 7.3*108 N/m2 (yellow zone), while for the worn inserts (b) it is much bigger 1.1x109 N/m2 (red zone). combination (or any other wear mechanism), where the cutting tool breakage can be foreseen. У I.IOOe I.OOBe 9.167e 8.250e 7.333e 6.417e 5.500e 4.583e 3.667e 2.750e 1.833e 9.167e З.ОООе Fig. 10. FEM simulation results for a) new (unworn) and b) worn cutting insert The presented results showed that FEM simulations are in good correlation with the experimental data presented in Fig. 8 (total volume VWear) and Fig. 9 (KT). The stress at the cutting tool edge increases drastically when crater wear becomes significant (from 5 min to the end), as can be seen in Fig. 11. This happens due to the notch effect. By increasing the depth of crater wear, the cross-section of the cutting edge became smaller, and therefore the stress in the material increases. This leads to catastrophic tool fracture (breakage) when the stresses became higher than the material yield strength. The cutting time limit before the breakage occurs defines the cutting tool life. 1.15E+09 1.10E+09 1.05E+09 1.00E+09 „ E 9.50E+08 ž 9.00E+08 » 8.50E+08 8.00E+08 7.50E+08 7.00E+08 E .1,0.08 0.00 12 ! * t—'-J Wr 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 Cutting timet [min] -»-Crater depth KT [mm] -»-Stress Fig. 11. FEM simulated stress and crater depth KT as function of cutting time t for cutting tool edge On the other hand, maximal crater depth KT can be predicted from the yield strength and actual stress. With maximal KT thus defined the presented spatial measurement system can be used for the robust control of the process and its efficiency. Going a step further, this procedure can also cover the KT and VB 4 CONCLUSIONS This paper presents an innovative, robust and reliable direct measuring procedure for determining spatial cutting tool wear in-line on machine tools. The technique allows for the determination of 3D wear profiles, an advantage over the currently used 2D subjective techniques (microscopes, etc.). The use of the proposed measurement system removes the need for manual inspection and minimizes the time used for wear measurement. The measurement data thus obtained are actually realistic 3D models of the cutting tool/wear, where we can clearly see what the existing conditions of the cutting tool are. With information on three-dimensional deviations, this method outperforms traditional 2D deviation methods in accuracy, efficiency and reliability. Besides tool wear, this work also upgrades existing tool life methodologies with spatial tool wear information. A new methodology for tool life prediction, based on FEM stress simulation is presented. Because it is known that the geometry of the cutting tool directly affects tool life, cross-section profiles of the worn cutting inserts have been used to build an accurate 3D model of the cutting edge. For this, the measurement system has been experimentally tested on a case study machining 100Cr6 bearing steel on a CNC lathe. The cutting tool wear has been measured directly on the machine tool, without the need to remove the cutting insert from the machine tool and analyse it under the microscope. In addition to the presented conventional results (VB, KT), this new spatial tool wear measuring system, allows tool wear volume VWear to be measured and analysed. A constructed FEM model was run to predict tool life. This model compared unworn and worn insert geometries and used stress analysis for the FEM. The mechanical loading was same for worn and unworn inserts and was experimentally defined. From the obtained stress results, the cutting time limit before breakage occurs (tool life) can be more accurately defined/predicted. These predictions for the case study are also correlated with the fractures occurring on the cutting tool. The results show that FEM simulations can be used to determine the time limit before fracturing of the cutting tool/tool life. The results of this work show that novel tool wear and tool life diagnostics represent an objective and robust method of tracking tool wear progress, without interrupting the machining process, as well as the performances of the machining process itself. Additionally, the tool life/tool fracturing can be robustly predicted in-line. To summarize, our work shows that the automation of tool wear diagnostics inline on a machine tool can beneficially influence the productivity, efficiency and quality of the machining process. 5 REFERENCES [1] Cheng, K. (2009). Machining Dynamics, Fundamentals, Applications and Practices. Springer, London. [2] Kopač, J. (2008). Machining: Theoretical Basis and Technological Information. Faculty of Mechanical Engineering, Ljubljana. (in Slovene) [3] ISO 3685 (1993): Tool-Life Testing with Single-Point Turning Tools. International Organization for Standardization, Geneva. [4] Grzesik, W. (2008). Advanced Machining Processes of Metallic Materials. Elsevier B.V., Amsterdam. [5] Ee, K.C., Balaji, A. K., Jawahir, I.S. (2003). Progressive tool wear mechanisms and their effects on chip-curl/chip-form in machining with grooved tools: an extended application of the equivalent toolface (ET) model. Wear, vol. 255, no. 7-12, p. 1404-1413, D0l:10.1016/S0043-I648(03)00112-1. [6] Kurada, S., Bradley, C. (1997). A review of machine vision sensors for tool condition monitoring. Computers in Industry, vol. 34, no. 1, p. 55-72, D0I:10.1016/S0166-3615(96)00075-9. [7] Dutta, S., Pal, S. K., Mukhopadhyay, S., Sen, R. (2013). Application of digital image processing in tool condition monitoring: A review. CIRP Journal of Manufacturing Science and Technology, vol. 6, no. 3, p. 212-232, DOI:10.1016/j. cirpj.2013.02.005. [8] Siddhpura, A., Paurobally, R. (2013). A review of flank wear prediction methods for tool condition monitoring in a turning process. The International Journal of Advanced Manufacturing Technology, vol. 65, no. 1-4, p. 371-393, D0I:10.1007/s00170-012-4177-1. [9] Dimla Snr, D.E. (2000). Sensor signals for tool wear monitoring in metal cutting operations-a review of methods. International Journal of Machine Tools and Manufacture, vol. 40, no. 8, p. 1073-1098, D0I:10.1016/s0890-6955(99)00122-4. [10] Govekar, E., Gradišek, J., Grabec, I. (2000). Analysis of acoustic emission signals and monitoring of machining processes. Ultrasonics, vol. 38, no. 1-8, p. 598-603, D0I:10.1016/S0041-624X(99)00126-2. [11] Teti, R., Jemielniak, K., O'Donnell, G., Dornfeld, D. (2010). Advanced monitoring of machining operations. CIRP Annals - Manufacturing Technology, vol. 59, no. 2, p. 717-739, D0I:10.1016/j.cirp.2010.05.010. [12] Venkata Rao, K., Murthy, B.S.N., Mohan Rao, N. (2013). Cutting tool condition monitoring by analyzing surface roughness, work piece vibration and volume of metal removed for AISI 1040 steel in boring. Measurement, vol. 46, no. 10, p. 4075-4084, DOI:10.1016/j.measurement.2013.07.021. [13] Kurada, S., Bradley, C. (1997). A machine vision system for tool wear assessment. Tribology International, vol. 30, no 4, p. 295-304, DOI:10.1016/S0301-679X(96)00058-8. [14] Lanzetta, M. (2001). A new flexible high-resolution vision sensor for tool condition monitoring. Journal of Materials Processing Technology, vol. 119, no. 1-3, p. 73-82, DOI:10.1016/S0924-0136(01)00878-0. [15] Niranjan Prasad, K., Ramamoorthy, B. (2001). Tool wear evaluation by stereo vision and prediction by artificial neural network. Journal of Materials Processing Technology, vol. 112, no. 1, p. 43-52, DOI:10.1016/S0924-0136(00)00896-7. [16] Jurkovic, J., Korosec, M., Kopac, J. (2005). New approach in tool wear measuring technique using CCD vision system. International Journal of Machine Tools and Manufacture, vol. 45, no. 9, p. 1023-1030, DOI:10.1016/j. ijmachtools.2004.11.030. [17] Devillez, A., Lesko, S., Mozer W. (2004). Cutting tool crater wear measurement with white light interferometry. Wear, vol. 256, no. 1-2, p. 56-65, DOI:10.1016/S0043-1648(03)00384-3. [18] Dawson, T.G., Kurfess, T.R. (2005). Quantification of tool wear using white light interferometry and three-dimensional computational metrology. International Journal of Machine Tools and Manufacture, vol. 45, no. 4-5, p. 591-596, DOI:10.1016/j.ijmachtools.2004.08.022. [19] Wang, W.H., Wong, Y.S., Hong, G.S. (2006). 3D measurement of crater wear by phase shifting method. Wear, vol. 261, no. 2, p. 164-171, DOI:10.1016/j.wear.2005.09.036. [20] Weckenmann, A., Nalbantic, K. (2003). Precision measurement of cutting tools with two matched optical 3D-sensors. CIRP Annals - Manufacturing Technology, vol. 52, no. 1, p. 443-446, DOI:10.1016/S0007-8506(07)60621-0. [21] Astakhov, V.P. (2004). The assessment of cutting tool wear. International Journal of Machine Tools & Manufacture, vol. 44, no. 6, p. 637-647, DOI:10.1016/j.ijmachtools.2003.11.006. [22] Keyence (2010). High-accuracy 2D Laser Displacement Sensor. Keyence. [23] Standa. Translation Stage with DC motor, from http://www. standa.lt, accessed on 2014-21-2 [24] Cerce, L., Pusavec, F., Dugar, J., Kopac, J. (2013). Further development of the spatial cutting tool wear measurement system. Journal of Production Engineering, vol. 16, no. 2, p. 5-10. [25] Tamizharasan T., Senthil Kumar N. (2012). Optimization of Cutting Insert Geometry Using DEFORM-3D: Numerical Simulation and Experimental Validation. International Journal of Simulation Modelling, vol. 11, no. 2, p. 65-76, DOI:10.2507/ IJSIMM11(2)1.200. [26] Zorev, N.N. (1966). Metal Cutting Mechanics. Pergamon Press, Oxford. [27] Zhao, Z., Hong, S.Y. (1992). Cryogenic properties of some cutting tool materials. Journal of Materials Engineering and Performance, vol. 1, no. 5, p. 705-714, DOI:10.1007/ BF02649252. [28] Sandvik. Cemented Carbide, from http://www.sandvik.com, accessed on 2015-2-10. Strojniški vestnik - Journal of Mechanical Engineering 61(2015)9, 498-506 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2453 Original Scientific Paper Influence of the Operating and Geometric Characteristics of a Bottom-hinged Flap Gate on the Discharge Coefficient of a Side Weir Matej Müller1* - Gorazd Novak2 - Franc Steinman3 - Gašper Rak3 - Tom Bajcar4 1 DHD Ltd, Digital Hydrodynamics, Slovenia 2 Hidroinštitut, Institute for Hydraulic Research, Slovenia 3 University of Ljubljana, Faculty for Civil and Geodetic Engineering, Slovenia 4 University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Bottom-hinged flap gates on side weirs are often used for the regulation of flow diversion in case of water abstraction for a variety of needs. In this study a new equation for the discharge coefficient of a bottom-hinged flap gate on a side weir was proposed on the basis of discharge measurements. The equation was divided into two parts. The first part covers the impact of the sharp-crested side weir and the second the influence of the position and the width of the flap gate. In this manner, the discharge coefficient can be calculated with other authors' equations for a sharp-crested side weir, which then must be multiplied by the new proposed coefficient. Very good agreement was found between the obtained discharge coefficients and the one calculated with the proposed equation. Furthermore, the results were compared with the equations of other authors for the discharge coefficients of sharp- and broad crested side weirs. The agreement was found to be very good. Additionally, measurements of water levels along the edge of the flap gate and measurements of the velocity field were carried out with a computer-aided visualization method. From these measurements, it was possible to show that the contraction of the water jet varies with the gate-opening angle. It was also found that the side weir with a flap gate has the most favorable hydrodynamic shape around the gate-opening angle of 33°, where the discharge coefficient reaches its maximum. Keywords: side weir, discharge coefficient, flap gate, overflow gate, visualization method Highlights • A new coefficient is proposed for the correction of the discharge coefficient of a sharp-crested side weir due to the opening gate angle of a bottom-hinged flap gate on the side weir. • Contraction of the water jet and the discharge coefficient of a side weir with a flap gate are dependent on the gate-opening angle of the flap gate on a side weir. • The side weir with a flap gate has the most favorable hydrodynamic shape around the gate-opening angle of 33°, where the discharge coefficient reaches its maximum. • The discharge coefficient of a flap gate on a side weir has its minimum at the gate-opening angle of 0°, where it can be considered as a broad-crested side weir. 0 INTRODUCTION Side weirs are structures used for flow diversion in case of water abstraction for a variety of needs such as irrigation or for flood water diversion. The latter structures are also known as flood overflows and have proven to be highly useful for the diversion of flood water in wastewater treatment plants as well as for the diversion of water from river channels, for example, when excess water has to be diverted from a river channel to a retention area. For the proper design of such structures, it is important to know the flow capacity of side weirs, which has already been investigated by various authors. If we define dQ as the flow through an infinitely small stripe ds along the side weir, the general equation for the specific flow over the side weir can be written as follows [1]: q=-f=(h - p)A, (1) ds 3 where (h -p) represents the overflow height over the side weir, g represents the gravitational acceleration, and Cd is an empirical coefficient. Several equations have been proposed for the calculation of the coefficient Cd for the sharp-crested side weir and for the broad-crested weir. Most of the authors have considered the impact of the inflow Froude number in their equation [2] and [3], some have also considered the impact of the overflow height [4] and [5], but only a few have taken into account the impact of the length of the side weir. The impact of the width of a broad-crested side weir has also been investigated [6] and [7]. Fig. 1 shows the geometric and hydraulic parameters that have been proven to have the greatest impact on the coefficient Cd [2]. Fig. 1. Sketch of a) a sharp-crested rectangular side weir, and b) a broad-crested side weir In Fig. 1, h1 represents the water level of the approaching flow in the main channel, h2 the water level of the outflow in the main channel downstream of the side weir, W the width of the flap gate or of the broad-crested side weir, B the width of the main channel, L the length of the side weir, Q1 the discharge of the approaching flow in the main channel, Q2 the discharge of the outflow in the main channel downstream of the side weir, and Qs the discharge over the side weir. The need for the regulation of the lateral discharge over the side weir often appears, for example, with the regulation of the water level in the reservoir in the emergence of a flood wave. Therefore, bottom-hinged flap gates are often used on side weir structures to regulate the amount of the lateral discharge over the side weir at different reservoir water levels. Good knowledge about the discharge capacity of such structures is necessary for the appropriate design and operation of flap gates. Several hydraulic analyses of such gates on weirs have been made by various authors [8] and [9], although available studies of the discharge capacity of flap gates are virtually non-existent. Furthermore, the hydraulic conditions at such structures lack deep analyses, for example of the three-dimensional velocity field. In this paper, a new equation for the coefficient Cd of a flap gate on a side weir is proposed; in addition, measurements of the velocity field in the horizontal plane have been made to explain the obtained dependence of the Cd from the gate-opening angle ф. Fig. 2 shows the characteristic parameters of the flap gate on the side weir used in this study. W denotes the width of the flap gate and ф the gate-opening angle. Fig. 2. Sketch of a rectangular side weir with a flat bottom - hinged flap gate 1 METHODS Within the experimental work, discharge, water surface and velocity field measurements were performed. The discharge measurements were performed with the help of a triangular Thomson weir on the upstream end of the model (inflow of water into the model) and a rectangular weir that was placed inside of the channel, on the downstream end of the model (outflow of water out of the main channel). Since such a weir for the measurement of discharges, should be wider than 0.3 m, according to the ISO standard [10], a rating curve of the weir was obtained with the help of the upstream Thomson weir (which was in accordance with the ISO standard) and with the side weir closed. It was found that the correlation between the obtained rating curve and the one according to the ISO standard is very good, since the maximum difference in flow rates does not exceed the value of 1.4 % (for spill heights greater than 3.5 cm, for which the measurements were conducted). In addition to the discharge measurements, the upstream water levels were also measured to determine a new equation for Cd. The water levels were measured with a point gauge of 0.1 mm precision. Furthermore, measurements of the velocity field in a horizontal plane located just above the edge of the flap gate were carried out. The measurements were made with the computer-aided visualization method based on the advection-diffusion equation. This method is implemented in ADMflow software, developed jointly by University of Ljubljana, Faculty of Mechanical Engineering and Abelium [11]. The method has been verified with a number of cases; it has already been used for the measurements of the velocity fields of a free surface flow specific for a case of a side weir. It was found that the most suitable tracer for the measurements in such cases are hydrogen bubbles, which can be produced with galvanization [12]. It has been shown that the uncertainty of the method is ±5 % [11]. 2 EXPERIMENTAL SET-UP A physical model of the side weir with a bottom-hinged sharp edged flat flap gate was constructed and placed in a 6 m long rectangular channel with the width of 0.2 m. The flap gate was hinged on the sharp-crested side weir and supported with a threaded metal rod, which was then used to fine-tune the position of the flap gate (angle ф). A 0.25 mm thick copper wire was placed approximately 0.15 m upstream of the upstream edge of the side weir. The wire was attached to a specially designed fork, which was constructed out of a 2 mm thick electrically isolated wire. Both the thin wire and the fork were placed at an appropriate distance from the side weir; in this way the water flow in the area of the side weir was not disturbed. The thin copper wire could then be adjusted to the height at which the velocity field was measured at that particular time (the height position of the wire was adjusted to the height of the illuminated plane). At the same location (about 0.15 m upstream of the upstream edge of the side weir), a copper plate of negligible thickness (< 1 mm) was placed on the bottom of the channel. The copper wire and the plate were connected to a power supply; in this way, a galvanization process started and hydrogen bubbles were produced, which served as a tracer. Fig. 3. A sketch of the established physical model and the measurement apparatus The velocities were measured in the illuminated plane, 5 mm above the overflow edge of the gate. A height-adjustable weir was placed in the main channel about 1.5 m downstream of the downstream edge of the side weir, which served for the regulation of the downstream water in the main channel. The position of this weir influenced the upstream conditions and was required to obtain a variety of different Froude numbers and water levels, and thus affected the size of the lateral outflow over the side weir. In actuality, such a regulation would be represented with a gated weir in the river channel. The fact that the impact of the surface tension is negligible at such model dimensions and that the main channel was already shown in previous studies [12], in which the measurements of the hydraulic conditions at a side weir with a similar length as used in this study made in a 150 mm wide channel were additionally verified with the measurements made in a wider channel (0.5 m) with a longer side weir. The research showed a good correlation between the measurements in the narrow and wider channel. In this study, various lengths L of the side weir (0.15 m, 0.2 m and 0.25 m) and various heights of the weir p (0.05 m, 0.075 m and 0.1 m) were considered. Furthermore, all lateral discharge measurements were made for different upstream Froude numbers F1 (0.1 to 0.41) and different overflow heights h1 - p (0.03 m to 0.1 m), in which both were varied randomly. In all described variants, 7 different positions of the gate (different ф) were considered (0°, 10°, 20°, 33°, 45°, 70° and 90°). A total of 380 measurements with different variants of the above-described parameters were carried out. In every measurement, the inflow and outflow discharges as well as the upstream and downstream water levels were obtained. To confirm and explain the resulting dependence of Сф from ф, additional velocity measurements were carried out. In order to investigate this dependence, all hydraulic parameters were fixed for all the additional measured variants. At all 7 variants, F1 and h1 -p and the geometric properties (L, B, p) remained the same, only ф was varied (0°, 10°, 20°, 33°, 45°, 70°, and 90°). In addition to the velocity measurements in the horizontal plane, the upstream and downstream water levels and the water levels along the edge of the flap gate were measured for all 7 variants. 3 RESULTS The analysis of the discharge coefficient Cd of the flap gate on a side weir has shown that Cd is highly dependent of the opening gate angle ф. The approximate trend of the variation of Cd with ф can be recognized if all the measured Cd against ф are placed in a graph (Fig. 4). Fig. 4. From the measurements obtained Cd in dependence of ф The trend approximately follows a third-degree polynomial function. In this case, the correlation factor r2 is equal to 0.809. The fact that this trend follows a polynomial function was already shown [8] and [9]; however, their studies were made for flap gates on the weir that were placed perpendicular to the water flow and not on the side weir. Fig. 4 also shows the range of Cd values at each ф, which is between 10 % and 15 %. This range is the result of other influential hydraulic parameters of the water flow, such as Fj and hj, and other influential geometric parameters, such as L and p. All of these parameters were also considered in the new equation for Cd of the flap gate on a side weir. 3.1 The New Equation for Cd 3.1.1 Equation Derivation It was already shown that the Cd is mostly affected by the following ratios [2]: с, = f h p, LB, F ). (2) For the equation of Cd, a function was chosen in which the ratios from Eq. (2) were considered as the products of several power functions. Similarly, it was already shown that the Cd of a broad-crested side weir is, in addition to those already mentioned ratios from Eq. (2), also affected by the ratio [7]: = f (hi -p/W). (3) Because of the finding that the Cd in the case of a flap gate on a side weir is also affected by the value of the opening gate angle ф (Fig. 4), the latter was also included in the new equation as a polynomial function. A side weir with a flat, bottom-hinged flap gate can be considered to be a sharp-crested side weir without a gate, when the gate is fully closed (i.e. when ф = 90°). Therefore, the new equation for Cd was divided into two parts. In the first part (Cd0), the influence of the ratios of the parameters from the Eq. (2) were captured, and in the second part (C9) the influence of ф and the ratios from Eq. (3) were captured. The divided equation reads: Cd = Cd 0 • C4>, (4) where Cd0 denotes the discharge coefficient of the sharp-crested side weir and Cф denotes the influence of the flap gate. Before the new equation was developed, an analysis of the dependence of Cd from the parameters in Eq. (3) was conducted. For this purpose, trend lines were produced, which are shown in Fig. 5. It was found that the dependence of Cd on (hj -p/W) is relative to ф, as the trend line in Fig. 4 for ф = 0° indicates a downward trend, while the trend line for ф = 33° indicates an upward trend. The downward trend at broad-crested weirs has already been shown by other authors ([6] and [7]), while the upward trend at ф = 33° is quite logical, since the effect of the gate width W decreases with the increase of the overflow height (h1 -p) and is at some point negligible (at W/ (h1 - p)<<1). The trend line for ф = 90° is approximately horizontal, as in this case this ratio has no effect on the Cd. The different trends of Cd in dependence on the ratio W / (h j -p) at different values of ф were captured in the new equation for Cd with an additional exponent in Cr Fig. 5. Measured Cd in dependence of the ratio W/(ht-p) A similar equation was already produced [6] for the broad-crested weir, where the ratio W/ (h1 -p) was used as an exponent. However, in our case a logarithmic function was used for the additional exponent instead of a power function, as it was found to suit better the measured data, as well as such a function is quite logical, because it causes the coefficient Cd to stabilize from a certain value of the ratio W/ (hj -p) on. The proposed form of the new equation for the Cd of a bottom-hinged flap gate on a side weir is as follows: c-=-( I ) ( * ■(f ■Ф' + g Фк + h ф' + i ) i4n(W/ ( hj- p))+n (5) where ф is given in radians and the coefficients a, b, c, d, f, g, h, i, j, k, l, m, n represent unknowns, which were determined on the basis of all 380 measurements and with the use of the generalized reduced gradient (GRG) method. Table 1 shows the obtained values for all unknowns. With the parameters a to n, the correlation factor between the values that were obtained from the measurements and the values that were calculated from Eq. (5) is Cd, is r2 = 0.856. Table 1. Values of the unknowns In Eq. (5) a) Unknowns a to g A b c d f g 0.482 0.023 -0.080 0.150 0.3086 -1.7499 b) Unknowns h to n h i j k l m n 1.5359 0.9423 2.0996 1.4039 1.1314 2.7 2.0 , measured Fig. 6. Correlation between the measured and the calculated values of Cd Furthermore, the average difference between these values does not exceed 3.2 %, the maximal difference does not exceed 13.9 %, and the standard deviation is 0.06, which is similay to what other authors found in their studies [2]. Fig. 6 shows the correlation between the measured and the calculated values of Cd. According to Eq. (3), Eq. (4) can be divided into two parts: c' • = " ( p J Г ( f сф=( ■ФJ + g ■Фк + h■ф' + i ) i-ln(W/ ( h-p ))+n (6) (7) when the value ф is 1.5708 rad, which is equivalent to the angle 90°, Eq. (7) reduces to the value of 1, regardless to the values W, h j or p. Thus, Eq. (3) reduces itself to Cd=Cđh which is entirely logical, as the side weir with a flap gate in a fully closed gate position (i.e. at ф = 90°) is equivalent to a sharp-crested side weir and the ratio W/ (hj -p) does not have any impact on the value of Cd. equations developed from other authors showed very good agreement. Furthermore, the comparison of our measurements of Cd for the gate angle ф = 90° with those calculated after Eq. (6) also showed very good agreement, as the average difference does not exceed 2.8 %. The best agreement of our measurements with the equations from other authors was found by the comparison with the study from [4], in which the side walls in the side channel were also considered, as was done in our case. Furthermore, the measured values are in good agreement with the calculated values according to the equation in [13], where the average difference does not exceed 3.7 % and the maximum difference 8 %. The calculated values from the equations of other authors and our measured values of Cd are mainly located within the 10 % range [14] and [15]; however, it must be noted that other authors did not consider the side walls in the side channel. The agreement of our measured values of Cd and the calculated values by the equation in [5] and [16] is out of the 10 % range. However, it should be noted that the trends are very similar. Furthermore, the values calculated from the equations of other authors differ; for example, the mean difference between [5] and [15] is 35 %, and the difference between [3] and [5] is more than 50 % (Fig. 7). 1.0 0.9 0.8 0.7 0.6 0.5 О 0.4 0.3 o Eq. 5 Borghei (1999) X Jalili and Borghei (1996) Д Mohammed (2013) o Subramany and Awasthy (1972) + Yu-Tech, 1972 * Swamee et al (1988) Novak et al. 2012 , ...AM... ÀtÙ ___ ^ r? —Z f T x X * Г-^Ло"'0 : x 0.4 0.5 Cd, measured [-] 0.6 0.7 Fig. 7. Comparison of the measured Cd at ф = 90° and calculated from the Eq. (6) and from equations of other authors 3.1.3 Verification Using Other Experimental Data for Broad-Crested Weirs 3.1.2 Verification Using Other Experimental Data for Sharp-Crested Weirs The comparison of our measurements of Cd for the gate angle ф = 90° with the calculated values from As already mentioned, the side weir with the position of the flap gate at ф = 0° can be considered to be a broad-crested side weir, for which the width of the gate W represents the width of the broad-crested weir. Therefore, at values of ф = 0°, Eq. (5) reduces to an equation for a broad-crested side weir, and the coefficient Сф is in this case equal to: сф (ф = 0°) = CW = imMW/ ( р ))+и. (8) where CW denotes the correction factor for Cd for the case of a broad-crested side weir. The comparison of the calculated coefficients CW after Eq. (8) for the measured variants where the gate was fully opened (ф = 0°) with the calculated values after the equations for a broad-crested weir from other authors, showed good agreement. From the correlation of the values calculated with Eq. (8) and values calculated after the equation in [7], a correlation factor r2 of 0.97 was obtained and from the correlation with [6] a value of r2 = 0.95. Furthermore, in the first case, the average difference is less than 1.3 % and in the second, even lower than 1 %, while the maximum difference does not exceed 5 % in either one. It should be noted that the equation in [7] was made for a broad-crested side weir, and the equation in [6] was made for a broad-crested weir. Fig. 8 shows the correlation between the calculated values according to Eq. (8) and according to other authors for all measured variants where the gate was fully opened (ф = 0°). the gate W). Therefore, in this case, the impact of the flap gate on Cd is negligible. Fig. 8. Correlation between Сф according to the Eq. (8) and the calculated CW according to other authors for all measurements with the flap gate fully opened (ф = 0°) 3.1.4 Verification Using Other Experimental Data for Bottom-Hinged Flap Gates on Weirs Fig. 9 shows the dependence of C in dependence of W / (hj -p) turns at ф = 7.5°. According to Eq. (7), C(f is approximately 1 when the W / (hj -p) < 0.5 (i.e. in the case when the overflow height (hj -p) is two times greater than the width of Fig. 9. Dependence of Cv on W /(hj -p) according to the Eq. (7) and the comparison of Cv with the equations from other authors A much greater impact of the flap gate on Cd was found to be at higher values of the ratio W / (hj -p) (i.e. in the case when the overflow height (hj -p) is smaller than the width of the flap gate W), where Cv reaches a value of 1.3 at the angle ф = 33°. Fig. 9 also shows that the calculated values of C according to Eq. (7) and the equation in [9] was found to be less good, but still satisfactory. However, the authors did not consider the impact of the ratio W / (h1 -p) on C^ moreover, those studies were made for weirs and not for side weirs. 3.2 Verification with Additional Measurements In order to verify the obtained equation for C4>, additional measurements were carried out, in which (h1 -p) and F1 were fixed for all measured variants and only ф was varied for 7 different opening gate angles ф. In this way, it was possible to assess the effect of ф on Cф and remove the influence of other parameters. First, the coefficients Cd were obtained from the measurements. Then, using the value obtained for Cd at the angle of 90°, the coefficients Cф were calculated. All geometric and other measured hydraulic parameters from the conducted additional measurements are listed in Table 2. The agreement of the measured values of Cd and the calculated with Eqs. (5) and (7) was found to be good, as the maximum relative difference between them is 4 % for Cd and 2.8 % for Cr The calculated Cф according to Eq. (7), and the obtained Cф from the additional measurements are given in Fig. 10. Table 2. Geometric and other measured hydraulic parameters for the additional measurements while at angles from ф = 33° to ф = 90° it increases only slightly (less than 10 %). ф Ö! Ös P0 h - p L B F, [°] [l/s] [l/s] [cm] [cm] [cm] [cm] [-] 0 5.92 1.9 7.5 3.54 20 20 0.26 10 6.85 2.36 7.5 3.52 20 20 0.26 20 7.54 2.68 7.5 3.51 20 20 0.26 33 8.77 2.83 7.5 3.54 20 20 0.26 45 9.70 2.74 7.5 3.50 20 20 0.26 70 11.02 2.47 7.5 3.53 20 20 0.26 90 11.40 2.42 7.5 3.68 20 20 0.24 Fig. 10. Comparison of the values obtained from measurements and with the Eq. (7) calculated values of Сф From the comparison of the values obtained from the measurements and the values calculated from Eq. (7) for the examples of the additional measured variants, a correlation factor r2 = 0.99 was found. Therefore, we can conclude that Eqs. (5) and (7) give adequate results. 4 DISCUSSION OF THE OBTAINED Cd 4.1 Water Surface and Flow Area The lower values of Сф in the case when the flap gate is fully opened can be explained by the fact that such a position is similar to a broad-crested side weir, where the flow rate over the weir is smaller than by sharp-crested side weirs, due to a different contraction of the water jet. By closing the gates (by increasing ф), the contraction of the water jet changes and it becomes much more similar to the water jet that occurs at sharp-crested side weirs. An example of the different contractions of the water jet at sharp- and broad-crested side weirs is shown in Fig. 1. To confirm this hypothesis, the flow area of the cross-section at the edge of the flap gate was obtained from the additional measured variants, where the water levels were measured using a point gauge. Fig. 11 shows that the flow area A, in spite of the same inflow conditions (h1 -p) and F1, increases rapidly at angles from ф = 0° to ф = 33° (over 50 %>), 40 50 , states that the compressive stress does not contribute to crack onset. Whenever the damage initiation criterion of Eq. (3) for the cohesive element is satisfied, the stiffness of the element is declined according to the corresponding constitutive law that is illustrated in Fig. 1. 1.2 Damage Propagation Criterion After damage initiation, the softening procedure occurs. This procedure is governed by the corresponding traction separation law as below. 1.2.1 Bilinear Traction-Separation Model As shown in Fig. 1a, in the mixed-mode constitutive law, cmc and 5mc represent the cohesive strength and critical separation (separation at which crack initiates), respectively. Moreover, 5mu refers to the separation at which failure occurs. The softening relation of cohesive elements that are governed by bilinear constitutive law can be expressed as: a, = (1 - d)K ß,, (4) where d is a variable that relates damage condition, which has the magnitude d = 0 when the interface is undamaged, and the magnitude d = 1 when the interface is completely fractured. In numerical analysis of damage evolution, there are two crack evolution criteria: energy- and displacement-based. In the analysis presented herein, the energy-based Benzeggagh and Kenane (BK) [14] damage evolution criterion was used, as expressed in Eq. (5). G = Glc + (Gnc - Glc) G G (5) In Eq. (5), n refereed to the BK material parameter, GShear = GII +GIII and GT = GI + GII + GIII- 1.2.2 Modified PPR Model In this model, the traction force in the interface is obtained by differentiating a potential function with respect to the interface separation. Fig. 1b shows the typical modified PPR traction-separation law of cohesive zone modelling. Park et al. proposed this potential based model for mixed mode fracture [15]. This law can be used for a wide variety of ductile and brittle fractures. Since it behaves nonlinearly prior to damage, it is required to develop a user defined element in ABAQUS for this model. Bhattacharjee et al. [13] developed a modified version of the PPR model for analysis of tearing in thin soft materials. In their model, as with the triangular constitutive law, a linear elastic response was assumed before damage initiation (ò< òc). Therefore, this modified version of the PPR model allows for a straightforward implementation of the model in the commercial finite element code ABAQUS, using tabular capability. After damage initiation, the softening relationship in this model can be expressed as [13]: 8 4„,w 8 л w-i g = A[w(1 -f- T(- + )w-1 - 8u « 8u -a(1 -88 )a-1(w + 8 )w ], 8 a 5, (6) where ò is normal displacement, and A, w, a, òu are PPR parameters that can be determined by satisfying all boundary conditions, including: a = a at 8=8 , a = 0 at 8=8, lim ^ = o, 8^80 08 G = J ad8, (7) where G is the total strain energy release rate. By applying the above-mentioned boundary conditions, it can be expressed as: w = - a(a - 1)A2 1 -a!2 A = C ' C = [w(1 - !)a (W + !)w-1 - a (1 -!)a-1(W + X)w], (8) a a where k = 8, òu can be determined using the following equation as: G = - KSS + f adS. n c J (10) In numerical damage simulation, the corresponding the damage variables can be defined using Eq. (4) and the corresponding traction in Eq. (6). After generating the damage variable at all the displacements, the modified PPR model can be directly implemented in ABAQUS using the tabular capability as a function of relative displacement. 1.3 Cohesive Zone Length As shown in Fig. (2), the length of the cohesive zone lcz is introduced as the distance between the crack tip and the point where the maximum cohesive traction is achieved. Fig. 2. Cohesive zone length of DCB specimen The length of the cohesive zone at the initiation of crack growth is independent of applied load and crack length. This means that the cohesive zone length at the crack initiation is a material property. There are several models in the literature that estimated the length of the cohesive zone [4] to [6], [16] and [17]. All of the models for estimation of the cohesive zone length have similar forms as mentioned in Eq. (11). G' (11) L ■■RE-2-, where GIc is the critical strain energy release rate, oIc is the normal cohesive strength, and R is a parameter n 0 that depends on each cohesive zone model. For instance, Irwin's model [17] carried out with R = 0.31 or in Dugdale [4] and Barenblatt [5] models, the value for R is considered to be 0.4. 2 FINITE ELEMENT SIMULATION Three-dimensional finite element models are developed in ABAQUS 6.12 [18] to model the delamination onset and growth, using two different constitutive laws. The DCB model was composed of two layers of eight-nodded solid elements, C3D8R (adherent layers), which were connected with a layer of zero thickness eight-node cohesive element, COH3D8, (cohesive layer). Adherent layers were connected to the cohesive layer by surface-to-surface tie constraints. Fig. 3 illustrates the deformed shape of the DCB specimen during crack propagation. A more refined mesh near the crack tip, the outer surface of specimen and in the damage propagation region was used. Boundary conditions included applying a vertical displacement and horizontally restraining the upper and lower edge node of the arms (Fig. 3). In order to predict an accurate propagation of delamination, it is necessary to have an adequate fine finite element discretization in the cohesive zone length to achieve a good estimation of the interlaminar stress fields. When the cohesive zone is discretized by too few numbers of elements, the distribution of tractions ahead of the crack tip is not represented accurately. Thus, in order to achieve successful FEM results, it is necessary to have a minimum number of elements in the cohesive zone length. The number of cohesive elements in the cohesive zone is: minimum number of elements in the cohesive zone length. However, this number is not well established. N 1, l (12) where le is the length of the cohesive elements in the direction of crack propagation and lcz is the cohesive zone length. There are several guidelines for the Fig. 3. Deformed shape of simulated DCB with 3D elements 3 EXPERIMENTAL APPROACHES 3.1 Mode I Fracture Toughness Measurement 3.1.1 The DCB Test Fig. 4a shows the geometry and dimensions of the DCB specimens. In the process of specimen fabrication, the E-glass/epoxy composite plate with a thickness of 2h = 4.2 mm was first prepared. The E-glass fibres were impregnated with a ML506 epoxy resin and HA11 hardener. The laminates were fabricated with the hand lay-up technique, and the pre-crack length was produced by positioning a 13 цт thick Teflon insert at the mid-plane of the plate. In order to produce plates with the desired fibre volume fraction, special pressure tool was applied in order to squeeze out excess resin. Then, the plate was left at room temperature for 24 h to dry. After that, the plate was trimmed with a water jet machine along the longitudinal direction in order to obtain narrow specimens with the desired dimension and initial crack length. After trimming, the nominal length (l) and the nominal width (b) of the DCB specimen were 108 and 25 mm, respectively. The initial crack length (a) was 40 mm. Fibre volume fractions, Vf , measured using the resin burn-off method. Table 1 lists the mechanical property of E-glass/epoxy woven fabrication with Vf = 49% that was used in this research. All of the tests on DCB specimen were mm bj cjl Fig. 4. a) DCB specimen, b) typical DCB test, c) crack length measurement during propagation completed on a ZWICK electro-mechanical loading frame at room ambient temperature. Fig. 4b shows a typical DCB test. All DCB tests were carried out using displacement control at the crosshead speed of 1 mm/min according to ASTM D5528 standard [19]. Three specimens with a = 40 mm crack length were tested. A load-displacement response was recorded for each specimen during the test. In this work, the crack length was monitored by bonding a strip of paper with the graduations printed on it to the specimen's edge and by taking photos during the experiments with 5 s intervals using a 5 megapixel digital camera. The experimental magnitudes of P-S-a as a function of time were determined. The time of each P-S data point was computed using the applied displacement and the corresponding loading rate. The time for each value of crack length is the one at which the corresponding photo was taken. The photo in Fig. 4c was taken during a test and shows the crack tip, allowing the crack length measurement. These experimental results were then used to verify the adequacy of the three-dimensional finite element scheme utilized to obtain GIc. 3.1.2 Data reduction Method for GIc In order to calculate the mode-I critical strain energy release rate, there are many commonly used data reduction methods, including compliance calibration method (CCM) that is based on the Irwin-Kies principle, direct beam theory (DBT) and the Corrected Beam Theory. In the present work, the corrected beam theory proposed by de Moura et al. [20] was used. Only one specimen is sufficient to obtain the resistance curve (R-curve) of the specimen, which is the main advantage of the presented method. According to this theory, mode-I critical strain energy release rate can be computed as: GIC = 3PÖ 2b(a + |Д| )' (13) where P, a and b are load, crack length and specimen width, respectively. The parameter A is the crack length correction to account the crack tip rotation and deflection. According to the beam theory, the relationship between the compliance and the crack length is expressed by: C1/S = 2(a~ H ) (14) h(E1b)in ' The crack length correction can be obtained using the linear regression of C1/3 versus crack length data. 3.2 Normal Cohesive Strength Measurement The objective of the NCS test is to determine the Mode I cohesive strength (oIc) as an essential parameter for description of the traction separation law of cohesive elements. The fabrication of an NCS specimen is similar to that of a double cantilever beam except for the delaminated area. In the process of NCS specimen fabrication, first a 14-layer woven rectangular plate [0/90] 14 was produced. After the 7th layer of fabrication, a 13 p,m thick Teflon layer was inserted at all sides of the plate so that a 10 mm x 10 mm square area at the middle of plate was released. This area is the cohesive area of NCS specimen, as shown in Fig. 5a. After that, the plate was trimmed to a 80 x 50 rectangular dimension so that the cohesive area was located at the middle of rectangular specimen. The specimen was bonded to fixture surfaces. Prior to bonding, the surfaces of both fixture and specimen, were lightly roughened with the sandpaper on the bonding face and cleaned with acetone. The area of bond is large enough compared to the cohesive area so that debonding would not occur between the specimen and fixture surfaces the testing procedure. Fig. 5a shows the NCS specimen after testing, and Fig. 5b shows a typical NCS test. ajl b) Fig. 5. a) NCS specimen after test, b) NCS test As shown in Fig. 5a, the dimensions of rectangular plates B1 and B2 are 80 mm and 50 mm, respectively, and the width of cohesive the square area (C) is 10 mm. All of the NCS tests were carried out using displacement control at the crosshead speed of 0.5 mm/min until fracture. In this process, three NCS failure tests were completed. In order to prevent slippage during the test, the specimen was accurately balanced and a very little clamping force was required. Decohesion between the fibre and matrix in the NCS cohesive area is the dominant failure mode. Thus, it is obvious that the bulk matrix behaves differently than the thin cohesive layer due to constraint effects induced by the adhesion between the fibre and matrix. As a result, the bulk matrix properties could not be used to determine normal cohesive strength, and this parameter should be determined using an NCS test method. For the purposes of data reduction, all specimens were assumed to have failed instantly. The specimen failure is assumed at its maximum load value. The normal cohesive strength values (oIc) were computed using Eq. (15) as: P — max/ (15) where Pmax is the maximum load in which failure occurs and ANCS is the cohesive area of the NCS specimen. 4 RESULTS AND DISCUSSION 4.1 Experimental Results Fig. 6 shows the load-displacement response of three NCS experiments. The measured load is initially negligible, corresponding to the clearance of the fixture. The test was preceded until a maximum load was achieved and followed by a sudden load drop, indicating specimen failure. Displacement [mm] Fig. 6. Load-displacement response of NCS test The mean value of maximum loads of Fig. 6 was considered as Pmax. Using Eq. (15) for data reduction and substituting this value for Pmax, the cohesive strength of composite material computed and is equal to 12.42 MPa. To calculate the experimental resistance curve, the numerical values of PS-a parameters were recorded during crack propagation and were used to obtain the critical fracture energies versus crack length. Fig. 7 shows the experimentally obtained R-curve of the material. To simulate the crack propagation using cohesive elements, the mean value of GIc was considered as the fracture toughness of the material. Table 2 lists the cohesive properties of E-glass/epoxy woven composite laminate. Fig. 7. R-curve of E-glass/epoxy specimen Table 1. Mechanical properties of E-glass/epoxy E1 E2 E3 [GPa] [GPa] [GPa] Vi2 G12 [GPa] G13 [GPa] G23 [GPa] 18.43 18.43 3.57 0.15 2.85 2.85 2.85 Table 2. Interfacial property of E-glass/epoxy °Ic TIIc TIIIc [MPa] [MPa] [MPa] GIc [J/m2] GIIc [J/m2] GIIIc [J/m2] n 12.42 22.64 22.64 604 720 720 1.8 4.2 Cohesive Zone Model Results FEM models of each specimen were carried out using the Ply elastic properties of adherent layers that are given in Table 1 and the interfacial properties obtained previously and listed in Table 2. It should be mentioned that when using Eq. (2) for interfacial penalty stiffness, the value of K = 85 x106 N/mm3 is used for all DCB simulations. 4.2.1 Determination of Cohesive Zone Length Cohesive zone length was previously introduced as the distance between the maximum traction and the crack front. Therefore, in order to calculate the distribution of traction in the cohesive layer of model and the corresponding cohesive zone length, a very refined mesh using element size of 0.125 mm was used in the area near the crack tip of the DCB specimen. The distribution of tractions ahead of the crack tip at the delamination onset of the DCB specimen is illustrated in Fig. 8. At the crack initiation point, traction at the crack tip vanished as expected from the cohesive zone theory. According to this analysis, the cohesive zone length of material is 3.76 mm, as shown in Fig. 8. Using the material property that shown in Table 1 and 2 and Eq. (11), the parameter R is computed as 0.27. This value is closest to the Irwin [17] (0.31) model. The cohesive zone length is a material property that has a high order of importance regarding obtaining a successful prediction of delamination onset. This parameter was previously introduced as a function of normal cohesive strength by Turon et al. [11]. Thus, using an absolute value for normal cohesive strength, this parameter as a material property is determined. due to the magnitude of tractions ahead of the crack tip. J-1 Distance from crack tip [mm] Fig. 8. Distribution of traction ahead of crack tip 4.2.2 Investigation of Mesh Refinement In order to investigate the effect of mesh refinement in the cohesive zone length on numerical prediction of delamination onset, several DCB specimens were simulated with different sizes of cohesive elements in the cohesive zone length. The predicted load-displacement responses obtained using DCB models are compared to the experimental solution in Fig. 9. In this analysis, cohesive element sizes (in the direction of crack propagation) range from 0.125 mm to 2 mm. The results illustrate that for all mesh sizes a converged solution was obtained but it is necessary to apply a mesh size, le, less than 1 mm to accurately predict delamination initiation. The analysis using coarser meshes significantly overpredicts the strength of the DCB specimen, and the response does not follow the experimental results. Cohesive zone length, lcz, for the material described in Tables 1 and 2 was determined as 3.76 mm. Therefore, for a mesh size smaller than 1 mm, more than three elements would span the cohesive zone length, which is enough for an accurate prediction of the fracture onset. There are several guidelines for the minimum required element in cohesive zone length. For example, Moes and Belytschko [21] proposed using more than 10 elements, while Falk et al. [16] suggested between 2 and 5 elements. In the work of Dàvila et al. [22], the minimum required element length to predict the delamination in a DCB model was 1 mm, and using more than 3 elements in cohesive zone length of simulated DCB specimen was recommended. The difference in predictions from using a coarse mesh in the modelling of delamination in a DCB specimen is 4 6 Displacement [nun] Fig. 9. Damage simulation using different mesh refinement 4.2.3 Comparison of Damage Models A study also was conducted to investigate the adequacy of the two mentioned traction separation laws used to simulate damage propagation. The objective was to determine how the used models reproduce crack initiation and propagation. Fig. 10 shows the load-displacement of the cohesive zone model and experimental work on a DCB specimen. Fig. 10. Load displacement response of experimental and damage models In the current study, the value of shape parameter a in the PPR model varied from 1 to 4; the result was that increasing the value of the shape parameter a increased the rate at which material loses its stiffness once damage was initiated. In other words, increasing parameter a decreases the fracture process zone effect ahead of the crack tip. In the case of a < 2, there is a gradual fall in the load, but in the case of a > 2 a sudden drop in the load-displacement response is achieved, which means a large number of cohesive elements failed at the same time. For clarity, the results are not shown in this figure. In this study, a value of a = 1.7 was found to be the optimum value for the numerical prediction of damage propagation. As shown in Fig. 10, the modified PPR model was found to be adequate to reproduce the experimentally observed behaviour of the composite specimen, and reproduced approximate smooth crack initiation and propagation while the bilinear model depicted sudden damage propagation. The maximum difference between the experimental and bilinear models is 8.8 % while for PPR it is 2.6 %. This means the modified PPR model accounted fracture process zone which created ahead of crack tip. 5 CONCLUSION A methodology for the delamination characterization of composite laminates under pure Mode I was proposed. • An NCS test has been proposed to compute the Mode-I cohesive strength as a cohesive parameter. • The Mode-I critical strain energy release rate versus the crack length of E-glass/epoxy composite laminate was computed using corrected beam theory for data reduction. • A mixed-mode triangular constitutive relationship between stress (c) and relative displacements (5) of cohesive elements and modified PPR damage model were considered to simulate delamination onset and propagation. • The results of the three-dimensional finite element analysis with cohesive parameters (cIc, GIc) enclosed the adequacy of cohesive parameters. • Accurate damage prediction was achieved using the modified PPR model, and it was considered by the authors to be an accurate model for damage characterization of material. • Modified PPR models accurately described the fracture process zone, which was created ahead of the crack tip as compared to bilinear model. • To ensure the sufficient dissipation of energy, cohesive zone length as a material property was determined. • Numerical analysis with different discretizations of the cohesive zone length showed that numerical predicted responses correlate well with the experimental solutions when at least 3 elements span the cohesive zone length. 6 ACKNOWLEDGMENT This work has been funded by University of Tabriz, and its authors would like to thank the University of Tabriz for the grant. 7 REFERENCES [1] Karpiliski, A. (2006). An Introduction to the diagnosis of the delamination process for glass/epoxy com posites during high-pressure abrasive water-jet cutting. Strojniški vestnik -Journal of Mechanical Engineering, vol. 52, no. 7, p. 532-538. [2] Hasiotis, T., Badogiannis, E., Tsouvalis, N.G. (2011). Application of ultrasonic C-scan techniques for tracing defects in laminated composite materials. Strojniški vestnik - Journal of Mechanical Engineering, vol. 57, no. 3, p. 192-203, D0l:10.5545/sv-jme.2010.170. [3] Brunner, A., Blackman, B., Davies, P. (2008). A status report on delamination resistance testing of polymer-matrix composites. Engineering Fracture Mechanics, vol. 75, no. 9, 2779-2794, DOI:10.1016/j.engfracmech.2007.03.012. [4] Dugdale, D. (1960). Yielding of steel sheets containing slits. Journal of the Mechanics and Physics of Solids, vol. 8, no. 2, p. 100-104, DOI:10.1016/0022-5096(60)90013-2. [5] Barenblatt, G. (1962). The mathematical theory of equilibrium cracks in brittle fracture. Advances in Applied Mechanics, vol. 7, no. 1, p. 55-129, DOI: 10.1016/j.ech.2007.03.012. [6] Hillerborg, A., Modéer, M., Petersson, P.-E. (1976). Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research, vol. 6, no. 6, p. 773-781, DOI:10.1016/0008-8846(76)90007-7. [7] Song, S.H., Paulino, G.H., Buttlar, W.G. (2006). Simulation of crack propagation in asphalt concrete using an intrinsic cohesive zone model. Journal of Engineering Mechanics, vol. 132, no. 11, p. 1215-1223, DOI:10.1061/(ASCE)0733-9399(2006)132:11(1215). [8] Zhang, Z.J., Paulino, G.H. (2005). Cohesive zone modelling of dynamic failure in homogeneous and functionally graded materials. International Journal of Plasticity, vol. 21, no. 6, p. 1195-1254, DOI:10.1016/j.ijplas.2004.06.009. [9] Khoshravan, M.R., Moslemi, M. (2014). Investigation on mode III interlaminar fracture of glass/epoxy laminates using a modified split cantilever beam test. Engineering Fracture Mechanics, vol. 127, p. 267-279, DOI:10.1016/j. engfracmech.2014.06.013. [10] Song, S.H., Paulino, G.H., Buttlar, W.G. (2006). A bilinear cohesive zone model tailored for fracture of asphalt concrete considering viscoelastic bulk material. Engineering Fracture Mechanics, vol. 73, no. 18, p. 2829-2848, DOI:10.1016/j. engfracmech.2006.04.030. [11] Turon, A., Davila, C.G., Camanho, P.P., Costa, J. (2007). An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Engineering Fracture Mechanics, vol. 74, no. 10, p. 1665-1682, DOI:10.1016/j.engfracmech.2006.08.025. [12] Yuan, H., Li, X. (2014). Effects of the cohesive law on ductile crack propagation simulation by using cohesive zone models. Engineering Fracture Mechanics, vol. 126, p. 1-11, DOI:10.1016/j.engfracmech.2014.04.019. [13] Bhattacharjee, T., Barlingay, M., Tasneem, H., Roan, E., Vemaganti, K. (2013). Cohesive zone modelling of mode I tearing in thin soft materials. Journal of the Mechanical Behavior of Biomedical Materials, vol. 28, p. 37-46, D0l:10.1016/j.jmbbm.2013.07.015. [14] Benzeggagh, M., Kenane, M. (1996). Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Composites Science and Technology, vol. 56, no. 4, p. 439449, D0I:10.1016/0266-3538(96)00005-X. [15] Park, K., Paulino, G.H., Roesler, J.R. (2009). A unified potential-based cohesive model of mixed-mode fracture. Journal of the Mechanics and Physics of Solids, vol. 57, no. 6, p. 891-908, DOI:10.1016/j.jmps.2008.10.003. [16] Falk, M.L., Needleman, A., Rice, J.R. (2001). A critical evaluation of cohesive zone models of dynamic fracture. Le Journal de Physique IV, 11(PR5), Pr5-43-Pr45-50. [17] Irwin, G. (1997). Plastic zone near a crack and fracture toughness. [18] Abaqus/CAE user manual. Abaqus 6.12 Documentation. [19] ASTM D5528:2007. Standard Test Method for Mode I Interlaminar Fracture Toughness of Unidirectional Fiber-Reinforced Polymer Matrix Composites, ASTM International West Conshohocken. [20] De Moura, M., Campilho, R., Gongalves, J. (2008). Crack equivalent concept applied to the fracture characterization of bonded joints under pure mode I loading. Composites Science and Technology, vol. 68, no. 10-11, p. 2224-2230, DOI:10.1016/j.compscitech.2008.04.003. [21] Moes, N., Belytschko, T. (2002). Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics, vol. 69, no. 7, p. 813-833, DOI:10.1016/S0013-7944(01)00128-X. [22] Davila, C.G., Camanho, P.P., de Moura, M. F. (2001). Mixed-mode decohesion elements for analyses of progressive delamination. Proceedings of the 42nd AIAA/ASME/ASC, DOI:10.2514/6.2001-1486. Strojniški vestnik - Journal of Mechanical Engineering 61(2015)9, 517-522 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2553 Original Scientific Paper Modelling of an Electrohydraulic Proportional Valve with a Synchronous Motor Andrzej Milecki* - Dominik Rybarczyk Poznan University of Technology, Institute of Mechanical Technology, Poland This paper presents the design of proportional valve with a Permanent Magnet Synchronous Motor (PMSM). The proposed valve is described, and its history is briefly reviewed. Basic equations describing the valve are formulated and its simulation model is implemented in MATLAB-Simulink software. Selected nonlinearities are included in this model. In order to determine the basic parameters of the discussed proportional valve, a test stand is built, which enables valve investigations. In this test stand, a valve control system based on programmable logic controller (PLC) with a touch panel and inverter module is implemented and used for investigations. The valve flow characteristics and step responses obtained in simulations are presented. These characteristics are compared with results obtained in experimental investigations. As a result, the valve simulation model is modified and improved. Keywords: proportional valve, electrohydraulic, synchronous motor Highlights • Proposing a new type of hydraulic proportional valve in which the synchronous motor is used. • Proposing a theoretical description and simulation model of the valve with a synchronous motor. • Presenting the proposed laboratory investigation results of the valve. • Improving the simulation model of the valve based on the outcomes of the investigations. 0 INTRODUCTION Electrohydraulic servo-drives can be controlled by two types of electro-valves: servo- and proportional valves. The servo-valves are used in high-accuracy applications. Proportional valves respond acceptably to the requirements stated in most industrial cases and are much cheaper than servo-valves. Therefore, these elements are commonly used in many industrial electrohydraulic drives. Proportional valves are activated by solenoids and can provide a smooth and continuous variation in flow or pressure in response to an electrical input valve. The basics of the design of proportional valves were established approximately 30 years ago and no significant progress in this area has been visible since then [1] and [2]. Many investigations and related publications about electrohydraulic servo drives have been focused on the improvement of the properties of these drives, via the implementation of modern forms and types of control [3]. Only a few works deal with improvements of the parameters of electrohydraulic drive elements, as well as with finding new ways to provide very accurate movement of valve parts, such as a valve's spool. Thus far, in proportional valves, proportional electromagnets are commonly used as the slider driving part [4]. Only in a few solutions can stepping or DC motors be found, even though they have been known for many years. Based on that, we decided to apply a new type of drive device, a permanent magnet synchronous motor (PMSM) as electromechanical actuator in the proportional valve (this motor is characterized by very good properties and can assure very good positioning accuracy). Over the previous decade, only a few papers focused on the applications of different types of motors in hydraulic proportional valves have been published. Murrenhoff [5] described the cross-cutting trends in the design and the development of electrohydraulic valves. He presented an interesting solution, using the direct drive in a proportional valve. In another solution, the use of a stepping motor to transfer the ball screw, which was moving the mechanism with four independent flow sliders, was proposed. The design and investigations of valve in which a piezo actuator was used were described in [6]. Such a valve is characterized by very good dynamic parameters in comparison with valves controlled by proportional electromagnets. Boes at al. [7] pointed out the advantages of valves with integrated control electronics. The use of powerful processors and a decentralized control system allowed the easy integration of electrohydraulic valve with other elements. Myszkowski and Milecki [8] described the use of a stepping motor in the valve to obtain a very low velocity of the electrohydraulic drive. It is worth emphasizing that the drive was able to move with a very low constant speed equal to 1 mm/s. However, the drive maximum speed was only 0.125 mm/s, which was rather low and reduced the number of possible applications. In another study [9], an unconventional electro-hydraulic proportional flow control valve based on a switching solenoid and a fuzzy-logic controller was proposed. The switching solenoid's non-linear force/ stroke characteristic was linearized by a fuzzy-logic controller. The obtained laboratory investigations of this valve were quite satisfactory, assuring low cost. Šimic et al. [10] presented a new approach to the modelling and simulation of hydraulic spool valves by using simple mathematical expressions to describe the geometry of the sliding spool metering edge. Consequently, different shapes of spool-metering edges in combination with other functional elements could be used in design of hydraulic valves. The literature overview has shown that only in one paper [11] was the use of an electric servo motor for control of the proportional valve described. However, this solution was significantly different from commonly used designs. The servo-motor that was used in a valve was able to produce torque of as much as 15 Nm. The cam mechanism used allowed the motor shaft to move in the range of only a few degrees of arc. The maximum valve flow was as much as 400 dmVmin; therefore that valve was much bigger than the one considered in this paper. 1 VALVE CONSTRUCTION 7 1 4 5 6 Fig. 1. Scheme of proportional valve with synchronous motor Fig. 2. View of proportional valve with PMSM The rotor position is measured by an absolute encoder type EnDat 2.2 (7), providing a current position, even after a power failure, with a very high positioning accuracy, equal to 262144 pulses on revolution. As a result, the drive is able to assure the linear resolution of 0.5 mm [12] and [13]. The photo of the valve is presented in Fig. 2. The valve is connected to the cylinder by a specially designed manifold. The schematic of the valve proposed in this paper is presented in Fig. 1. The valve spool (3) is actuated by a permanent magnet synchronous motor (1). In the proposed valve the motor type 8LVA22 (B&R company) is used. Its basic parameters are power 105 W, current 2.9 A, nominal torque: 0.65 Nm. This motor is connected to the spool (3) by a flexible bellows coupling (2). The second ending of the spool is directly connected to a ball screw (5), whose nut (6) is fixed to the valve body. The rotation of the motor causes the rotation of the ball screw and the axial movement of the spool and its control edges. This movement is proportional to the angular motor displacement and the pitch of the thread used. The direction of rotation determines the direction of spool translation and the opening or closing of valve gaps. It results in the flow of oil to and from the actuator chambers and the displacement of the piston. The spool diameter is 10 mm. There are three rectangular gaps in the body with dimensions of 2.5 mm x 2 mm. The valve size is 10 (below 64 dmVmin). 2 VALVE MODELING AND BASIC INVESTIGATIONS In order to investigate the designed and built valve, a test stand was arranged. Its schematic diagram is presented in Fig. 3. A valve control system is based on PLC working under a real time operating system Automation Runtime with a Power Panel 500 touch panel and ACOPOS servo-controller used for PMSM control [14]. A HySense QG100 flow meter is used, whose parameters are max. flow Q = 30 dmVmin, max. pressure p = 30 MPa, non-linearity ± 0.5 %, 1640 pulses per dm3. It is connected to the counter module in the PLC unit. As a communication interface between the PLC and the servo controller, a Powerlink interface is used. The view of the control system is shown in Fig. 4. In the synchronous motor controller, two feedback loops are used. The inner one is a velocity feedback in which a P type regulator is used. The second one is the positioning loop with a PI regulator. Their parameters are set using autotuning mode, based on the detection of one or two points of the process frequency response using relay excitation described by Äström and Hägglund [15]. The following values of parameters are used: kv = 0.4 1/s for the velocity loop and kP = 220, integrative time T = 0.0004 s for positioning loop. Fig. 3. The scheme of control stand The valve flow through the sharp-edged nozzles is turbulent [2] and [9], which can be described using the following general equation: Q = KQy[Äp • X, (1) where Kq is a flow coefficient, Ap pressure drop on a valve, and x a spool displacement (±2.5 mm). Fig. 4. A view of proportional valve with a synchronous motor The flow coefficient Kq depends on the following parameters: where is the discharge coefficient, p oil density: 780 dm3/min, w width of the gap (6 mm). The value of the discharge coefficient depends mainly on the geometry of the slot. In the case of the proportional valve considered in this paper, this value is usually taken as 0.64 [2]. Using Eq. (1), the theoretical calculations of valve flow Q for spool displacements in a range of ±2.5 mm and for pressures 4 MPa, 6 MPa, 12 MPa and 14 MPa are made and shown in Fig. 5. In order to verify these calculations, the real valve is investigated. A flow Q changes for two supply pressures, 6 MPa and 12 MPa, are measured. Due to the applied flow meter, it was possible to measure flow only until 30 dm3/min. It can be noted that the characteristics of a real valve are curved when the input signal increases. Most probably, this is caused by the valve spool's non-symmetry. In the test stand (Fig. 3), both gaps are connected serially. If one of these gaps is smaller than the other, the pressure drop on the first one is bigger. As a result, the flow saturation occurs in this gap, which is determined by pressure square root characteristic, described by Eq. (1). This is visible in the presented laboratory results. -1.25 0 1.25 Fig. 5. The flow characteristics Kq 2/P • W (2) Fig. 6. A four-edge hydraulic amplifier The hydraulic part of proportional valve is a four-edge amplifier, shown in Fig. 6. Its static behaviours can be characterized using so-called square root equations, which describe the oil flow through the hydraulic nozzles. These flows can be described by following equations: _ . Qa it ) = Кел/Ро - p a it ) ■ x(t ) for x > 0: ,__, (3) Qbit) = Ke4P() ■ x(t) Qa (t) = Ke4PJi) - x(t) for x < 0 :_ , (4) Qb(t) = Keylp0 - Pb(t) - x(t) where s is a spool displacement [mm], Qa flow through the gap A [dm3/min], Qb flow through the gap B [dm3/min], Kq flow coefficient, p0 supply pressure [Pa], andpa,pb pressures in chambers A and B [Pa]. The spool dynamics can be described by the following PMSM motion equation: d2x(t) ~dx(t) 2n „ . ^ . m-— + D—— = — KTi(t ) = Kj(t ), (5) dt dt P where m is the mass reduced in motor axis (0.2 kg), D viscous friction coefficient (120 Ns/m), P ball screw pitch (2.5 mm/rev.), T PMS motor torque, KT torque coefficient (0.23 Nm/A). In the first approach, the electric circuit of PMS motor is described by the second order transfer function: L^ + Ri(t ) = U (t), dt (6) where L is 4.1 mH, R = 2 Ohm coil inductance and resistance, i current, U supply voltage. After the Laplace transformation of Eqs. (5) and (6), and their combination we obtain: * < •> = f^-^I U (s ). Ts + s t2s + I The drive transfer function is as follows: G(s) = k X is) =_ U (s) s(Tls + 1)(T2 s +1) ' (7) (8) where T1 = m/D and T2 = L/R are time constants. Using Eqs. (3), (4) and (7), a simulation model of the proportional valve with PMSM is prepared in Matlab-Simulink software (Fig. 7). It includes square root flow non-linearity, which is modelled using the function Sqrt. A simplified model of the PMSM is also included. The rotor position and the valve spool position are simultaneously measured by element characterized by coefficient km and fed back to the PI type, position controller. The measured position is derived in order to obtain a velocity signal, used in inner feedback loop. a) b) Fig. 7. Matlab-Simulink simulation models of the: a) valve, b) PMSM (simplified) 2.5 2 1.5 -0 -2.5 0 5 10 15 20 25 t [ms] Fig. 8. The step responses of the valve spool 30 The valve step responses obtained by simulation (spool displacement S) are shown in Fig. 8. The time constant Tdom, which represents the time it takes the system's step to reach 0.63 of its final value, is equal to about 11 ms. 3 MODEL IMPROVEMENT In order to verify the correctness of the proposed simulation model, the real valve step responses are measured. The results are shown in Fig. 9. The recorded curves indicate that the valve is characterized by a significant time delay, which is caused by a PLC controller and communication interface, which is used to control the PMSM. Curves collected during several experimental investigations showed that the delay time constant T0 is equal to 10 ms. This delay is included into the following transfer function: 0 1 G (s) = X(s) = kd • esT° U (s) (T, • s +1) • (T2 • s +1)' (9) In this transfer function, the transport delay T0 and time constants T and T2 are used. The valve of dominant time constant Tdom depends on the assumed step signal values and varies from 7 to 1г ms (Fig. 9). These changes are most probably caused by a maximum current limitation of the motor. This current saturation is included in the model, shown in Fig. 10. a a 0 -1 Y Time delay T0 Dominant time constant: Tdom -2 - 0 5 10 15 20 25 Fig. 9. Step responsesofthe valve spool 30 U CH 1 ■ T T Saturation of Transport currnt slope T0 r Ti +T2 f T T T T2 Fig. 10. Non-linear model of the PMS motor 2 1 ^ A J= 0 * -1 -2 5 10 15 20 25 30 t [ms] Fig. 11. The model and real valve step response In Fig. 11, the step responses obtained in simulation with the use of model shown in Fig. 10 are presented. For comparison, curves obtained in experimental investigations are also shown in this figure. The maximum difference between the model and the real object curves emax is about 9 %. Consequently, one may conclude that the simulation results fit quite well with the results obtained in laboratory tests. The improved simulation model can then be regarded as suitable for the modelling and simulation of the proportional valve with PMSM drive. 4 CONCLUSIONS In this paper, the theoretical description of a proportional valve with a PM synchronous motor is presented. Based on this description, a computer simulation model of a proposed valve with PMSM is proposed and implemented in MATLAB-Simulink software. This model is then compared with a real valve and improved. The study includes the laboratory examination of basic characteristics, such as valve flow and step response. The chosen simulation results are shown and discussed. The laboratory investigations of the proposed valve with PMSM are presented and compared with the simulation results. The simulation model proposed in the paper fits reflected the real situation very well. The proportional valve with a permanent magnet synchronous motor that is presented in this paper represents an interesting alternative to servo-valves. It may assure better properties (especially dynamic and accuracy ones) than electrohydraulic standard proportional valves with solenoids do. 5 ACKNOWLEDGMENT This paper was supported by the Polish Ministry of Science and Education, grant no. 02/23/DSPB/1208, Poland. 6 REFERENCES [1] Chapple, P. (2003). Principles of Hydraulic System Design. Coxmoor Publishing Company, Oxford. [2] Cundiff, S.J. (2000). Fluid Power Circuits and Control Fundamental and Applications. CRC Press, Boca Raton, London, New York, Washington. [3] Ming, X., Jin, B., Chen, G., Ni, J. (2013). Speed-Control of Energy Regulation Based Variable-Speed Electrohydraulic Drive. Strojniški vestnik - Journal of Mechanical Engineering, vol. 59, no. 7-8, p. 433-442, D0l:10.5545/sv-jme.2012.911. X [4] Muraru, V., Muraru, C. (2000). Optimization of the proportional solenoids for electrohydraulic control systems. Proceedings of the 1st FPNI-PhD Symposium, p. 157-166. [5] Murrenhoff, H. (2003). Trends in valve development. Ölhydraulik und Pneumatik, vol. 46, no. 4, p. 1-36. [6] Herakovič, N. (1995). Piezoactuator for a single-stage servovalve with high dynamic response (Piezoaktorbetatigung fur ein einstufiges hochdynamisches Servoventil), Olhydraulik & Pneumatik, vol. 39, no. 8, p. 601-605. (in German) [7] Boes, Ch., Lenz, W., Mueller, J. (2003). Digital servo valves with fieldbus interface in closed loop applications. The 8th Scandinavian International Conference on Fluid Power, Tampere. [8] Myszkowski, A., Milecki, A. (2009). Modelling of electrohydraulic servo drive used in very low velocity applications. International Journal of Modelling, Identification and Control, vol. 7, no. 3, p. 246-254, D0l:10.1504/ IJMIC.2009.027211. [9] Renn, J.C., Tsai, C. (2005). Development of an unconventional electro-hydraulic proportional valve with fuzzy-logic controller for hydraulic presses. The International Journal of Advanced Manufacturing Technology, vol. 26, no. 1-2, p. 10-16, D0I:10.1007/s00170-003-1973-7. [10] Šimic, M., Debevec, M., Herakovič, N. (2014). Modelling of hydraulic spool-valves with specially designed metering edges. Strojniški vestnik - Journal of Mechanical Engineering, vol. 60, no. 2, p. 77-83, D0I:10.5545/sv-jme.2013.1104. [11] Wiegandt, W. (2010). Development of a servomotor driven proportional valve. 7th International Fluid Conference, Aachen. [12] Rybarczyk, D., Milecki, A. (2015). Modeling and control of proportional valve with synchronous motor. Solid State Phenomena, vol. 220-221, p. 457-462, DOI:10.4028/www. scientific.net/SSP.220-221.457. [13] SQdziak D., Regulski R. (2015). Design and investigations into piezobender controlled servovalve. Solid State Phenomena, vol. 220-221, p. 520-526, D0I:10.4028/www.scientific.net/ SSP.220-221.520. [14] B&R, Motion control, from http://www.br-automation.com/en-gb/products/motion-control/, accessed on 2015-02-20. [15] Aström K., Hägglund T. (1995). PID controllers: Theory, Design and Tuning. Instrument Society of America, Research Triangle Park. Strojniški vestnik - Journal of Mechanical Engineering 61(2015)9, 523-532 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2627 Original Scientific Paper Computerised Design and Simulation of Meshing and Contact of Formate Hypoid Gears Generated with a Duplex Helical Method Yu Zhang1 - Hongzhi Yan1* - Tao Zeng12 1 Central South University, State Key Laboratory of High-Performance Complex Manufacturing, China 2 Changsha Haliang Kaishuai Precision Machinery Co., China The duplex helical method has higher machining efficiency for face-milling spiral bevel and hypoid gears. An accurate and practical approach for calculating the basic machine-tool settings of spiral bevel and hypoid gears manufactured with the duplex helical method is proposed in the present work. The gear tooth surface vector functions and curvature parameters based on basic machine, and head-cutter settings are calculated. Two types of curvature parameters for pinion tooth surfaces are obtained by utilizing the conjugation of gear and pinion tooth surfaces, and the conjugation of pinion tooth surfaces and pinion head-cutter surfaces at the reference point. Next, the basic machine settings for generating pinions in accordance with the two types of curvature parameters for pinion tooth surfaces are determined. Finally, a numerical examples using the duplex helical method are performed and validated by comparing with experimental results. New hypoid gear software has been developed using the new approach. Keywords: duplex helical method, spiral bevel and hypoid gears, face milling, basic machine settings, curvature parameter Highlights • A new exact calculation approach for basic machine-tool settings is proposed. • Three reference points are used to determine the optimal machine settings. • Based on the optimal approach tooth bearings and functions of transmission error on the both sides are favorable. • Hypoid gear design software has been developed using the new approach. 0 INTRODUCTION Spiral bevel and hypoid gear drives are widely employed as transmission elements in vehicles, aircrafts, ships, and other gear reducers. There are currently two main methods of producing spiral bevel [1] and hypoid gears [2] in the production environment: the single indexing method referred to as "face milling" [3], and the continuous indexing method referred to as "face hobbing" [4]. In both face milling and face hobbing, the gear may be cut using either a generating method or a non-generating (formate) method [4]. However, the pinion of a pair of mating hypoid gears is always cut using the generating method to satisfy the required contact characteristics. The non-generating method offers higher productivity than the generating method because the generating roll is eliminated in the former method. The manufacturing of face-milled spiral bevel and hypoid gear sets can be accomplished by using the five-cut process [5] or by using the duplex helical method or completing process [6]. The five-cut process consists of five independent operations: two operations to finish the gear, and three operations to finish the pinion [5]. Generalized theory and application of bevel and hypoid gears generated by the five-cut process have been comprehensively presented by several gear scientists [7] to [9]. The generating method and the formate method of the five-cut process for face-milled spiral bevel and hypoid gears have been described in detail in [10] and [11]. Litvin et al. developed the principle and the calculation processes for the five-cut process independently in a manner that is different from Gleason's technology described in [12] and [13]. Astoul et al. [14] presented a new design method of spiral bevel gears based on an optimization process to reduce their quasi-static transmission error. Cao et al. [15] proposed a new method to design pinion machine tool-settings for spiral bevel gears by controlling contact path and transmission errors based on the satisfaction of contact the condition of three given control points on the tooth surface. Su et al. [16] proposed a new approach to designing and implementing a seventh-order polynomial function of transmission error for spiral bevel gears in order to reduce the running vibration and noise of gear drive and improve the loaded distribution of the tooth. Computerised design, manufacturing and simulation of meshing, and contact stress analysis of spiral bevel and hypoid gears are the subjects of research performed by many scientists [17] and [18]. Lin et al. [19] developed a numerical tooth contact analysis technique for simulating the single flank test of the gear geometry data measured on a gear measuring centre. Fan developed a new generalized tooth surface generation algorithm and a tooth contact analysis (TCA) approach [3], and presented a new generic model of generating spiral bevel and hypoid gears, this model is applicable to both face-milling and face-hobbing processes based on the universal motion concept (UMC) [20]. Simon [21] presented a new method for computer-aided tooth contact analysis in mismatched spiral bevel gears. In addition, Tamizharasan and Senthil Kumar [22] proposed an attempt to minimize flank wear of uncoated carbide inserts while machining AISI 1045 steel by finite element analysis, this simulation method can provide a reference for the finite element simulation of spiral bevel and hypoid gears. The principle and machining character of the duplex helical method are evidently different from that of the five-cut process. The cutters used for the duplex helical method have alternate (inside and outside) blades. The head cutter is mounted on the cradle that has a helical motion with respect to the gear blank; the work spindle is mounted on the sliding base that have an infeed motion with the rotation of the cradle. When a single cutter is used in one operation, both sides of the tooth slot are finished from a solid blank during machining. The advantages of using this method are as follow [6]: (1) the higher machining efficiency, (2) the assurance of uniform gears and, therefore, greater accuracy, since the size of the teeth is not dependent upon the manual controls of the operator, (3) the reduction of spoilage by manual mistake, and (4) the smoother bending of the bottom and sides of the teeth. The duplex helical method was invented several decades earlier by Gleason [5], although only some formulas and calculating instructions for the duplex helical method for hypoid gears have been published. However, as Gleason's technology is confidential, the public has little knowledge of its principle and method in detail thus far. To the best knowledge of the authors of this paper, Gonzalez-Perez et al. dealt with conversion of the specific machine-tool settings of a given hypoid generator to the neutral machinetool settings and adjustment the contact pattern by considering parabolic profiles on the blades of the head-cutter [23]. Fong proposed a mathematical model of a universal hypoid generator and applied it to simulate virtually all primary spiral bevel and hypoid cutting methods, including the duplex helical method, the supplemental kinematic flank correction motions, such as modified generating roll ratio, helical motion, and cutter tilt were included in the proposed mathematical model [24]. In this paper, the authors present a new method that is used to accurately calculate basic machine-tool settings for formate hypoid gears. The new method aims to: (i) present the generalized theory of the duplex helical method in detail, (ii) obtain the desired meshing quality for the duplex helical method by using precise calculation. 1 CALCULATION OF BASIC MACHINE-TOOL SETTINGS FOR GENERATING GEAR The coordinate system Sm{Xm, Ym, Zm} is rigidly connected to the cutting machine (Fig. 1). The top, bottom (A-A) and right (B-B) of Fig. 1 are the machine front view, the machine bottom view and the side view (the projection of the head cutter). The cradle rotates about the Ym axis; the p axis and g axis are projections of gear and pinion axes in the XmOmYm plane, respectively. The points M, Mj and O0 are the reference points of the tooth surface and the projection of M on the cutting edge and the centre of the head cutter, respectively. O2 is the cross point of the gear, and Om is the machine centre. The process for calculation of the machine-tool settings of the duplex helical method is the same as that of the five-cut process, details of which can be found in [10]. The unit vectors of the Ym axis, g axis, and p axis can be represented in the coordinate system Sm by the following equations: e = [0 1 (1) g = [- C0SYm 2 0 Sin Ym 2], (2) p = [- cos (X-rm2 ) 0 - sin (l-ym2 ). (3) The unit normal n0, the unit vector t0 and the position vector a0 of a point on the outside blade edge of the gear to the head cutter-generating surface (drive or convex side) can be defined in the coordinate system Sm by the following equations: n0 =[-cosa21sin ß02 - cosa21cos ß02 - sina21],(4) 10 =[-sina21smß02 -sina21cosß02 cosa21], (5) , (u ) = ut0 + H - Xg cosYm2 + rCG1Sin ß02 rcG1 C0S ß02 - V XG sin Ym 2 (6) Here, u is a profile parameter. The position vector of point M1 on the blade edge can be represented by a0(5G1) when u is equal to sG1 as denoted in Fig. 1. 2 CALCULATION OF BASIC MACHINE-TOOL SETTINGS FOR GENERATING PINION 2.1 Calculation of Pinion Curvature Parameters Based on Gear Tooth Surfaces To facilitate the description of this paper, the following definitions are made: E2 is the gear tooth surface. E1 is the pinion tooth surface. is the gear head-cutter surface or the generating tooth surface of the gear. is the pinion head-cutter surface or the generating tooth surface of the pinion. M1 is the reference point of the drive side. M2 is the reference point of the coast side. The pinion is the driving wheel and left hand, gear is the driven wheel and right hand. The driving side is the convex side of gear, the concave side of pinion, the inside blades of the gear head-cutter and the outside blades of the pinion head-cutter. The coast side is the concave side of the gear, the convex side of the pinion, the outside blades of the gear head-cutter and the inside blades of the pinion head-cutter. The formate-cut gear tooth surface is a copy of the surface of the head-cutter, which is a surface of revolution. Therefore, the vectors a0(sG1), t0, n0 of generating gear are the same as the vectors of the gear tooth surface. A rotation angle в1 of the gear about the g-axis is necessary for tangency at M1 of the gear and pinion tooth surfaces E2 and E1. The position vector of the reference point (M1) r1dr(61), the unit normal n1dr(61) and the unit vector t1dr(01) in the pinion tooth surface can be represented as: rWr (0,) = Ee + a1(01), a1(01) = a0 (^ )• R [g,0! ], «1* (A ) = n0 • R [gA ], t^ (0! ) = t о • R [g,ei ]. (7) (8) (9) (10= Here, R [g, 01] is a transformation matrix that denotes the rotation angle в1 about the vector g. Unit vector e is given in Eq. (1). The vectors aj(0j), e, g, p and are determined in the coordinate system S2{tWrxnWr, tw„ nWr} by the following equations: a2 (01 ) = a, (Ö, )• M 2 e2 = e • M2m , g2 = g • M2m , P2 = P • M2m. (11) (12) (13) (14) Here, matrix M2m represents the coordinate transformation from Sm to S2. The meshing equation for hypoid gears can be represented as: fi2 (01 ) = v12 (0, )• nUr = 0. (15) Here, the relative velocity v12 of the gear and pinion tooth surfaces E2 and E1 at M1 can be obtained from the abovementioned vectors. The parameter 01 can be obtained by solving Eq. (15). The related parameters can be obtained by inserting the value of 01 into Eqs. (7) to (11). Uldr •tldrXnld -UldrXnidr ^Contact pattern Fig. 2. The principal directions on the tangent plane Fig. 2 shows the principal directions on the tangent plane at the reference point. The so-called first principal direction denotes the direction of maximum curvature of tooth surface; the second principal direction denotes the direction of minimum curvature of the tooth surface. In general, the two directions are perpendicular to each other. The first or second principal curvature or torsion is the curvature or torsion in the corresponding first or second principal direction. The first principal curvature A2 in the gear tooth surface is represented as: 1 A =—. (16) Rn1 Here, Rn1 denotes the radius of curvature of the inside blade of the gear head-cutter. The first principal torsion B2 and the second principal curvature C2 of the gear tooth surface are equal to 0. According to the Baxter method [25], the induced normal curvatures AA, AC, and torsion AC of the two conjugate tooth surfaces (E2 and Ej) at M1 can be represented as: 0 = ev (0* ), AC = AC (0; ), AA = AC tan2 ev, AB = -AC iand,,. (17) (18) (19) (20) Here, the 0v is the direction angle of the contact line that is formed between the tooth surfaces E2 and Eb ev and AC can be obtained from the relative angular velocity ю12, the relative velocity v12, and the relative acceleration a12 of the gear and pinion tooth surfaces E2 and E1 at M1. The first principal curvature A0, the first principal torsion B0, and the second principal curvature C0 of the pinion tooth surface (drive side) at M1 for the two principal directions of the gear tooth surface (Fig. 2) can be represented as: A = a -AA B0 =-AB . (21) C = -AC Based on the generalised Euler and Bertrand formulas [13], the curvature parameters of the pinion tooth surface (drive side), along the two principal directions of the pinion tooth surface, can be represented as: A1dr = Ai cos2 A- 2B0 sin A cos A + C0 sin2 A. (22) C1dr = C0 cos2 A + 2B0 sin A cos A + ^sin2 A, (23) B1dr = B0 (cos2 A- sin2 A) + (A - C0 ) sin A cos A. (24) Here, A is the angle between the first principal direction of the gear tooth surface and the pinion tooth surface on the tangent plane (Fig. 2). u1dr can be represented as: U1dr = ( X nidr )sin Л + ti* C0s Л. (25) In the same way, all of the related parameters and vectors (including the curvature parameters A1co, B1co, CJco) of the coast side of the pinion tooth surface can be obtained. 2.2 Calculation of Pinion Curvature Parameters Based on Pinion-Generating Surfaces The configuration in Fig. 3 is the same as that in Fig. 1. The coordinate system Sm{Xm, Ym, Zm} is rigidly connected to the cutting machine (Fig. 3). The top (A-A), bottom and middle of Fig. 3 are the machine front view, top view and side view (projection of the head cutter), respectively. The cradle rotates about the G axis. The p axis is the projection of the pinion spindle in the XmOmYm plane. The points Om, O0 and O1 are the machine centre, the centre of the head cutter, and the cross point of the pinion, respectively. ep Omr Й AA Hp . epxG JMI Head-cutter for pinion A L Om m "mT h d^JOi- lit lit G —Pinion Fig. 3. Coordinate system applied for pinion generation The manufacturing coordinate system of the gears and pinions and the installation coordinate system of the hypoid gear set are represented in the coordinate system Sm. Taking backlash and other factors into account, for tangency at the reference point M2 of the gear and pinion tooth surfaces E2 and E1, the position vector r2 of the reference point, the unit normal n2 and the unit vector u2 of the coast side of the pinion tooth surface at the reference point M2 can be represented as: r2 = rico ■R M2 ] n2 = nico ■ RM2 ] U2 = «h» R[P ^2 ]. (26) (27) (28) Here, the matrix R [p, 02] is a transformation matrix that denotes the rotation angle в2 (known parameter) about the vector p. The pinion cradle spindle G does not always coincide with the Ym axis (Fig. 3), and it can be represented in the coordinate system Sm as: G = G(ouYml ). (29) The position vector of the reference point M1 and the unit vector of the blank offset direction of the generating gear and the pinion can be represented in the coordinate system Sm as: pi (Xp, Emi ) = ridr + XpP + E, e P = e P KlYml ) e mi p ' (30) (31) The relative velocity vp1 of the generating gear and pinion tooth surfaces and E1 at M1 can be obtained using the vectors G, ap1, p, r1dr . This can be represented as: 'pi = Vpi (i , Ymi, Xp , Emi, Hl, Rai). (32) Here, aGi, у mi, Xp, Emi, H RaX, are unknown parameters. The meshing equation for the generating gear and pinion tooth surfaces at M1 may be represented in the coordinate system Sm as: fpi = fpi (оиУти , EmV H,, Ral ) = v pi ■ nldr = (33) Based on Eq. (33), Ra1 can be obtained using the following equation: Ra1 = fp1 (оиУти Xp , EmV Hl) • (34) To perform tangency at M1 or M2 of the pinion tooth surface E1 and the pinion head-cutter surface Ep, a rotation angle в3 of the pinion about the p axis and a rotation angle djRal of generating gear about the G axis are necessary (Fig. 3). The position vector from the reference point to the crossing point r4, the position vector from the reference point to the machine centre a4, the unit normal n4 and the unit vector u4 of the coast side of the pinion tooth surface at the reference point M2 can be represented as: r4 = r2 •R [PA ] n4 = n2 • R [PA ] u4 = u2 • R[рД ], a4 = r4 + ^pP + Emiep . (35) (36) (37) (38) Here, matrix R [p, в3] is a transformation matrix that denotes the rotation angle в3 (unknown parameter) about the vector p. Based on Eqs. (33) and (34), 03* can be obtained by solving the meshing equation for the coast sides of the generating gear and pinion tooth surfaces at M2. The angle 03* can be obtained using the following equation: 03* = 0* {(X^Y^ Xp, E^ H, ). (39) Given a rotation angle 03*/Ra1 of the pinion-generating surface about the G-axis, the position vector a6 and the unit normal n6 can be represented as: a, (/ml, ' EmU Hl ) = = Hl G + a4 • R Ra1 G n6 (ovYmi,Xp,Eml,H1 ) = n4 •R G ^3 G' t: (40) (41) The unit vector of the pinion head-cutter axis can be obtained using the vectors a6, n6, ap1, n1, and can be represented as: c = c (иУти Xp, Emi,H,, i-m, rn 2). (42) The unit vector of the first principal directions of the pinion-generating surface at M1 can be represented as: c X n,. c X n,. (43) Fig. 4. Projection of pinion head-cutter in the coordinate system Sm O As in Eqs. (11) to (14), the parameters and vectors of the gear and pinion tooth surfaces are replaced by the parameters and vectors of the pinion-generating surface and the pinion tooth surface. Thus, the vectors a3, e3, G3and p3 are obtained in the coordinate system S1 (ij, njxij, nj} from the vectors ap1, ep, G, p. Relative angular velocity юрЬ relative velocity vp1, relative acceleration apj of generating gear and pinion can be obtained from the abovementioned vectors. The first principal curvature of the pinion-generating surface are represented as: A - -1 Apdr ~ . rnl (44) The first principal torsion Bpdr and the second principal curvature Cpdr of the pinion-generating surface are equal to 0. As in Eqs. (17) to (21), the curvature parameters of the pinion tooth surface, along the principal directions of the pinion-generating surface can be obtained. The angle Д1 that is formed between the first principal directions of the pinion-generating surface and the pinion tooth surface on the tangent plane can be represented as: sinAj = Uj • ij. (45) As in Eqs. (22) to (24), the curvature parameters A\dr, B\dr, C\dr of the pinion tooth surface, along the principal directions of the pinion tooth surface, can be obtained using the generalised Euler and Bertrand formulas. In the same way, all of the related parameters and vectors (including the curvature parameters A\co , B'1co , C\co) of the coast side of the pinion tooth surface can be obtained. 2.3 Determination of Basic Machine-Tool Settings for Generating the Pinion The theoretical outside and inside blade angles for the pinion can be represented as: (see Fig. 4) (46) (47) Theoretically, the curvature parameters of the pinion tooth surfaces along the principal directions of the pinion tooth surface should be the same as that determined by the two abovementioned methods; the theoretical outside and inside blade angles for the pinion should also be equal to the actual blade angles of the head-cutter for the generating pinion. Therefore, the seven equations can be written as follows: fb A1dr + A1co = fb A1 B1dr = B1dr C1dr = C1dr B1co = B1co C1co = C1co а 1 + A 1dr л1со Here, fb denotes the tooth-bearing length unbalancing factor. The seven unknown parameters (aG1, ym1, Xp, Em1, Hi, rn1, r„2) can be obtained by solving Eq. (48). The position vector of the pinion head-cutter generating surface at the blade tip midpoint P can be represented as: Г = hfmlC- 0.5 (r (49) Here, hfm1 denotes the pinion mean dedendum in the hypoid gear dimensions. The position vector of the pinion head-cutter centre O0 can be represented as: s = rninidr + ridr + H1C - XpP - Emiep , (50) H1 = c - rninidr - О . (51) Here, H1 denotes the projection distance of M1 on the pinion head-cutter axis. A rotation angle в4 of the pinion-generating surface about the G axis is necessary for meshing the contact of the pinion-generating surface and the tooth surface at P. The value of 04* can be obtained by solving the meshing equation. The root angle of the pinion Sr can be represented as: sin^r = P •n pm ). (52) Here, npm (04*) denotes the unit normal vector after a rotation angle of the pinion-generating surface about the G axis 04*. Generally, Sr is not equal to the root angle d^ of the hypoid gears' blank dimensions. Therefore, the resulting new mean dedendum bm1 of the pinion is different from the mean dedendum Afm1 of the hypoid gears' blank dimensions. Using the modified mean dedendum bm1, the vectors rt, s and H1 should be recalculated. The outside and inside cutter point radii for the pinion can be obtained as follows: cosa. - - Hj tana4 (53) cosa. - H2 tan ab2. (54) After a rotation angle 04* of the pinion-generating surface about the G axis, the position vector of the pinion head-cutter centre and the unit vector of the pinion head-cutter axis c1 can be redefined as: s x = s • R [g,04* ] + HfilG, C = c • R [g,04* ]. (55) (56) The rest of the basic machine settings for the pinion can be obtained as follow: Xb1 - -sx • G, (57) H = sX-(ep XG), (58) V =-s x • e p, (59) sin / = c X G . (60) 3 NUMERICAL EXAMPLE In this section, a hypoid gear design software based on the abovementioned calculation strategy for basic machine-tool settings was developed; the software development flow chart is shown in Fig. 5. Fig. 5. Spiral bevel and hypoid gear design software development flow chart The theoretical analysis is performed based on a hypoid gear set generated by using the duplex helical method. The design parameters for the face-milled hypoid gear set are listed in Table 1. The basic machine settings are listed in Table 2. Table 1. Design data Design features Pinion Gear Number of teeth 7 43 Module [mm] 6.861 Face width [mm] 43.73 40.00 Pinion offset [mm] 25.4 Shaft angle [°] 90 Mean spiral angle [°] 45 33.75 Hand of spiral LH RH Cutter radius [mm] 114.3543 114.3 Table 2. Basic machine settings for the duplex helical method Applied settings Pinion Gear Radial distance [mm] 114.2545 117.4921 Tilt angle [°] 15.7363 0.0000 Swivel angle [°] -31.6295 0.0000° Blank offset [mm] 25.0224 0.0000 Machine root angle [°] -9.0996 70.2509 Machine centre to cross point [mm] 0.3431 9.6518 Sliding base [mm] 23.8256 0.0000 Ratio of roll 5.9651 0.0000 Cradle angle [°] 66.8700 70.9771 Helical motion velocity coeff [mm/rad] 11.5478 0.0000 b) Fig. 6. Tooth bearings for the hypoid gear drive 7*43; a) drive side, and b) coast side Figs. 6 and 7 show the tooth bearings and the transmission error functions on the drive and coast side by the duplex helical method for the hypoid gear a) drive 7x43 (number of pinion teethxnumber of gear teeth), respectively. In Fig. 6, the red point is the first contact point or reference point. It is obsevered that a continuous and negative function is obtained for transmission errors, with a maximum level on the drive side for the method of approximately 15.5" (0.004306°), and on the coast side of approximately 15.9" (0.004417°). applied to the pinion are 500 Nm and 1000 r/min, respectively. a) Pillion rotation angle [°] -10 -5 0 5 10 X Pitch £ \ \ E / V j -A j \ O j 5i j l / -10 \ / \ 1 / \ ' ..-.15 ... V \ 1 / / \ / / \ -20 1 \ \ b) Fig. 7. Function of transmission errors for the hypoid gear drive 7*43; a) drive side, and b) coast side Fig. 8. Finite element model for the hypoid gear drive 7*43 The load history, especially the transfer of load between neighbouring gear pairs, is very helpful for understanding the gear mesh characteristics. Fig. 8 shows a finite element mesh; five teeth of the gear are used to save costs. Fig. 9 show the evolution of contact stresses of the pinion and gear by the method for the hypoid gear drive. The torque and the rotational speed 1000 rt cx Ш) o X 600 « fU H 400 1Л a 200 C o U 0 —Gear contact stress Pinion contact stress i i i i о 15 5 30 Time X 10 3[s] Fig. 9. Evolution of contact stresses for the hypoid gear drive 7*43 20 a) b) I Fig. 10. Tooth bearings on the hypoid gear test machine for the hypoid gear drive 7*43; a) drive side, and b) coast side The analysis of the comparison of contact stress evolution is as follows: (1) There is higher contact stress due to impact in the initial phase; with the smoothing of the rotational speed, contact stress is gradually stabilizing in Fig. 9. (2) The maximum gear tooth contact stress for the hypoid gear drive is about 300 MPa to 400 MPa; the maximum pinion tooth contact stress is about 150 MPa to 200 MPa after 4.5 ms. There is no problem of the appearance of areas of severe contact stress and edge contact throughout the process. To verify the effectiveness of the generalized theory of the duplex helical method, the cutting and the rolling test experiments are done. Fig. 10 shows the tooth bearings on the hypoid gear test machine for the hypoid gear drive 7*43 by the duplex helical method. The results of the rolling test and the results of TCA in Fig. 6 are carried out under light load. The shape and location of the tooth bearings in Fig. 10 are basically consistent with those in Fig. 6. The experiment results achieve the desired effect. Finally, all of the results proved the correctness of the generalized theory of the duplex helical method described in this paper. 4 CONCLUSIONS The duplex helical method is an advanced and primary manufacturing method of face-milled spiral bevel and hypoid gears. It is a completing process in which the concave and convex tooth surfaces are generated simultaneously under a single set of machine settings, and it is diffcult to obtain a set of optimal machine settings that can ensure both sides of TCA with good characteristics. To slove the problem, some valid conclusions are obtained through this investigation: (1) A general calculation method of the basic machine settings for all spiral bevel and hypoid cutting methods, including face-milling and face-hobbing, is proposed. (2) Three reference points (Mb M2, P) are used to calculate basic machine-tool settings for formating hypoid gears manufactured by the duplex helical method; they can accurately control the position and movement relationships between the generating gear and pinion. (3) The TCA and rolling test experiment results of the hypoid gear set manufactured by the duplex helical method show that the new methodology for calculating the basic machine settings achieves the desired effect. Both sides of TCA results have good characteristics. 5 NOMENCLATURE £ shaft angle [°] E pinion offset (mm) ym2 machine root angle of gear [°] XG machine centre to back of gear [mm] H horizontal setting of gear head-cutter [mm] V vertical setting of gear head-cutter (mm) q2 cradle angle of gear [°] Sr2 radial distance of gear [mm] a21 inside blade pressure angle (drive side or convex side) for gear [°] a outside blade pressure angle (drive side or concave side) for pinion [°] a12 inside blade pressure angle (coast side or convex side) for pinion [°] ß02 spiral angle of the gear generating surface [°] rcG1 cutter point radius for the convex side of the gear [mm] AA the first induced normal curvature of two conjugate tooth surfaces (£2 and £0 at M1 AB the first induced normal torsion of two conjugate tooth surfaces (£2 and £1) at M1 AC the second induced normal curvature of two conjugate tooth surfaces (£2 and £1) at M1 aG1 pressure angle of drive side of pinion-generating surface [°] Hp horizontal setting for pinion head-cutter [mm] Vp vertical setting for pinion head-cutter [mm] Em1 blank offset for the pinion [mm] Xp machine centre to back for pinion [mm] Xb1 sliding base for pinion [mm] yn1 machine root angle for pinion [°] Sr1 radial distance for pinion head-cutter [mm] q1 cradle angle for pinion [°] I1 tilt angle for pinion head-cutter [°] J1 swivel angle for pinion head-cutter [°] Ra1 ratio-of-roll of pinion and generating gear Hi helical motion velocity coefficient represents a displacement of the pinion blank along the axis of the cradle for a rotational angle of 1 radian of the cradle [mm/rad] rn1 curvature radius of the outside blade of the pinion head-cutter (drive side or concave side) [mm] rn2 curvature radius of the inside blade of the pinion head-cutter (coast side or convex side) [mm] Subscripts 1 pinion 2 gear G generating surface of gear p genarating surface of pinion dr drive side co coast side 6 ACKNOWLEDGEMENTS The authors express their deep gratitude to the National Natural Science Foundation of China (NSFC) (No. 51575533), National Key Basic Research Program of China (973 Program) (NKBRP) (No. 2011CB706800-G), Specialized Research Fund for the Doctoral Program of Higher Education (SRFDP) (No. 20120162110004) and National Natural Science Foundation of China (NSFC) (No. 51375159) for the financial support of the research. 7 REFERENCES [1] Stadtfeld, H. (2011). Tribology aspects in angular transmission systems part IV: spiral bevel gears. Gear Technology, no. January/February, p. 66-72. [2] Stadtfeld, H. (2011). Tribology aspects in angular transmission systems part VII: hypoid gears. Gear Technology, no. June/July, p. 66-72. [3] Fan, Q., Wilcox, L. (2007). New developments in tooth contact analysis (TCA) and loaded TCA for spiral bevel and hypoid gear drives. Gear Technology, no. May, p. 26-35. [4] Fan, Q. (2006). Computerized modeling and simulation of spiral bevel and hypoid gears manufactured by gleason face hobbing process. ASME Journal of Mechanical Design, vol. 128, no. 6, p. 1315-1327, D0l:10.1115/1.2337316. [5] Gleason Works (2004). Applied gear engineering - CAGEWin software customer/dealer tranining center. The Gleason Works, Rochester. [6] Gleason Works (1971). Calculation instructions generated hypoid gears duplex helical method. The Gleason Works, New York. [7] Fan, Q. (2011). Optimization of face cone element for spiral bevel and hypoid gears. Journal of Mechanical Design, Transactions of the ASME, vol. 133, no. 9, p. 091002-091001091002-091007, D0l:10.1115/1.4004546. [8] Zhang, Y., Litvin, F.L. (1995). Computerized design of low-noise face-milled spiral bevel gears. Mechanism Machine Theory, vol. 30, no. 8, p. 1171-1178, D0I:10.1016/0094-114X(95)00052-Z. [9] Wang, J., Kong, L., Liu, B., Hu, X., Yu, X., Kong, W. (2014). The mathematical model of spiral bevel gears - A review. Strojniški vestnik - Journal of Mechanical Engineering, vol. 60, no. 2, p. 93-105, DOI:10.5545/sv-jme.2013.1357. [10] Shtipelman, B.A. (1978). Design and Manufacture of Hypoid Gears. A Wiley-Interscience Publication, New York. [11] Fuentes, A., Gonzalez-Perez, I., Litvin, F.L., Hayasaka, K., Yukishima, K. (2007). Determination of basic machine-tool settings for generation of spiral bevel gears from blank data. Proceedings of the ASME International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, p. 57-68, DOI:10.1115/detc2007-34038. [12] Litvin, F.L., Zhang, Y., Lundy, M., Heine, C. (1988). Determination of settings of a tilted head cutter for generation of hypoid and spiral bevel gears. ASME Journal of Mechanisms, Transmissions, and Automation in Design, vol. 110, no. 4, p. 495-500, DOI:10.1115/1.3258950. [13] Litvin, F.L., Fuentes, A. (2004). Gear geometry and Applied Theory (2nd edition). Cambridge University Press, New York, DOI:10.1017/CBO9780511547126. [14] Astoul, J., Mermoz, E., Sartor, M., Linares, J.M., Bernard, A. (2014). New methodology to reduce the transmission error of the spiral bevel gears. CIRP Annals - Manufacturing Technology, vol. 63, no. 1, p. 165-168, DOI:10.1016/j. cirp.2014.03.124. [15] Cao, X.M., Fang, Z.D., Xu, H., Su, J.Z. (2008). Design of pinion machine tool-settings for spiral bevel gears by controlling contact path and transmission errors. Chinese Journal of Aeronautics, vol. 21, no. 2, p. 179-186, DOI:10.1016/S1000-9361(08)60023-0. [16] Su, J.Z., Fang, Z.D., Cai, X.W. (2013). Design and analysis of spiral bevel gears with seventh-order function of transmission error. Chinese Journal of Aeronautics, vol. 26, no. 5, p. 13101316, DOI:10.1016/j.cja.2013.07.012. [17] Sobolewski, B., Marciniec, A. (2013). Method of spiral bevel gear tooth contact analysis performed in CAD environment. Aircraft Engineering and Aerospace Technology, vol. 85, no. 6, p. 467-474, DOI:10.1108/AEAT-11-2012-0207. [18] Litvin, F.L., Fuentes, A., Hayasaka, K. (2006). Design, manufacture, stress analysis, and experimental tests of low-noise high endurance spiral bevel gears. Mechanism and Machine Theory, vol. 41, p. 83-118, DOI:10.1016/j. mechmachtheory.2005.03.001. [19] Lin, C.H., Fong, Z.H. (2015). Numerical tooth contact analysis of a bevel gear set by using measured tooth geometry data. Mechanism and Machine Theory, vol. 84, p. 1-24, DOI:10.1016/j.mechmachtheory.2014.09.010. [20] Fan, Q. (2007). Enhanced algorithms of contact simulation for hypoid gear drives produced by face-milling and face-hobbing processes. ASME Journal of Mechanical Design, vol. 129, no. 1, p. 31-37, DOI:10.1115/1.2359475. [21] Simon, V. (2007). Computer simulation of tooth contact analysis of mismatched spiral bevel gears. Mechanism and Machine Theory, vol. 42, no. 3, p. 365-381, DOI:10.1016/j. mechmachtheory.2006.02.010. [22] Tamizharasan, T., Senthil Kumar, N. (2012). Optimization of cutting insert geometry using DEFORM-3D: Numerical simulation and experimental validation. International Journal of Simulation Modelling, vol. 11, no. 2, p. 65-75, DOI:10.2507/ IJSIMM11(2)1.200. [23] Gonzalez-Perez, I., Fuentes, A., Hayasaka, K. (2011). Computerized design and tooth contact analysis of spiral bevel gears generated by the duplex helical method. ASME International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, p. 149-158, DOI:10.1115/detc2011-47108. [24] Fong, Z.H. (2000). Mathematical model of universal hypoid generator with supplemental kinematic flank correction motions. ASME Journal of Mechanical Design, vol. 122, no. 3, p. 136-142, DOI:10.1115/1.533552. [25] Baxter, M.L. (1973). Second-order surface generation. Industrial Mathematics, vol. 23, no. 2, p. 85-106. Strojniški vestnik - Journal of Mechanical Engineering 61(2015)9, 533-542 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2015.2731 Original Scientific Paper High-Frequency Calibration of Piezoelectric Displacement Sensors Using Elastic Waves Induced by Light Pressure Jernej Laloš* - Tomaž Požar - Janez Možina University of Ljubljana, Faculty of Mechanical Engineering, Slovenia In the study of ultrasound propagation in matter, displacement sensors are indispensable and of these, the most sensitive are piezoelectric sensors. In order to eliminate the intrinsic effects of the sensor from the measurements, the sensor has to be properly calibrated, which means that its transfer function has to be evaluated from a known sensor input signal and a measured sensor output signal. This has usually been done by comparing the sensor response signal to a known input signal, namely, an ultrasonic waveform, which can be theoretically calculated using mathematical models and numerical algorithms. Until now, the point-source-point-sensor model has been primarily used, while ultrasonic waves were induced mechanically either by a dropped ball or a capillary fracture. In this paper, a real-source-real-sensor model is presented. It provides a more faithful waveform construction and it enables the removal of the aperture effect from the calculated sensor transfer function, thus giving correct and universal sensor response characteristics. This was corroborated by high-frequency calibration measurements of the output signal of a Glaser-type conical sensor in two positions on both surfaces of a glass plate, while ultrasonic waves were induced by the radiation pressure of a nanosecond laser pulse. Keywords: absolute sensor calibration, piezoelectric sensor, laser pulse, aperture effect, elastic waves in a plate, optodynamics, Green's function Highlights • An absolute, high-frequency calibration method of a Glaser-type piezoelectric displacement sensor. • A way of incorporating the real dimensions of the sensor contact area and the stimulant signal origin area. • Elimination of the sensor aperture effect and the source force distribution effect from the calculated transfer function. • Detailed separation of spectral characteristics by their individual contributors. 0 INTRODUCTION Displacement measurements are indispensable in the study of ultrasound propagation in matter. Such measurements may be carried out using sensors that operate through transduction mechanisms based on piezoelectric, electrostatic, electromagnetic or optical (mainly light interferometric) principles [1]. Of these, piezoelectric sensors are the most sensitive [2], making them especially applicable for low-amplitude measurements (~ 1 pm). Piezoelectric sensors are being used in mechanical engineering, as well as for scientific and technological research. They are used for acoustic emission testing [3] and [4] and nondestructive examination of materials [4] and [5], in microseismology [6] and [7], wave propagation studies [8], light-matter interaction studies [9] and [10], optodynamics [11], and many other applications. For the measurements to be accurate, it is essential that the sensor is properly calibrated [12]. Its transfer function has to be known so that the sensor's intrinsic effects can be removed from its output signal and displacement measurements exclusively obtained. Sensor calibration should be carried out uniquely for the material on which the sensor is intended to perform future measurements. While each material has its specific acoustic impedance, an equal mechanical disturbance may produce a different sensor displacement due to their impedance mismatch. Such calibrations have usually been carried out by comparing the response of the sensor under consideration to a known stimulant signal, namely, a surface waveform, which has been either measured by another, precalibrated sensor or calculated theoretically [12] and [13]. With the use of a precalibrated sensor, both calibrations have to be carried out using the same materials, the same sensor apertures and at the same relative positions. While this is also true for the theoretical calculation, the latter is more universal and less laborious as the parameters are easily changed in the model. Until now, the theoretical waveform calculation has been mostly carried out using mathematical approximations which consider the sensor contact area and the stimulant signal origin area to be mere points and not of real size. This was done, because such approximations were valid as the frequencies of interest were sufficiently low. For the stimulant waveform origin, a dropped ball [13] and [14], a glass capillary fracture [13] and [15] or a pencil lead break [5] has usually been used [1], and [3]. In these cases, the contact times are rather long (about 1 ^s for a dropped ball and 200 ns for a capillary fracture) and their force distributions, while concentrated in a relatively small impact area (with an estimated order of magnitude of about 10 дт to 100 дт [16]), may not be uniform over time and both may vary slightly with each repetition. All of this makes such calibrations credible only for signal wavelengths larger than the sizes of the signal origin and the sensor contact areas, which, quite undesirably, caps the calibration frequencies and results in inaccurate spectral transfer functions. Thus, the sensor aperture effect and the source force distribution effect are not addressed entirely as they are mostly avoided [13]. In this paper, an absolute high-frequency calibration of the piezoelectric Glaser-type conical sensor [1] and [3] is presented. An absolute calibration is possible, because such sensors operate without resonance and have the flattest response function among piezoelectric sensors [1]. Although this technique follows the general guidelines presented in the standard [12] and in the paper by McLaskey and Glaser [13], it improves upon them with several new features: an optodynamic interaction as the source force of the stimulant ultrasonic signal with an expanded mathematical model that incorporates real dimensions of the sensor contact area and the source force impact area and can allocate the signal's spectral characteristics to each individual contributor. laser pulse f(f,t) sensor signal J (r,t) PLATE g jr,r,t) V elastic waves u (r,t) SENSOR // Fig. 1. Schematics of mathematical transition from the source force impulse through the plate transfer function, surface waveform and sensor transfer function to the sensor voltage output signal The main advantage of using a laser pulse to induce the stimulant signal is that it has a very short temporal distribution and a known spatial intensity profile - both of which are independent from each other - which produces a well-determined force impulse while being consistently repeatable. The use of a laser pulse also enables a controlled size variation of the source force impact area, which in turn allows variation in the stimulant signal frequency range. In this calibration, for example, the size of the impact area is relatively large in order to demonstrate the incorporation of the macroscopic spatial extent of the source in the expanded model. It can be, however, reduced to only a few wavelengths of laser light, which is about an order of magnitude less than from any known mechanical device, thus approximating a 5-source as closely as possible. Such concentration of laser light has to be used with caution, though, as the high fluence may exceed the laser-induced-damage threshold and the specimen, along with any surface coating on it, may become damaged. The expanded mathematical model enables the proper stimulant waveform to be constructed and the sensor aperture effect [17] to be correctly accounted for, resulting in a proper and more accurate sensor spectral transfer function. The sensor can thus be correctly calibrated, even for higher frequencies than before, and used for measurement and identification of the individual wave-arrivals of ultrasonic waves in acoustic emission and laser ultrasound. 1 METHOD In order for a sensor to be considered calibrated, the spectral characteristics of its transfer function have to be known. The essence of this method is, therefore, to theoretically calculate the displacement waveform that is detected by the sensor and compare it to the measured sensor output, as outlined in the standard [12] and in paper [13]. Due to the complicated nature of wave propagation in matter [8] and sensor transductivity mechanisms [15], it is convenient to introduce certain simplifications and idealizations to their mathematical description, such as the transfer function concept and the Green's function formalism. Thus the plate and the sensor are each considered to have their own transfer functions, which transform a certain time t dependent input signal into a certain different output signal. Therefore, the plate has a transfer function g (r, r, t ) which transforms an input source force signal f (r , t ) at a position r into an output displacement signal u (r, t) at another position r. Similarly, the sensor has a transfer function i (t) which transforms an input displacement signal u (r, t) to an output voltage signal 5 (r, t). The schematics of this linear transform chain are shown in Fig. 1. Both transfer functions are considered to be algebraically linear and time invariant. This is significant as it allows for the output signal to be expressed as a time convolution of the input signal and the appropriate transfer function: u (r, t ) = f (r, t ) * g (r, r, t ), s(r, t ) = u (r, t ) * i (t ). (1) (2) Since s (r, t) is measured and u (r, t) is theoretically calculated, it is useful to perform a Fourier transform deconvolution when searching for i (t). Fourier transform of Eq. (2) gives: 5 (r ,m) = U (r ,m) I (m). (3) To obtain the spectrum of the sensor transfer function I (m), one must simply divide the other two transforms: IИ = S (r ,m) U (r,a>y (4) This is permissible because U(r, m) is not zero at any frequency m=2 n v where it is defined. The plate transfer function g (r, r, t ) is a superposition of individually weighted Green's functions based upon the distance distribution between points in the source force area and the sensor contact area on the plate's surface. Green's functions, in general, are the solutions of wave equations for a 5-function source impulse and are highly specific for each material, shape and distance. Each Green's function here is calculated using the modified numerical algorithm developed by Hsu [18]. It calculates Green's functions in the infinite planeparallel plate approximation, which means that it accounts for the direct waves and the waves multiple reflected from the top or the bottom surface of the plate but not for those reflected from the sides of the finite plate. Four distance distributions and four plate transfer functions are used and compared in the calibration discussion further on. One is a simple point-source-point-sensor (PP) model gPP(r, r, t) where both the source force and the sensor are considered to act upon and from only one point on plate's surface - such as it has been known and used until now [13]. The real-source-real-sensor (RR) model g^Cr,r,t) takes faithfully into account real dimensions of the source force impact area and the sensor contact area. The third and the fourth one are the point-source-real-sensor (PR) model gPR(r,r,t) and the real-source-point-sensor (RP) model gN(r, r, t), which are in between the former two and are used to garner insight into the spectral characteristics and their contributors. Each of these plate transfer functions, of course, produces its own waveform model and from each of those a different sensor transfer function is calculated. Since the arrangement is such that the source force acts only normally on the plate's surface and the sensor is assumed to detect only out-of-plane displacements, only the z direction is of concern and only the gzz component of the elastodynamic Green's tensor is used (this direction notation, however, is omitted for simplicity). 2 SENSOR CALIBRATION The sensor used here is the piezoelectric Glaser-type conical sensor: 'SteveCo ' KRNBB-PC Point Contact Sensor, produced by KRN Services, Inc. [19]. This sensor is designed to detect small-scale, absolute out-of-plane displacements. The detection element in this kind of a sensor is a piezoelectric crystal PZT-5A in the shape of a truncated cone covered with a nickel electrode and backed by a heavy brass mass of irregular shape surrounded by rubber and encased in a steel case [1]. The sensor is assumed to have a uniform sensitivity over its contact area, which is circular in shape with a radius of rS0 = 0.5 mm. The spatial distribution of the waveform under the sensor's aperture it therefore averaged out and combined into one output signal by the sensor. The geometric properties and electrical wiring of the sensor are specifically designed to minimize the ultrasonic and electromagnetic distorting effects. In this paper, all electrical support equipment, such as amplifiers, cables, and oscilloscope, are considered as part of the sensor. 0.06 0.05 0.04 0.03 0.02 0.01 0.00 n(t) [1/ns] A ; I \ ^ \ f \ t [ns] 0 20 40 60 80 100 120 140 Fig. 2. Normalized temporal distribution of the incident laser pulse The source force f (r, t ) is the result of a laser pulse light pressure on the highly reflective mirror surface [9] and [10]. The laser pulse is produced by Nd:YAG laser of wavelength ÄL = 1064 nm and energy El = 200 mJ. The pulse, which was measured with a fast photodiode, has a very short temporal distribution n(t) with a full-width-at-half-maximum of AtL = 17 ns, as shown in Fig. 2. It is circularly symmetric and has a top-hat spatial profile [10] with a radius of rL0 = 1.75 mm on the plate's surface. The substrate of the mirror is a plane-parallel disk-shaped glass plate made of UV-grade fused silica (SiO2). It has a diameter of d = 50 mm and a thickness of h = 12 mm. Its elastic properties are: mass density p = 2200 kg/m3, Young's modulus Y = 72 GPa and Poisson's ratio f = 0.17. The top side of the plate is coated with a highly reflective (HR) layer with reflectivity of ^HR > 99.8 % at wavelength ÄL. This layer does not affect the reflection of the mechanical waves as its thickness is much smaller than the wavelength of even the shortest detected waves. a) LASER PULSE Fig. 3. Experimental setup schematic showing both positions of the piezoelectric sensor relative to the laser pulse as well as the most significant geometric properties in a) side view and in b) top view The calibration arrangement is shown in Fig. 3. The laser pulse is incident on the top of the plate and is mostly reflected by the HR layer. The small amount of light that does pass through the plate is absorbed only insignificantly. In this manner, the pulse delivers a force impulse of J = 2EL/c0 = 1.33 nNs normally to the plate's surface; where c0 is the speed of light in vacuum. This impulse, in turn, generates elastic waves that propagate through the plate. These waves are therefore only light-pressure-induced and not thermoelastic or ablation-induced at all [9]. The sensor is deployed at two different positions during the calibration process. In position 1, at r1, it is placed on the top surface of the plate next to the laser pulse impact area at r with their respective centers r0 = 11.5 mm apart. In position 2, at r2, the sensor is placed on the opposite side of the plate, directly beneath position 1. In such an arrangement, waves reflected from the sides of the plate may reach the sensor as well: the earliest after about 8.0 (is and the others well after 9.0 (s. It was experimentally found that the first one, the surface skimming P-wave, has such a low amplitude that it can be disregarded in this case. Therefore, for all practical purposes, the plate can be considered as infinitely large in the time period of at least the first 9.0 (is after the laser pulse illumination. Even a small amount of absorbed light by the PZT element greatly disturbs the elastic wave displacement measurement. For that reason, the sensor cannot be placed directly in the path of the laser pulse. To additionally minimize light absorption by the PZT element, a 20-(m-thick gold foil with some couplant was inserted between the sensor and the plate. This was done because gold has higher reflectivity (RAu = 99 % [20] and [21]) at laser light wavelength 4 than the nickel (RNi = 73 % [20] and [21]) that covers the sensor's tip. The calibration measurements of the sensor output signal are averaged out of N = 200 repetitions in each of the two sensor positions to reduce the stochastic noise and improve the signal-to-noise ratio. Each measurement has its start time set to zero by a photodiode with a rise time of 1 ns, which is triggered when the laser pulse illuminates the plate's surface. The sampling period is 5t=2 ns while an individual measurement lasts tn=9.0 (is, thus ending well before any significant waves reflected from the sides of the plate reach the sensor. Since the waveform amplitudes are small, it is essencial to reduce the noise in the measurements and thus improve the signal-to-noise ratio. To illustrate the importance of this, the absolute values of the fast Fourier transform (FFT) spectra of the sensor outpust signals for both sensor positions, \Sr (rb v)| and \Sr (r2, v)|, averaged out over 200 repeated measurements, and the absolute values of the FFT spectra of the noise in one measurement \E1 (v)| and the noise averaged out over 200 repeated measurements |£200 (v)| , are shown in Fig. 4. The noise measurements were carried out under the same circumstances as the signal measurements except that, obviously, no stimulant force was introduced. 101 10' 10 10-Г 10-3r 10 10' 10" ^ 10 Fig. 4. A spectral comparison between the absolute values of the FFTs of the sensor output signal in both sensor positions and the sensor output noise averaged over 200 measurements and another measured only once Comparing both noise spectra in Fig. 4, it is evident that averaging does, in fact, reduce the noise by about a decade throughout most of the frequency range. This seems consistent with signal processing theory which states that when averaging the signal over N repeated measurements, its stochastic noise sould decrease by a factor of -JN . The figure also shows that the useful sensor response signal for both sensor positions is well above the averaged noise levels for frequencies smaller than about vz = 5 MHz. At frequencies greater than vz, however, the signal decreases below the averaged sensor output noise and becomes unsuitable for further processing. Any sensor calibration measurement must therefore be regarded as inaccurate at frequencies above vz and cannot be used for credible sensor calibration at those frequencies. It has to be noted here that much of the output noise is most likely due to the electronics of both amplifications stages of the sensor as they were not optimized for such high-frequency measurements. The signal-to-noise ratio could be improved even more with larger laser-pulse energies and with more numerous measurement repetitions. The sensor output measurements and calculated input waveforms are presented in Figs. 5 and 6. Fig. 5a shows a comparison of the two theoretically calculated sensor input waveform models in picometres: a PP model uPP (r1, t) calculated from gPP(r,rpг) and an RR model uRR(r1,t) calculated from g^Cr,r1,t), with the measured sensor output voltage 5 (r1, t) in millivolts for the sensor in position 1. In addition, Fig. 5b shows a spectral comparison of the absolute values of their respective FFTs | Upp (ri, v) |, | Urr (ri, v) | and | S (ri, v) |. a) 10 10' 10 10 10 10 10' ''r 10" 'V 10 Fig. 5. a) A comparison of the two theoretically calculated sensor input waveform models and the sensor output voltage in sensor position 1 (same side as the source force) as well as b) a spectral comparison between the absolute values of their respective FFTs Similarly, Fig. 6 shows all of these for the sensor in position 2. The distinctive peaks and abrupt changes in slope in the waveforms in Figs. 5a and 6a indicate wave-arrivals. The sensor transfer functions are calculated as described in section 1. The sensor input waveform models and the sensor output voltage measurements for both sensor positions are each fast Fourier transformed to obtain their spectra. The output measurement spectrum is then divided by the corresponding input waveform spectrum, as in Eq. (4), to obtain their sensor transfer functon spectrum. Before the input models and the output measurements are fast Fourier transformed, they are all multiplied by the appropriate tapered cosine (Tukey) window function [22] to avoid any spurious frequencies in the FFTs [12] due to the possible artificial steps between the first and last points of the waveform models and the voltage measurements. a) -н 2 3 6 7 8 9 ' 101 10° 10101010- b) [pm, mV] v — U',<,<(r:,v)| ■ — u;,p(r,,v)i ПНшм — S(r,,v)l ; v [Hz]! 105 vr 106 10' Fig. 6. a) A comparison of the two theoretically calculated sensor input waveform models and the sensor output voltage in sensor position 2 (opposite side as the source force) as well as b) a spectral comparison between the absolute values of their respective FFTs 102 101 100 101010-: 10- [mV/pm]yr „ ^-J Л II VMMJI1 l \ s. гтшпл. » i* ( « 1 в • ......W-Jv -■— A