.. Since 1955 Strojniški vestnik Journal of Mechanical Engineering Strojniški vestnik - Journal of Mechanical Engineering (SV-JME) Aim and Scope The international journal publishes original and (mini)review articles covering the concepts of materials science, mechanics, kinematics, thermodynamics, energy and environment, mechatronics and robotics, fluid mechanics, tribology, cybernetics, industrial engineering and structural analysis. The journal follows new trends and progress proven practice in the mechanical engineering and also in the closely related sciences as are electrical, civil and process engineering, medicine, microbiology, ecology, agriculture, transport systems, aviation, and others, thus creating a unique forum for interdisciplinary or multidisciplinary dialogue. The international conferences selected papers are welcome for publishing as a special issue of SV-JME with invited co-editor(s). Editor in Chief Vincenc Butala University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Technical Editor Pika Škraba University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Founding Editor Bojan Kraut University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Editorial Office University of Ljubljana, Faculty of Mechanical Engineering SV-JME, Aškerčeva 6, SI-1000 Ljubljana, Slovenia Phone: 386 (0)1 4771 137 Fax: 386 (0)1 2518 567 info@sv-jme.eu, http://www.sv-jme.eu Print: Abografika, printed in 300 copies Founders and Publishers University of Ljubljana, Faculty of Mechanical Engineering, Slovenia University of Maribor, Faculty of Mechanical Engineering, Slovenia Association of Mechanical Engineers of Slovenia Chamber of Commerce and Industry of Slovenia, Metal Processing Industry Association President of Publishing Council Mitjan Kalin University of Ljubljana, Faculty of Mechanical Engineering, Slovenia International Editorial Board Kamil Arslan, Karabuk University, Turkey Hafiz Muhammad Ali, University of Engineering and Technology, Pakistan Josep M. Bergada, Politechnical University of Catalonia, Spain Anton Bergant, Litostroj Power, Slovenia Miha Boltežar, UL, Faculty of Mechanical Engineering, Slovenia Franci Čuš, UM, Faculty of Mechanical Engineering, Slovenia Janez Diaci, UL, Faculty of Mechanical Engineering, Slovenia Anselmo Eduardo Diniz, State University of Campinas, Brazil Igor Emri, UL, Faculty of Mechanical Engineering, Slovenia Imre Felde, Obuda University, Faculty of Informatics, Hungary Janez Grum, UL, Faculty of Mechanical Engineering, Slovenia Imre Horvath, Delft University of Technology, The Netherlands Aleš Hribernik, UM, Faculty of Mechanical Engineering, Slovenia Soichi Ibaraki, Kyoto University, Department of Micro Eng., Japan Julius Kaplunov, Brunel University, West London, UK Iyas Khader, Fraunhofer Institute for Mechanics of Materials, Germany Jernej Klemenc, UL, Faculty of Mechanical Engineering, Slovenia Milan Kljajin, J.J. Strossmayer University of Osijek, Croatia Peter Krajnik, Chalmers University of Technology, Sweden Janez Kušar, UL, Faculty of Mechanical Engineering, Slovenia Gorazd Lojen, UM, Faculty of Mechanical Engineering, Slovenia Thomas Lubben, University of Bremen, Germany George K. Nikas, KADMOS Engineering, UK José L. Ocana, Technical University of Madrid, Spain 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 Branko Vasič, University of Belgrade, Faculty of Mechanical Eng., Serbia Arkady Voloshin, Lehigh University, Bethlehem, USA Vice-President of Publishing Council Jože Balič University of Maribor, Faculty of Mechanical Engineering, Slovenia Cover: Experimental investigation of thrust force (Fz) and torque (Mz) values during the drilling tests of Al7075 workpiece material are presented. A full factorial set of tests was performed for all the combinations of different cutting speeds, feed rates and solid carbide tool diameters. Experiments were performed in cooperation with the University of Zaragoza and Centre of Professional Training ''Corona de Aragon". Image Courtesy: Nikolaos Efkolidis, Department of Design and Manufacturing Engineering, University of Zaragoza, Spain ISSN 0039-2480 General information Strojniški vestnik - Journal of Mechanical Engineering is published in 11 issues per year (July and August is a double issue). Institutional prices include print & online access: institutional subscription price and foreign subscription €100,00 (the price of a single issue is €10,00); general public subscription and student subscription €50,00 (the price of a single issue is €5,00). Prices are exclusive of tax. Delivery is included in the price. The recipient is responsible for paying any import duties or taxes. Legal title passes to the customer on dispatch by our distributor. Single issues from current and recent volumes are available at the current single-issue price. To order the journal, please complete the form on our website. For submissions, subscriptions and all other information please visit: http://en.sv-jme.eu/. You can advertise on the inner and outer side of the back cover of the journal. The authors of the published papers are invited to send photos or pictures with short explanation for cover content. We would like to thank the reviewers who have taken part in the peerreview process. © 2018 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. Strojniški vestnik - Journal of Mechanical Engineering 64(2018)6 Contents Contents Strojniški vestnik - Journal of Mechanical Engineering volume 64, (2018), number 6 Ljubljana, June 2018 ISSN 0039-2480 Published monthly Papers Nikolaos Efkolidis, César García Hernández, José Luis Huertas Talón, Panagiotis Kyratsis: Modelling and Prediction of Thrust Force and Torque in Drilling Operations of Al7075 Using ANN and RSM Methodologies 351 Yangzhi Chen, Zheng Li, Xiongdun Xie, Yueling Lyu: Design Methodology for Coplanar Axes Line Gear with Controllable Sliding Rate 362 Márton Gróza, Károly Váradi: Fatigue Design of Ferritic-Pearlitic Nodular Cast Iron Components with Surface Discontinuities 373 Ioana Petre, Andrea Deaconescu, Flavius Sârbu, Tudor Deaconescu: Pneumatic Muscle Actuated Wrist Rehabilitation Equipment Based on the Fin Ray Principle 383 Matej Babic, Karolj Skala, Dookhitram Kumar, Roman Sturm: New Hybrid System of Machine Learning and Statistical Pattern Recognition for a 3D Visibility Network 393 Saroj Yadav, Krishna Murari Pandey: A Parametric Thermal Analysis of Triangular Fins for Improved Heat Transfer in Forced Convection 401 Eduardo Esquivel González, José Pablo Rodríguez Arzate, Giuseppe Carbone, Marco Ceccarelli, Juan Carlos Jáuregui: A Trajectory Compensation Model for Roll Hemming Applications 412 Strojniški vestnik - Journal of Mechanical Engineering 64(2018)6, 351-361 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2017.5188 Original Scientific Paper Received for review: 2017-12-29 Received revised form: 2018-03-15 Accepted for publication: 2018-04-11 Modelling and Prediction of Thrust Force and Torque in Drilling Operations of Al7075 Using ANN and RSM Methodologies Nikolaos Efkolidis1* - César García Hernández1 - José Luis Huertas Talón1 - Panagiotis Kyratsis2 íüniversity of Zaragoza, Department of Design and Manufacturing Engineering, Spain 2Western Macedonia University of Applied Sciences, Department of Mechanical Engineering & Industrial Design, Greece Many developed approaches for the improvement of sustainability during machining operations; one of which is the optimized utilization of cutting tools. Increasing the efficient use of cutting tool results in better product quality and longer tool life. Drilling is one of the most popular manufacturing processes in the metal-cutting industry. It is usually carried out at the final steps of the production process. In this study, the effects of cutting parameters (cutting velocity, feed rate) and tool diameter on thrust force (Fz) and torque (Mz) are investigated in the drilling of an Al7075 workpiece using solid carbide tools. The full factorial experimental design is implemented in order to increase the confidence limit and reliability of the experimental data. Artificial neural networks (ANN) and response surface methodology (RSM) approaches are used to acquire mathematical models for both the thrust force (Fz) and torque (Mz) related to the drilling process. RSM- and ANN-based models are compared, and it is clearly determined that the proposed models are capable of predicting the thrust force (Fz) and torque (Mz). Nevertheless, the ANN models estimate in a more accurate way the outputs used in comparison to the RSM models. Keywords: sustainable manufacturing, AI7075, artificial neural networks, response surface methodology, thrust force, torque Highlights • Investigations of the effects of cutting parameters (cutting velocity, feed rate) and tool diameter on thrust force (Fz) and torque (Mz) in the drilling of an Al7075 workpiece using solid carbide tools are performed. • Artificial neural networks (ANN) and response surface methodology (RSM) approaches are used to acquire mathematical models. • RSM and ANN models are capable of predicting the thrust force (Fz) and torque (Mz) to a great extent. • ANN models estimate in a more accurate way the outputs used in comparison to the RSM models. 0 INTRODUCTION Drilling is an essential part of processes used in subtractive manufacturing for achieving the desired products. Its massive usage requires that any parameter optimization in the process promotes the efficiency and the greener machining. Cutting tool technology is considered one of the most significant impact factors on the sustainability of machining processes and systems. Therefore, the best choice of cutting tools and parameters for the quality of the products and the tool life criterion are both critical for the sustainability metrics of a drilling process. The correct usage of the cutting tool can affect the reduction of needed resources and energy for the production of a new one. Aluminium alloy 7075 (Al7075) is suitable for a variety of specific applications in the aerospace and chemical industries because of its excellent mechanical properties. Its strength, which is comparable to many steels and its excellent corrosion resistance are the main reasons that researchers focus on it [1] to [3]. The chemical composition of A17075 is given in Table 1. The development of an appropriate cutting factor analytical model is challenging. Nevertheless, there is much research developed about the suitability of different empirical models to predict different factors during several cutting processes. A number of methods, such as artificial neural networks (ANN) and response surface methodology (RSM), have been developed for modelling many manufacturing parameters. The successful application of ANN for prediction and modelling purposes in several science and engineering domains is substantial. ANN is utilized for modelling and prediction purposes due to its advantages in nonlinear response and when time variability occurs. It covers the difficulty in inferring input/output mapping [4] and [5] and the search algorithms for optimization, based on genetic and evolution principles [6] to [8]. RSM is a reliable statistical tool for the mathematical modelling of engineering systems and for the optimization of different manufacturing processes. Table 1. Chemical composition of Al7075 Elements_Zn_M_Cu_Cr_Fe_Sí_M_Tí_Al Percentage_6_3_2_0.3_0:®_0:®_0:4_0:3_Balance *Corr. Author's Address: University of Zaragoza, Department of Design and Manufacturing Engineering, Spain. efkolidi@unizar.es 3S1 It has been considered to be useful with minimum knowledge about the process under consideration, while simultaneously reducing the resources needed for experiments. It is an appropriate methodology when many factors and interactions affect the desired responses for a given process and an effective technique for evaluating the process parameters with the least number of experiments [9] to [12]. The present paper deals with the study of the drilling process of Al7075 when using solid carbide tools. A series of different diameters were used, together with a number of different cutting speeds and feed rates. A complete set of experimental work was implemented, and different mathematical models of the thrust force and torque developed were calculated. Both RSM and ANN methodologies were applied, and although their results were highly accurate, the ANN model was more effective. 1 METHODS As mentioned before, a number of researchers have focused on the development of empirical models to predict different factors during several cutting processes. Nouioua et al. [13] developed ANN and RSM models related to cutting force and surface roughness during the turning process of X210Cr12 steel under dry, wet, and MQL machining. Multilayer-coated tungsten carbide inserts with various nose radii were used for the experiment. Measurements were taken combining cutting speed, feed rate, and cutting depth. The comparison between ANN and RSM models showed that ANN models are more accurate for the prediction of surface roughness and cutting force. Kumar and Chauhan [14] used RSM and ANN methodologies in turning, for the validation of the results obtained during experimentation and the prediction of the behaviour of the system under any condition within the operating range. Al7075 hard ceramic composite and Al7075 hybrid composite were under investigation according to the effect of cutting parameters (cutting speed, feed rate, and approach angle) on roughness, using a polycrystalline diamond tool (PCD). Furthermore, RSM, ANN and support vector regression (SVR) in turning operations were used [15] for the development of empirical models for predicting surface roughness, tool wear, and power required. These response parameters were mainly dependent upon cutting velocity, feed, and cutting time. The result was that ANN and SVR models are much better than RSM models for predicting the three response parameters. ANN and multiple regression approaches were also used for the measurement of the surface roughness of AISI 1040 steel during turning at different cutting parameters, including speed, feed, and depth of cut [16]. The full factorial experimental design was implemented to increase the confidence limit and reliability of the experimental data. Finally, it was proved that the ANN model estimated the surface roughness more accurately than the multiple regression model did. Many researchers have inestigated the appropriateness of different empirical models to predict burr size during drilling. Karnik et al. [17] described the comparison of the burr size predictive models based on ANN and RSM on AISI 316L stainless steel workpiece with cutting speed, feed and point angle as the critical parameters. The comparison showed that the ANN models provide more accurate predictions than the RSM models do. The result was the same for Qigek et al. [18] for the investigation of the effects of cutting parameters (i.e. cutting speed, feed rate) and deep cryogenic treatment on thrust force in the drilling of AISI 316 stainless steel. Using the same techniques Mayyas et al. [19] investigated the influence of cutting speed, feed and volume fraction of the reinforcement particles used on the thrust force and torque in the drilling processes of self-lubricated hybrid composite materials. While Bajic et al. [20] examined the influence of three cutting parameters (cutting speed, feed per tooth and depth of cut) on surface roughness, tool wear and cutting force components in a face milling as part of the off-line process control. Alharthi et al. [21] developed ANN and regression analysis models for the prediction of surface roughness in a face milling of an AZ61 magnesium alloy workpiece, for different spindle speed [rpm], depth of cut [mm], and table feed [mm/ min]. The coefficient of determination was found to be sufficiently accurate for the best neural network and regression analysis model from the comparison of the models with thirteen experimental validation tests. Rooki et al. [22] described a simple and more reliable ANN method and multiple linear regressions (MLR) for the prediction of cutting concentration during foam drilling operation. The results indicated the high ability in the prediction of ANN methods. Kahraman [23] investigated the predictability of penetration rate for the diamond drilling from the operational variables and the rock properties, such as the uniaxial compressive strength, the tensile strength and the relative abrasiveness. Finally, Chavoshi [24] using methods, such as regression analysis (RA), ANN, and co-active neuro-fuzzy inference system (CANFIS) attempted to predict the surface roughness, material removal rate, and over-cut of SAE-XEV-F valve-steel during electrochemical drilling in NaCl and NaNo3 electrolytic processes. 2 EXPERIMENTAL 2.1 Experimental Settings In this study, an Al7075 plate was used as the workpiece material (150 mm x 150 mm x 15 mm). The drilling tests were performed using a HAAS VF1 CNC machining centre with continuous speed and feed control within their boundaries. The cutting forces were measured by utilizing a Kistler four components dynamometer (type 9123) with all the appropriate accessories. The dynamometer signals were processed via charge amplifiers and an A/D converter to a personal computer. The measured thrust force and torque were displayed and analysed to implement an early error detection strategy. During the drilling tests, thrust force (Fz) and torque (Mz) values were measured in 36 experiments which were performed by using solid carbide drill tools (Kennametal - multilayer TiAlN-PVD-coated universal fine-grain grade). Fig. 1 shows details of the cutting tool geometry. v\ 1 D 7 L4 max t V L3 LS L (a) Catalogue 1 [L4) max drill (L3) flute (LS) flute (D)shank number depth [mm] length [mm] length [mm] diameter [mm] B041A08000CPG 29 41 36 8 B041A10000CPG 35 47 40 10 B041A10000CPG 40 55 45 12 B041A10000CPG 43 60 45 14 (b) Fig. 1. Cutting tool; a) geometry details; and b) dimensions The cutting tools are not through coolant category tools; therefore, coolant fluid (90 % water, (a) (d) Fig. 2. Experimental workflow; a) The required hardware equipment, b) Drilling process, c) Development of mathematical models, d) Comparison between experimental and predicted values for thrust force and cutting torque 10 % KOOLRite 2270 coolant) was provided by the delivery system near the cutting tool. The feed rates of 0.2 mm/rev, 0.4 mm/rev, and 0.6 mm/rev were used together with cutting velocity values of 10 m/ min, 40 m/min, and 70 m/min. The constant depth of the holes drilled was 15 mm. The full factorial set of tests was performed for all the combinations of cutting speeds, feed rates, and tool diameters. The workflow of the research is depicted in Fig. 2 and the cutting parameters, units, and notations are listed in Table 2. Table 2. Cutting variables used in the experiments Parameters Values Cutting velocity, V [m/min] 10, 40 ,70 Feed rate, f [mm/rev] 0.2, 0.4, 0.6 Tool diameter, D [mm] 8, 10, 12, 14 Axial depth of cut, ap [mm] 15 Workpiece dimension [mm] 150 x 150 x 15 Fig. 3 expresses the evolution of thrust force and the cutting torque related to different feed rates and cutting speeds. According to the graph, it can be seen that when the tool diameter increases, both the thrust force and the cutting torque values are increased as expected. The same happens in the case of the feed rate. As feed rate values are increased, the thrust force and the cutting torque are then increased respectively. In contrast, the different values of cutting speed do not noticeably affect the experimental values. The importance of cutting tool diameter and feed rate is much greater than that of cutting speed related to thrust force and cutting torque. 2.2 RSM-Based Predictive Models The RSM is an accurate tool used to check the influence of a series of input variables on the response when studying a complex phenomenon. The models produced use the least square fitting in order to provide a reliable mathematical model. In this case, a full factorial strategy was followed, 36 drilling experiments were performed, and both the thrust force and torque were modelled using polynomial Fig. 3. Experimental values derived from Kistler 9123, a) Fz and b) Mz mathematical models. The following form was selected for this case: Y = b0 + b1 X1 + b2 X2+b3 X3 + b4X4 + bu X12+b22X22 + b33X32 + b44 X42 + b 12 Xj X2 + b j3 Xj X3 + bj4 Xj X4+ b23 X2 X3 + b24X2 X4 + b34 X3 X4, (1) where Y is the response, X stands for the coded values and b0, ..., b34 stand for the models' regression coefficients. Based on this mathematical model, the data acquired formed the following equations for the thrust forces [N] and the torque [Nm] respectively: Fz=-166 + 43.2 D + 0.405 V+ 1273/+ 2.10 DxD + 0.0163 V 2 + 86/2 -0.140 DxV + 64.2 Dx/+ 0.156 Vx/ (2) and Mz = 4.52 - 0.913 D - 0.00502 V- 7.31/+ 0.0548 DxD - 0.000046 V 2 - 1.39/2 + 0.000405 DxV+ 1.73 Dx/+ 0.0197 Vx/ (3) where D is the tool diameter [mm], / is the feed rate [mm/rev], and V is the cutting speed used [m/min] and solid carbide tools and Al7075 workpiece. The adequacy of the models is provided at a 5 % level of significance. ANOVA was used for establishing the validity of the developed models. The calculated values of the F-ratio of the developed models (Tables 3 and 4), are significantly increased compared to the tabulated value of the F-table at a 95 % confidence level (2212.84 for the Fz and 2908.77 for the Mz), while the P-values are 0.000, Table 3. ANOVA table for the Fz (thrust force) Source DF SS MS F P Regression 9 6239800 693311 2212.84 0.000 Residual Error 26 8146 313 Total 35 6247946 R-Sq(adj) = 99.8 % Predictor Coef. SE Coef. T P Constant -165.6 100.3 -1.65 0.111 D 43.17 16.74 2.58 0.016 V 0.4045 0.8729 0.46 0.647 f 1272.7 157.4 8.09 0.000 DxD 2.1042 0.7375 2.85 0.008 VxV 0.016343 0.006953 2.35 0.027 fxf 86.5 156.5 0.55 0.585 DxV -0.13972 0.05386 -2.59 0.015 Dxf 64.208 8.079 7.95 0.000 Vxf 0.1563 0.7375 0.21 0.834 which proves the highest correlation between data and model in each case. The validity of the models is also proved because the R-sq(adj) is very high in both cases (99.8 % for the Fz and 99.9 % for the Mz). In addition, the significant terms of the models, when a level of significance of 5 % is used, are those with a P-value less than 0.05. For the Fz these factors are: D (P = 0.016),/(P = 0.000), DxD (P = 0.008), VxV(P = 0.015) and Dx/(P = 0.000), while for the Mz the significant terms are: D (P = 0.000),/(P = 0.000), D2 (P = 0.000), Dx/ (P = 0.000), Vx/ (P = 0.000). Residual analysis was performed to test the models' accuracy; in both cases, the residuals follow the normal distribution. They follow straight lines (almost linear patterns) proving that the errors follow the normal distribution. All the scatter diagrams of the Fz and Mz residuals versus the fitted values depict that the residuals are evenly distributed on both sides of the centreline. The same is true for the residuals versus the order of the data (see Fig. 4). The accuracy achieved is very high when comparing the measured values and those calculated from the mathematical models (3 % and 5.6 %, respectively). The derived mathematical models can be considered to be very accurate and can be used directly for predicting both the thrust force and the cutting torque within the limits of the tool diameter, feed rate, and cutting speed used. Table 4. ANOVA table for the Mz (torque) Source DF SS MS F P Regression 9 322.187 35.799 2908.77 0.000 Residual Error 26 0.320 0.012 Total 35 322.507 R-Sq(adj) = 99.9 % Predictor Coef. SE Coef. T P Constant 4.5188 0.6287 7.19 0.000 D -0.9134 0.1049 -8.71 0.000 V -0.005017 0.005471 -0.92 0.368 f -7.3081 0.9862 -7.41 0.000 DxD 0.054812 0.004622 11.86 0.000 VxV -0.00004644 0.00004358 -1.07 0.296 fxf -1.3917 0.9806 -1.42 0.168 DxV 0.0004053 0.0003376 1.20 0.241 Dxf 1.73142 0.05064 34.19 0.000 Vxf 0.019729 0.004622 4.27 0.000 Fig. 4. Residuals analyses for the Fz and Mz a) Normal Probability Plot of the Residuals, b) Residuals Versus the Fitted Values, c) Residuals Versus the Order of the Data Fig. 5. Architecture of ANN with a) 3 inputs-6 to 12 hidden neurons -2 outputs, and b) 3 inputs-6 to 12 hidden neurons-1 output topology 2.3 ANN Modelling Artificial Neural Network ability to learn complex non-linear and multivariable relationships between process parameters makes them very useful in many applications. An ANN consists of a number of neurons, which are divided into the three basic layers: input, hidden, and output. The neurons between the layers are linked, having synaptic weights. One of the basic advantages of ANN is its ability to learn from the process. When the architecture of the network is defined, then, through a learning process, weights are calculated to present the desired output. The present research used a series of neural network pieces of software available for the development of a multilayer feedforward neural network. An elementary neuron with R inputs is shown in Figure 6a. Each input is weighted with an appropriate value (w). The sum of the weighted inputs and the bias forms the input to the transfer function f). Neurons can use any differentiable transfer function f) to generate their output. The standard network that is used for function fitting is a two-layer feedforward network, with a tan-sigmoid transfer function in the hidden layer and a linear transfer function in the output layer (see Fig. 6b). Assuming that the activation function used in the hidden and the outer layer is sigmoidal, the outputs of the hidden and outer layer were calculated using the following equation: f (net) = 1 1 + e" (4) Tan-sigmoid Tranfer Function a +1 -1 a =f(Wp+b) R = number of elements in input vector a ~ tansig(n) (a) (b) Fig. 6. Multilayer neuron network architecture, a) elementary neuron with R inputs, b) linear transfer function All of the original (36) experimental data were randomly divided into three data sets including training, validation and testing. The back-propagation training algorithms, the scaled conjugate gradient (SCG) and Levenberg-Marquardt (LM), were used for ANNs training. The training set used 70 % of the data to build the network, 15 % to measure network generalization and 15 % as a testing set of the neural network. Three-layer network architectures were used to predict the thrust force and torque as shown in Fig. 5. In the first case, a 3-(6~12)-2 ANN topology was used, which consists of three input nodes (tool Fig. 7. Neural network plot linear regressions for a) Fz and b) Mz diameter, cutting speed, feed rate), one hidden layer (6 to 12 neurons) and two outputs (thrust force and torque). In the second case, two different 3-(6~12)-1 ANN topologies were used as the output layer consists of one neuron corresponding to one output variable Fz and Mz, respectively. In consequence of trials, the best network architecture for the prediction of thrust force was the 3-10-1 topology. The comparison of predicted results from ANN model with experimental measurements shows that there is a very good correlation between them. It is obvious that a neural network is an excellent tool for predicted values of Fz according to experimentally measured ones. The correlation coefficient (R value) between the outputs and targets is a measure of model accuracy. The R value for the entire dataset (training, validation and testing) is 0.99967, and it represents high correlation. In addition, all the categorized R-values (training, validation and test) are very close to 1 (see Fig. 7). When comparing both the measured and the predicted Fz values the highest discrepancy observed is 2.18%. In the case of Mz, the best network architecture for the prediction of the experimental values was the 3-8-1 topology. The R value for all the dataset (training, validation and testing) is 0.99978 and in each one of them separately is close to 1. As a result, the network achieves high accuracy (see Fig. 7). When comparing the measured with the predicted Mz values the highest difference observed is 3.15%. From the output of both cases, it a) b) Fig. 8. Comparison of the experimental values with the predicted values for a) Fz and b) Mz is concluded that the training of the ANN with one output for each case offers greater accuracy than in the case of the network with two outputs. 3 COMPARISON BETWEEN RSM AND ANN MODELS A full factorial experimentation design is implemented to search for the effects of the cutting parameters (i.e. cutting speed, feed rate, and tool diameter) on the thrust force and torque in the case of drilling. After each hole was made, the measurements of thrust force and torque were documented. Artificial neural network and response surface methodology models were developed to predict the thrust force and torque using the experimental data. A comparison was established between the experimental values of Fz and Mz (see Fig. 8) and their values of ANN and RSM model, respectively. It is obvious that the predicted values from both models approximate the experimental values. Nevertheless, the predicted values from the ANN more accurately predict the experimental values than the regression analysis model. This can also be proved from the calculated percentage of error between experimental and predicted values of the models which are depicted in Fig. 9 for both the Fz and the Mz. In the case of a) b) Number of experiments Fig. 9 Error between experimental and predicted values for a) Fz and b) Mz the ANN models, the accuracy achieved was 2.18 % and 3.15 % for both the thrust force and torque, while in the case of the response surface methodology model the accuracy achieved was 3 % and 5.6 %, respectively. 4 CONCLUSIONS The optimization of a drilling process focused on cutting tools and parameters is a basic part of the entirely machining process, which influences sustainability. The aim of this study was the generation of mathematical models for the prediction of the thrust force (Fz) and torque (Mz) related to the cutting tools and other crucial cutting parameters as feed rate and cutting velocity during the drilling process. Two modelling techniques the RSM and ANN were used to predict the thrust force and torque in a series of drilling operations of Al7075. The developed models were considered to be very accurate for the prediction of the Fz and Mz within the range of the manufacturing parameters used. A number of different ANN architectures (3-(6~12)-2) and (3-(6~12)-1) have been tested to obtain the best neural network configuration in each case. Finally, the best architecture was the (3-10-1) for the case of the thrust force and the (3-8-1) for the case of torque. The outcomes proved that the accuracy achieved was 2.18 % and 3.15 % when using the ANN models for both the thrust force and torque measured, respectively. In the case of the RSM model, the accuracy achieved was 3 % and 5.6 %, respectively. As a result, both strategies are suitable for modelling and predicting thrust force and torque when drilling is concerned, but the application of ANN is slightly more accurate than the RSM one. 5 ACKNOWLEDGEMENTS The authors would like to thank the Direction Team of the C3 CAD/CAM/CAE) Laboratory in the University of Applied Sciences of Western Macedonia, for its valuable help, making available human and technical resources. 6 REFERENCES [1] Ucun, I. (2016). 3D finite element modelling of drilling process of AI7075-T6 alloy and experimental validation. Journal of Mechanical Science and Technology, vol. 30, no. 4, p. 18431850, D0l:10.1007/s12206-016-0341-0. [2] Suryaa, V.R., Kumara, K.M.V., Keshavamurthya, R., Ugrasenb G., Ravindrac H.V. (2017). Prediction of machining Characteristics using artificial neural network In wire EDM of AI7075 based in-situ composite. 5th International Conference of Materials Processing and Characterization (ICMPC) Proceedings, p. 203-212, D0l:10.1016/j.matpr.2017.01.014. [3] Karabulut, S., Karakoc, H. (2017) Investigation of surface roughness in the milling of AI7075 and open-cell SiC foam composite and optimization of machining parameters. Neural Computing & Applications, vol 28, p. 313-327, D0I:10.1007/ s00521-015-2058-x. [4] Osman, A., Mohd, N., Suziah, S., Wan F. (2009). Study of genetic algorithm to fully-automate the design and training of artificial neural network. International Journal of Computer Science and Network Security, vol. 9, no. 1, p. 217-226. [5] Leite, W.O., Campos Rubio, J.C., Duduch J.G., Maciel de Almeida P.E. (2015). Correcting geometric deviations of CNC Machine-Tools: An approach with Artificial Neural Networks. Applied Soft Computing, vol. 36, p. 114-124, D0I:10.1016/j. asoc.2015.07.014. [6] Patra, K., Jha, A.K., Szalay, T., Ranjan J., Monostori L. (2017). Artificial neural network based tool condition monitoring in micromechanical peck drilling using thrust force signals. Precision Engineering, vol. 48, p. 279-291, D0I:10.1016/j. precisioneng.2016.12.011. [7] Qin, X., Wang, B., Wang, B., Li, H., Jiang, Y., Zhang, X. (2014). Delamination analysis of the helical milling of carbon fiber-reinforced plastics by using the artificial neural network model. Journal of Mechanical Science and Technology, vol. 28, no. 2, p. 713-719, D0I:10.1007/s12206-013-1135-2. [8] Cornea, R., Natha, C., Mansorib, M., Kurfessa, T. (2016). Enhancing spindle power data application with neural network for real-time tool wear/breakage prediction during inconel drilling. Procedia Manufacturing, vol. 5. no. 1-14, D0I:10.1016/j.promfg.2016.08.004. [9] Qigek, A., Kivak, T., Samta§, G. (2012). Application of Taguchi method for surface roughness and roundness error in drilling of AISI 316 stainless steel. Strojniški vestnik -Journal of Mechanical Engineering, vol. 58, no. 3, p. 165-174, D0I:10.5545/sv-jme.2011.167. [10] Rao, S., Sethi, A., Das, A.K., Mandal, N., Kiran, P., Ghosh, R., Dixit, A.R., Mandal, A. (2017). Fiber laser cutting of CFRP composites and process optimization through response surface methodology. Materials and Manufacturing Processes, vol. 32, no. 14,p. 1612-1621, D0I:10.1080/1042 6914.2017.1279296. [11] Thakre, A.A., Soni, S. (2016). Modeling of burr size in drilling of aluminum silicon carbide composites using response surface methodology. Engineering Science and Technology, an International Journal, vol. 19, no. 3, p. 1199-1205, D0I:10.1016/j.jestch.2016.02.007. [12] Boyaci, A.i, Hatipoglu, T., Balci, E. (2017). Drilling process optimization by using fuzzy-based multi-response surface methodology. Advances in Production Engineering & Management, vol. 2, no. 2. p. 163-172, D0I:10.14743/ apem2017.2.248. [13] Nouioua, M., Yallese, M.A., Khettabi, R., Belhadi, S., Bouhalais, M. L., Girardin, F. (2017). Investigation of the performance of the MQL, dry, and wet turning by response surface methodology (RSM) and artificial neural network (ANN), International Journal of Advanced Manufacturing Technology, vol. 93, p. 2485-2504, DOI:10.1007/s00170-017-0589-2. [14] Kumar. R., Chauhan. S. (2015). Study on surface roughness measurement for turning of Al 7075/10/SiCp and Al 7075 hybrid composites by using response surface methodology (RSM) and artificial neural networking (ANN). Measurement, vol. 65, p. 166-180, DOI:10.1016/j.measurement.2015.01.003. [15] Gupta, A.K. (2010). Predictive modelling of turning operations using response surface methodology, artificial neural networks and support vector regression. International Journal of Product Research, vol. 48, no. 3, p. 763-778, DOI:10.1080/00207540802452132. [16] Asiltork, I., Hunkas, M. (2011). Modeling and prediction of surface roughness in turning operations using artificial neural network and multiple regression method. Expert Systems with Applications, vol. 38, no. 5, p. 5826-5832, DOI:10.1016/j. eswa.2010.11.041. [17] Karnik, S.R., Gaitonde, V.N., Davim, J.P.A. (2008). Comparative study of the ANN and RSM modeling approaches for predicting burr size in drilling. International Journal of Advanced Manufturing Technology, vol. 38, p. 868-883, DOI:10.1007/ s00170-007-1140-7. [18] Qigek, A., Kivak, T., Samta§, G., gay, Y. (2012). Modelling of thrust forces in drilling of AISI 316 stainless steel using artificial neural network and multiple regression analysis. Strojniški vestnik - Journal of Mechanical Engineering, vol. 58, no. 7-8, p. 492-498, DOI:10.5545/sv-jme.2011.297. [19] Mayyas, A., Qasaimeh, A., Alzoubi, K., Lu, S., Hayajneh, M., Hassan, A. (2012). Modeling the drilling process of aluminum composites using multiple regression analysis and artificial neural networks. Journal of Minerals and Materials Characterization and Engineering, volo 11, no. 10, p. 10391049, DOI:10.4236/jmmce.2012.1110108. [20] Bajic, D., Celent, L., Jozic, S., Bajic, D., Celent, L., Jozic, S. (2012). Modeling of the influence of cutting parameters on the surface roughness, tool wear and cutting force in face milling in off-line process control. Strojniški vestnik - Journal of Mechanical Engineering, vol. 58, no. 11, p. 673-682, DOI:10.5545/sv-jme.2012.456. [21] Alharthi, N.H., Bingol, S., Abbas, A.T., Ragab, A.E., El-Danaf, E.A., Alharbi, H.F. (2017). Optimizing cutting conditions and prediction of surface roughness in face milling of AZ61 using regression analysis and artificial neural network, Advances in Materials Science and Engineering, vol. 2017, art. ID 7560468, DOI:10.1155/2017/7560468. [22] Rooki, R., Ardejani, D.F., Moradzadeh, A. (2014). Hole cleaning prediction in foam drilling using artificial neural network and multiple linear regression. Geomaterials, vol. 4, no. 1, p. 4753, DOI:10.4236/gm.2014.41005. [23] Kahraman, S. (2016). Estimating the penetration rate in diamond drilling in laboratory works using the regression and artificial neural network analysis. Neural Processing Letters, vol. 43, no. 2, p. 523-535, DOI:10.1007/s11063-015-9424-7. [24] Chavoshi, S.Z. (2011). Analysis and predictive modeling of performance parameters in electrochemical drilling process. International Journal of Advanced Manufacturing Technologies, vol. 53, no. 9-12, p. 1081-1101, DOI:10.1007/ s00170-010-2897-7. Strojniški vestnik - Journal of Mechanical Engineering 64(2018)6, 362-372 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2017.5110 Original Scientific Paper Received for review: 2017-11-25 Received revised form: 2018-03-10 Accepted for publication: 2018-03-28 Design Methodology for Coplanar Axes Line Gear with Controllable Sliding Rate Yangzhi Chen - Zheng Li - Xiongdun Xie - Yueling Lyu* South China University of Technology, China Line gear (LG) pair is a novel gear mechanism via point contact meshing on the basis of space curve meshing theory. In this paper, aiming to control the sliding rate, related to the shape of the contact curves of LGs, the driving contact curve of a coplanar axes LG pair was extended from a circular helix to a conic helix with variable cone angle, and the driven contact curve was subsequently obtained based on the space curve meshing theory. A new type of line tooth, of which two contact curves lie on each side, was proposed for forward-reverse transmission without backlash. The tooth surface is formed by a normal tooth profile moving along the contact curve and the tooth thickness auxiliary curve. The formulae of the contact curves, tooth thickness auxiliary curves, and spatial cylindrical tooth surfaces of the new type of line tooth are derived, which theoretically provides a foundation for the standardized production of LG pairs in industry. The calculation examples of the LG pairs show that the sliding rate can be reduced by using a conic helix as the driving contact curve. The results of the kinematics experiment show that the LG pair designed by the methodology could be capable of achieving forward-reverse transmission without backlash and with a controllable sliding rate. Keywords: line gear, sliding rate, contact curve, tooth surface, non-backlash, forward-reverse transmission Highlights • Design methodology for coplanar axes line gear with controllable sliding rate is deduced. • The mechanism can transmit rotation at an arbitrary alternate angle within 0° to 180°. • A conic helix with variable cone angle was chosen as the driving contact curve. • A new type of line tooth, of which two contact curves lie on each side, was proposed for forward-reverse transmission without backlash. • An example is illustrated in detail and simulated in the Pro/E software. • An experiment was carried out to verified the LG pair designed by the methodology could be capable of achieving forward-reverse transmission without backlash and with controllable sliding rate. 0 INTRODUCTION In 1733, Camus' Theorem [1] was proposed for the conjugation of tooth profiles. In 1956, Beam [2] first proposed the so-called 'Beveloid gears', which can be adopted for power transmission between nearly-parallel shafts and allow for adjustment of the gearing backlash, following which many studies detailed their geometry, design and manufacturing [3] to [5]. In 1962, Baxter [6] analysed the basic geometric and tooth contact of hypoid gears. In the 1980s, Litvin et al. [7] systematically proposed the gear geometry theory. In recent decades, the design and manufacturing of gear pairs transmitting rotations between intersecting shafts have been a topic of research [8] to [10]. The line gear (LG) pair, also called a 'space curve meshing wheel' (SCMW) in studies before 2014, proposed by Chen [11], is a novel transmission mechanism via point contact meshing on the basis of space curve meshing theory. According to the basic theory of LG, the motion at the meshing point should satisfy the following equation [12]: v„x p =0, (1) where, v12 is the relative velocity at the meshing point between the driving contact curve and the driven contact curve, and p is the unit principal normal vector of the driving contact curve at the meshing point. Compared with traditional gear, LG has characteristics of compact space, large transmission ratio and flexible design theory [11]. There is no undercut phenomenon on the line gear, and the minimum number of the line teeth of an line gear equals 1. It could be flexible to design for the application of the LGs to transmissions with perpendicular shafts [13], intersecting shafts [14] or skew shafts [15]. The driving contact curve of the coplanar axes of the LG pair is a circular helix, satisfying the condition of contact ratio [16], and a conjugated driven contact curve can be obtained through the coordinate system transformation based on space curve meshing equation Eq. (1), which may be a circular helix, a conic helix or a planar Archimedean spiral. An LG pair consists of a driving LG and a driven LG. The driving line teeth and the driven line teeth are uniformly distributed around the basal wheels of the driving LG and the driven LG, respectively [15]. There are mainly two forms of the line-tooth of LG, one of which is cantilever structural 362 *Corr. Author's Address: School of Mechanical and Automotive Engineering, South China University of Technology, Guangzhou, China, lvyl88@scut.edu.cn form with a circular profile, and the other is that its line tooth is constructed by a convex or a concave arc profile radially attached to the wheel body. The cantilever structural tooth is the simplest structure in previously published papers. The line tooth radially attached to the wheel body is suitable for the manufacturing process of cutting materials such as CNC machining, and be of better bearing capacity [17] and [18]. The sliding rate of a gear pair, used for describing the relative movement trend of each gear at the meshing point, is a critical impacting factor on the transmission quality [19]. The relative sliding between two tooth surfaces might bring about negative influences such as surface wear and frictional loss, so the allowable value of the sliding rate depends on different working conditions such as load, speed, lubrication condition and so on [20] and [21]. If a good elasto-hydrodynamic-lubrication (EHL) oil film is formed at a certain relative sliding speed, the surface wear will be eliminated and transmission efficiency will be increased [22] and [23], at which circumstance the main failure form of gears is tooth surface fatigue pitting. Therefore, for the gear transmission, the research result on sliding rate is the basis of gear lubrication and failure analysis. The sliding rate of the LG pair depends on the relative position parameters and the meshing radius of the LG pair. Aiming to decrease the sliding rate and reduce the gear wear, the design process usually follows the position parameter selection criterion in which the sliding rate equals zero at the midpoint of the contact curves [24]. In this paper, with the aim of modifying the sliding rate of the LG pair, a conic helix with variable cone angle was chosen as the driving contact curve. Then, the equation forms of the driving contact curve and the driven contact curve in the respective follow-up coordinate systems were unified. A new type of the line-tooth of LG was subsequently designed, which could be capable of achieving forward-reverse transmission without backlash. There are two contact curves that lie on each side of one line tooth. Two normal convex arc tooth profiles are established in the normal sections of two contact curves at any meshing point to form a line tooth entity. 1 DESIGN METHODS 1.1 Design Formulae of a Pair of Coplanar Axes LGs The axes of the driving LG and the driven LG are in the same plane, covering the parallel axes and arbitrary angle intersecting axes. The coordinate systems are established, as shown in Fig. 1. The coordinate systems o0 - x0y0z0 and o1 - x1 y1z1 are the fixed coordinate system and the follow-up coordinate system of a driving LG, respectively. The coordinate systems op-xpypzp and o2 -x2y2 z2 are the fixed coordinate system and the follow-up coordinate system of a driven LG, respectively. The driving LG rotates around z0 (z1)-axis and the driven LG rotates around z2 (zp)-axis. The shaft angle between z0 (z1)-axis and z2 (zp)-axis is denoted as 60, 80 e [0, n). With the aim of avoiding overlap between two shafts of LG pair, define o0 op = a. Fig. 1. The coordinate systems of a pair of coplanar axes LGs The transformation matrixes for the coordinate systems in the design process of the coplanar axes LG are as Eqs. (2) to (4). cos0„ 0 - sin0„ a cos0„ Mpo = 0 1 0 0 sin0o 0 cos$0 a sin00 0 0 0 1 (2) M 01 = M 2 p = cos t - sin t 0 0" sin t cos t 0 0 0 0 1 0 0 0 0 1 t t 0 0 cos- - sm- i i . t t 0 0 sm- cos- i i 0 0 1 0 0 0 0 1 (3) (4) where, t is an independent parameter of the contact curves and i is the transmission ratio. The equation of a driving contact curve in o1 - x1 y1 z1 is given as Eq. (5). x1 =-(ml + n1t sin 01)cos t R = < y1 = (m1 + nxt sin d1 )sin t , (5) z1 = -n1t cosQ1 where, 6j is half cone angle, mj and nj are conical parameters. When 0j equals to 0, the driving contact curve is a circular helix. Under the following conditions for the value of the parameter 01, 6j e (0, 0O] and 6j 4 90°, the driving contact curve is a common conic helix. When 6j equals to 90°, the driving contact curve is a planar Archimedean spiral. From Eqs. (2) to (5), the equation of a driven contact curve in o2 - x2y2 z2, conjugated with the given driving contact curve above, is reflected as follows. R2 = M2 p • Mp o • M01 • R1 x2 =[(« - m1 )cos0„ + njtsin (0 -91 )] cos - y2 = [( - mj )cos0„ + njtsin (0-91 )] sin ^. (6) z2 = -njtcos(0 -dj) + (-mj)sin0o Similar to the driving contact curve, the driven contact curve may be a circular helix, a conic helix or a planar Archimedean spiral. Eqs. (5) and (6) with different parameters 80 and 6j can generally be used to describe a circular helix, a conic helix and a planar Archimedean spiral. Therefore, in geometry, an undifferentiated form of mathematical expressions of the contact curves of an LG pair exists, which can generally be used in the research of the line-tooth entity construction and machining method of the driving LG and the driven LG. Referring to R1, aiming to eliminate the constant (a - m1) sin 0O in R , the equation of the driven contact curve in a new follow-up coordinate system o2- - x2-y2- zr, which is formed by moving the coordinate system o2 - x2y2 z2 along z2-axis by the distance (a- m1) sin 0O, is expressed as Eq. (7). R2 = x2 =[(« - m )cos$0 + raisin (0O -dl cos -7 y2 = [(« - mj )cos0O + njtsin (0 sin -. (7) z2 = -njt cos (6>0 -dj ) The transformation matrix between two coordinate systems, from o2- -x2-y2- z2- to o2 - x2y2 z2, is as Eq. (8), which is formed by moving the coordinate system: M 22. = 0 0 (a - m1 )sin0o 1 (8) After unifying the coordinate systems O1 - x1 y1 z1 and o2- - x2-y2- z2- into the coordinate system om-xmymzm, which can generally be used to describe the driving contact curve and the driven contact curve, from Eqs. (5) and (7), the undifferentiated form of the contact curves of an LG pair can be written as Eq. (9). = kj ( yt = (m nm Rfl. = -n-t sindj) cos t,. -njtj sin sin tj Zj = -njtj cos dj (9) (( = !'2) where, Rjjj represents a driving contact curve and Ri represents a driven contact curve. tj are independent parameters, tj e [j t;J. tJs and tje are starting values and ending values for meshing points, respectively. 6j are half cone angles, mj and nj are conical parameters. kj represents a parameter of the hand of the spiral. When kj equals 1, the conic helix is left-handed, and when kj equals -1, it is right-handed. The equations of a driving contact curve and a driven contact curve are uniquely determined by the parameters 0j, mj, nj, kj. The parameters of a driving contact curve, a driven contact curve and the relative position parameters should meet the following mathematical conditions: ti==i t n ml + m2sec0„ = a. Ql +02 = 00 —k = k2 = l (10) A driving contact curve and a driven contact curve of an LG pair are uniquely determined by parameters required, including relative position parameters a and 0O, transmission ratio i, conical parameters m1 and n1, half cone angle 61, starting value t1s and ending value t1e for meshing point of a driving contact curve. The values of parameters t1s and t1e should satisfy the condition of contact ratio [16]. For a special condition of that 90 equals to 90°, the value of the parameter m2 can be selected on request. 1.2 Construction of Line Tooth The motion transmission for an LG pair exists via two contact curves on the driving LG and the driven LG, respectively, which are expressed as Eq. (9). The contact curves lie on the line teeth entity with a spatial cylindrical surface [18]. The driving line teeth and the driven line teeth are uniformly distributed around the cylindrical or conical wheel bodies of the driving LG and the driven LG. Aiming to realize the non-backlash forward-reverse transmission, the number of contact curves on each line tooth equals two, including the first contact curve R™ and the second contact curve Rm2. That is, the driving contact curves in one line tooth of the driving LG include the first contact curve of the driving LG Rjl and the second contact curve of the driving LG R12, and the driven contact curves in one line tooth of the driven LG include the first contact curve of the driven LG R1 and the second contact curve of the driven LG R2J2. When the driving LG rotates forward, R^ it engages with R^. And when it rotates in reverse, R'm, it engages with RJJ2. Rmi and R1 are expressed as Eq. (9). Rm is obtained by rotating Rm an angle of y1PR around the axis of the driving LG clockwise, and R2m2 is obtained by rotating R2mi an angle of y2FR around the axis of the driven LG counter-clockwise. Because the line teeth are uniformly distributed in the circumferential direction of the wheel body, Rmlne, the first contact curve that lies on the next line tooth adjacent to the second contact curve R"2, can be obtained by rotating R™1 an angle of 2n / Nj around the axis of the LG, where N1 is tooth number of the driving LG, and N2 is tooth number of the driven LG, as shown in Fig. 2. The transformation matrixes for coordinate systems involved are as Eqs. (11) and (12). M m lrl j2,jl C0SVjFR -kj smj 0 0 kj sin9jFR C0S9jFR 0 0 0 0 1 0 0 0 0 1 (11) Mm cos- , 2n N , • 2n -k sin— j N,- , • 2n k sin— j N 0 0 cos- 0 0 , 2n 'n,. 0 0 0 0 1 0 0 1. (12) Fig. 2. The schematic diagram of contact curves The values of the parameters ^1FR and y2FR determine the relative positions of the first contact curve and the second contact curve of the LGs. The second contact curves of the driving LG and the driven LG are also a couple of conjugated curves satisfying the meshing condition, which is expressed as Eq. (13). M22' • Rm = M2p • Mp0 • M01 • R2. (13) Furthermore, the values of the parameters y1FR and y2FR determine the tooth thickness and the space width, which can be defined in the transverse plane, normal plane. and axial section. According to the meshing theory of LG, while the driving LG rotates forward or reverses at the same position, the two meshing points located respectively at the first contact curve or the contact curve of a line tooth lie on the same axial section of the LG. Therefore, the line-tooth thickness and the space width of LG are defined in the axial section, i.e., the line-tooth thickness and the space width are along the direction of a straight generatrix of a suppositional conical surface, a cylindrical surface or a plane, where the contact curves lie on. The line tooth thickness and the space width are two constant values, and they are independent of the position of the meshing point. The line tooth thickness Sj of LG is defined as the distance between two intersection points, one of which is the intersection point of Rj1 and any one axial section, the other is the intersection point of R" and the axial section mentioned above. The space width ej is defined as the distance between two intersection points, one of which is the intersection point of R" and any one axial section, the other is the intersection point of Rm1ne and the axial section mentioned above, as shown in Fig.3. Fig. 3. The schematic diagram of line tooth thickness and the space width in axial section From Eqs. (8), (11) and (12), the parameters s,- and e,- are derived. = n, (In/ Nj -j ). (14) An LG pair with non-backlash forward-reverse transmission should satisfy the correct meshing condition, that is, the line-tooth thickness of the driving LG is equal to the space width of the driven LG, and the space width of the driving LG is equal to the line-tooth thickness of the driven LG. The correct meshing condition of an LG pair is as Eq. 15. (15) From Eqs. (8), (11) to (15), the second contact curves of the driving LG and the driven LG meeting the correct meshing condition were obtained nm Rj 2 = Xj = kjmj' cos ( \ n V j N V j ( y = mj sin n t+ — j N,. \ (j = 1,2), (16) zj = -njtj cosdj where m - = m- + n-t- sin 9-. A line tooth in this paper is composed of three surfaces, including S™, Sm2 and S™3 . S™ and Sm2 are two parts of spatial cylindrical surfaces located at each side of the line-tooth. The contact curves Rm and R"2 lie on Sj and Sm2, respectively. The normal tooth profile is located in the normal section of the contact curve at any meshing point. The normal tooth profile is a section of invariant arc with the radius of r, which intersects with the contact curve. While the normal tooth profile moves along the contact curve, which takes a role of a generatrix, the continuous spatial cylindrical tooth surfaces are formed, i.e. the spatial cylindrical tooth surfaces Smj1 is formed by a section of arc moving along the Rml, which takes a role of a generatrix, and the Sm2 is formed by a section of arc moving along the R" in the same principal. yffl S j 3 is the upper surface. If 9- equals to 0, E" is a j3 cylindrical surface. If 6- equals to 90°, Sj3 is a plane. Under following conditions for the value of the parameter 6j, 6j e (0, 60], S"3 is a conical surface. Because the contact curves are spatial curves, for the definition of the positions of S" and Sm2, which must satisfy the conditions of none meshing interference [18], the second constraint should be requested. Therefore, the tooth thickness auxiliary curves are defined as the trajectories of arc centres of normal tooth profiles. There are two functions of the tooth thickness auxiliary curves, one of which is to limit the tooth profile constantly perpendicular to a tangent vector of the contact curve, and the other is to keep the arc radius r. The tooth thickness auxiliary curves comprise the first tooth thickness auxiliary curve RJlc and the second tooth thickness auxiliary curve Rmlc, lying between the first contact curve and the second contact curve, and the first (or second) tooth thickness auxiliary curve is shifted MM1 (or MM2) from every point on the first (or second) contact curve, as shown in Fig. 4. The spatial cylindrical tooth surfaces S" and Sm2 are formed by the two sections of normal tooth profiles moving along the corresponding contact curves and the tooth thickness auxiliary curves respectively. Fig. 4. The schematic diagram of the construction of the line-tooth The directions of MM1 and MM2 are defined to be perpendicular to the tangent vector of the corresponding contact curves, respectively, so that the tooth thickness auxiliary curves can limit the normal tooth profile constantly perpendicular to the tangent vector of the contact curve. The coordinate system oq -xqyq zq, a fixed coordinate system, is an auxiliary coordinate system that concisely describes the coordinates of any meshing point. The zq-axis is coincident with the meshing line, the trajectory of a meshing point in the fixed coordinate system, and the Xq-axis and the shaft of LG intersect at a right angle. The transformation matrix from the coordinate system oq -xqyq zq to the coordinate system om -xmymzm is as Eq. (17). -kj cos Oj cos tj -kj sin Oj M = mq - cosOj sin tj cos tj - sinO,. -k; sinO,. costi km. cos t. J J J J J 1 - sin Oj sin tj mj sin t cos Oj 0 In the coordinate system oq (17) yq zq, MMj and MM2 are perpendicular to a tangent vector a of the corresponding contactcurves at the meshing point M, that is, MMj and MM2 both lie on the normal plane n of the corresponding contact curves at meshing point M. Fig. 5. shows the directions of MMj and MM2 , starting from the meshing point M on the first contact curve respectively. and the second contact curve, Rn (or ) ^ \ / / j^jk y i / 1 ¿fer i j y« //S In the coordinate system oq -xqyq zq, the tangent vector a at the meshing point is expressed as: a = 0 m ' V 2 '2 nj + mj V 2 . '2 nj + mj (18) The direction of the vector x = [1 0 0]T is utilized to indicate the radially extending of line tooth. The vector k in the normal plane is perpendicular to the vector x, which is utilized to indicate the circumferentially extending of line tooth and is expressed as: 0 n. k = xxa = i 2 . '2 nj + mj i 2 . '2 nj + mj (19) In the coordinate system oq -zqyZZqq, the analytical expressions of MM1 and MM2 can be derived. sin0 nj cos0 MM1(2) = r ( x sin 0 ± k cos 0) = r '2 2 cos< ±- j 7 '2 . 2 mj + nj (20) where Q is a parameter determining the positions of and Zj2. Usually Q is equal to 30°. The analytical expressions of the tooth surfaces and S jj2 are achieved. ô is a parameter to form the surfaces. Fig. 5. The schematic diagram of MM1 and MM2 n mj 2jl(j2) = r(xsin^± k cos^) + r(xsin8 ± k cos8) sin^ + sin 8 nj (cos^ + cos 8) I '2 , 2 mj + nj ± mj (cos^ + cos 8) '2 , 2 mj + nj (21) Back to the coordinate system om-xmym zm, by utilizing the transformation matrix Mmq, the analytical expressions of MM1 and MM2 are obtained. MMj" = - r sin^ k j cos 6 cos t j cosQj sin tj sin 6. r cos^ I '2 . 2 m , + n , -kj (nj sintj + m j sin6y costj) n , cos t j - m j sin 6j sin t j MM 2m =-r sin^ m- cos 6j kj cos 6 . cos tj cos 6 sin tj sin6 r cos^ n , 2 m- + n j -k- (nj sintj + m j sin6j cost j) n j cos t j - m j sin 6j sin t j m j cos6j (22) From Eqs. (9), (11) and (22), the first tooth thickness auxiliary curves and the second tooth thickness auxiliary curves were represented as: R j = R j + MM? Rj2c = Rj2 + Mjj • MM2j, (23) where RHc and R?2c represent the first tooth thickness auxiliary curve and the second tooth thickness auxiliary curve of the driving LG, respectively. R?1c and R2c represent the first tooth thickness auxiliary curve and the second tooth thickness auxiliary curve of the driven LG, respectively. The analytical expressions of the tooth surfaces E" and E" in the coordinate system om -xmymzm are also achieved. Zmfl = R; + MM? + rfl, z;2 = r;2 + M;2 fl. (MM?+^ where rfl = r sin S r cosS I n , 2 m j + n j k j cos6y cos tj cos6y- sin tj sin 6, -k-(n- sint- + m ' sin6, cost- j \ j j j j j j n j cos tj - m j sin 6. sin t. m j cos 6j rj2 = r sinS r cosS '2 2 m- + n j kj cos 6 cos tj cos 6. sin tj sin6 -k - (n. sin t. + m ' sin6, cos t. j \ j j j j j j n j cos t j - m j sin 6j sin tj m '' cos 6j (25) The 'sweep' function in 3D modelling software is usually used to model the line-tooth entity. Firstly, the first contact curve and the first tooth thickness auxiliary curve are plotted by submitting the required parameters into the Eqs. (9) and (16), and a convex arc tooth profile is plotted with a radius of r in a normal plane of the first contact curve. The centre of the arc is located at the first tooth thickness auxiliary curve. Then, the arc tooth profile sweeps along the first contact curve and the first tooth thickness auxiliary curve, taking a role of the trajectory line, by the setting of the tooth profile to be constantly perpendicular to the first contact curve. With the use of 'Swept Blend' to the modelling of line tooth entity in Pro/E, the contact curve should be selected as 'Origin trajectory', along which the tooth profile sweeps, while the corresponding tooth thickness auxiliary curve should be selected as the 'X vector'. Then, 'Perpendicular to the trajectory' should be selected. In another way, with the use of 'Sweep' in Solidworks, the contact curves should be selected as the 'Path', and the corresponding tooth thickness auxiliary curve should be selected as the '1st guide curve'. By selecting of 'Follow path and 1st guide curve', the twist of the intermediate sections is determined by the vector from the path to the first guide curve, so that the tooth profile remains r perpendicular to the contact curve. The process above forms a tooth surface containing the first contact curve. Finally, another tooth surface containing the second contact curve is formed by the same method, and the modelling of a line tooth entity is completed. 2 COMPARISON OF SLIDING RATES The relative sliding between a traditional gear pair might cause surface wear and fictional loss, and then the transmission efficiency and fatigue life of the gears are affected. In general, the sliding rates are one of the tribological parameters and subsequently defined and calculated in theory. Similar to the traditional gear pair, the sliding rates of the LG pair are defined as the limit value of the ratio of the lengths difference between two relative arcs divided by the length of the given arc [24], they are expressed as Eq. (26). a2 = 1 - V(x M )2- +(y M )2 +( M )2 M) )2 +(y M) ) +(zM) )2 ' !(* M) )2 +(y M )2 +( M) )2 M )2 +( y M )2 +( 1?) (26) where a1 and a2 represent the sliding rates of the driving LG and the driven LG, respectively, , dx, , dy, , dz _ x i(i) =__ y kj ) _ j j z i(_) =__ j = 12 m ,, , y m ^ , m j, , j i,2" at at at Taking three couples of LG with vertical cross axes as design examples, the values of various parameters of LG pairs are shown in Table 1. The driving contact curve of the LG pair No. 1 is a circular helix, and the driving contact curves of the LG pair No. 2 and the LG pair No. 3 are conic helixes with a half cone angle of 7°. The heights of the driving LGs are all 16 mm for three examples. The driving contact curves and the driven contact curves can be derived from Eq. (9). Table 1. The various parameters of designed the LG pair No. ml [mm] «j [mm] 0j [°] m2 [mm] i 00 [°] 1 2.5 4 0 12 8 90 2 1.5 4 7 12 8 90 3 1.8 4 7 12 8 90 a) b) Fig. 6. 3D solid model; a) the LG pair No. 1, b) the LG pair No. 2 From Eqs. (9) and (26), the sliding rates of contact curves of three LG pairs are obtained and plotted with Wolfram Mathematica, as shown in Fig. 7. oi o.O- a) ct, o.o -0.1 b) Fig. 7. Sliding rates of the LG pair No. 1-3; a) Oj vs t, b) c2 vs t fl-Cj of the LG pair No. 1; 2-oJ of the LG pair No. 2; 3-oJ of the LG pair No. 3; 4-o2 of the LG pair No. 1; 5-o2 of the LG pair No. 2; 6-o2 of the LG pair No. 3) The 3D solid models of the LG pairs No. 1 and No. 2 are shown in Fig. 6, with the different shapes of the driving LG. From Fig. 7, the sliding rates of the LG pairs, of which the relative position parameters are limited, can be controlled by changing the shape of the driving contact curve. When the driving contact curve is designed by a certain conic helix, the absolute value of sliding rates can be significantly decreased. 3 A KINEMATICS EXPERIMENT A 3D printing technology, stereolithography appearance (SLA), was adopted to manufacture the samples of the designed LG pair. The surface profile accuracy of the process is 0.025 mm. Fig. 8 shows the prototypes of the manufactured LG pair No. 2; the driving LG is on the left, and the driven LG is on the right. Fig. 8. The prototypes of the LG pair No. 2 To measure the transmission ratio of the designed LG pair, a kinematics experiment was conducted by use of a test rig with a 4-DOF positioning table, as shown in Fig. 9. Fig. 9. The test rig (1-the driving LG;2-the driven LG; 3-a DC motor with the encoder;4-a encoder) The driving LG and the driven LG were fixed with encoders and the Decklink SDI to collect their angular rotations in the same time interval. The sampling interval in this experiment was set as 150 ms. Then the instant angular velocity was calculated. The ratio of the instant angular velocity of the driven LG to the driving LG is defined as the instant transmission ratio of the prototype of the designed LGs, as shown in Fig. 10. .2 8.04 » CO o 8.03 ro 8.02 J 8.01 M— S 8.00 c 1 7.98 H M I7"97 2 7.96 30 a) 60 90 Sampling point 120 150 I7-97 2 7.96 150 60 90 b) Sampling point Fig. 10. The instant transmission ratio; a) under the forward rotation and b) under the reverse rotation As shown in Fig. 10, the instant transmission ratios are fluctuating around the theoretical design value of 8, which means the results of the kinematics experiment are consistent with the desired value. The relative error is less than 0.01 %. The maximum instant relative errors under the forward rotation and under the reverse rotation are 0.314 % and 0.389 %, respectively, which may come from the following aspects: manufacturing error of the contact curves, installing error and positioning error, and data errors in collecting and calculating covering original data error, truncation error and round-off error. 4 CONCLUSION This paper proposed the design methodology of coplanar axes line gear with the controllable sliding rate for non-backlash forward-reverse transmission. The main conclusions include: (1) The driving contact curve of an LG pair with coplanar axes was extended from a circular helix to a conic helix, and the driven contact curve was derived on the basis of a given driving contact curve. This process provides a new method to control the sliding rates of the contact curves. The design examples show that the LG pair, of which the driving contact curves are conic helixes, possess better sliding rates. (2) There are two contact curves lying on two sides of a line tooth, and the two couples of contact curves that lie along the driving line tooth and the driven line tooth are two couples of conjugate curves. A line tooth entity, attachment to the wheel body, is formed by a normal arc profile moving along the contact curves and the tooth thickness auxiliary curves. The analytical expressions of spatial cylindrical tooth surfaces are also achieved. The proposed LG pair could be capable of achieving forward-reverse transmission continuously and smoothly without backlash. However, some issues for the proposed LG pair remain to be studied, covering that selection criterion of the geometric parameters, an experimental measure method of the sliding rates, the failure criteria, bending and contact strength, the efficiency of transmission, etc. N, N2 tooth numbers of driving LG and driven LG, m1, n1, m2, n2 conical parameters of driving contact curves and driven contact curves, [mm] 61, 02 half cone angles of driving contact curves and driven contact curves, [°] s1, s2 line tooth thicknesses of driving LG and driven LG, [mm] e1, e2 space widths of driving LG and driven LG, [mm] 91Fr, y2FR phase differences between first contact curve and second contact curve, a unit tangent vector of the driving contact _ curve at the meshing point, MM1 , MM2 Contact vector from M to M1 and M2 , r arc radius of tooth profile, [mm] R™, R12, Rm1, R2 first contact curve and second contact curve of driving LG and driven LG, nm nm Rllc , R12 Âm first and second tooth thickness auxiliary curves of driving LG and driven LG, s™, sm , S™ , sm2 the analytical expressions of tooth surfaces of driving LG and driven LG, oi, ct2 sliding rates of driving contact curve and driven contact curve. 5 ACKNOWLEDGEMENT 7 REFERENCES The authors gratefully acknowledge the support from National Natural Science Foundation of China (No. 51575191), China Postdoctoral Science Foundation (No. 2017M612652) and Central university basic scientific research service fee (No. 2017B0071). 6 NOMENCLATURE v12 relative velocity at the meshing point between the driving line tooth and the driven line tooth, [m/s] p unit principal normal vector of the driving contact curve at the meshing point, o0-x0y0 z0, o1 -x1 y1 z1 fixed coordinate system and follow-up coordinate system of driving LG, op-xpypzp, o2 -x2y2z2 fixed coordinate system and follow-up coordinate system of driven LG, 80 shaft angle between driving LG and driven LG, [°] a distances from point o0 to op, [mm] t1, t2 parameter variable of driving contact curves and driven contact curves, i transmission ratio, [1] Litvin, F.L., Fuentes, A. (2004). Gear Geometry and Applied Theory (2nd ed.), Cambridge University Press, New York, DOI:10.1017/CBO9780511547126. [2] Osplack, J.J., Beam, A.S. (1956). Conical Involute Gearing. US, US2747424. [3] Brauer, J. (2002). Analytical geometry of straight conical involute gears. Mechanism and Machine Theory, vol. 37, no. 1, p. 127-141, DOI:10.1016/S0094-114X(01)00062-3. [4] Liu, C.C. Tsay, C.B. (2002). Contact characteristics of beveloid gears. Mechanism & Machine Theory, vol. 37, no. 4, p. 333350, DOI:10.1016/S0094-114X(01)00082-9. [5] Staniek, R. (2011). Shaping of face toothing in flat spiroid gears. Strojniški vestnik - Journal of Mechanical Engineering, vol. 57, no. 1, p. 47-54, D0I:10.5545/sv-jme.2010.093. [6] Sugimoto, M., Maruyama, N., Nakayama, A., Hitomi, N. (1991). Effect of tooth contact and gear dimensions on transmission errors of loaded hypoid gears. Journal of Mechanical Design, vol. 113, no. 2, p. 182-187, D0I:10.1115/1.2912767. [7] Litvin, F.L., Vecchiato, D., Gurovich, E., Fuentes, A., Gonzalez-Perez, I., Hayasaka, K., Yukishima, K. (2005). Computerized developments in design, generation, simulation of meshing, and stress analysis of gear drives. Meccanica, vol. 40, no. 3, p. 291-323, D0I:10.1007/s11012-005-4020-y. [8] Tsai, Y.C., Hsu, W.Y. (2008). The study on the design of spiral bevel gear sets with circular-arc contact paths and tooth profiles. Mechanism and Machine Theory, vol. 43, no. 9, p. 1158-1174, D0I:10.1016/j.mechmachtheory.2007.08.004. [9] Alves, J.T., Guingand, M., de Vaujany, J.P. (2013). Designing and manufacturing spiral bevel gears using 5-axis computer numerical control (CNC) milling machines. Journal of Mechanical Design, vol. 135, no. 2, p. 024502, D0I:10.1115/1.4023153. [10] Simon, V.V. (2015). Design and manufacture of spiral bevel gears with reduced transmission errors. Journal of Mechanical Design, vol. 131, no. 4, p. 041007-041017, D0I:10.1115/1.3087540. [11] Chen, Y.Z. (2014). Line Gear. Science Press, Beijing. (in Chinese) [12] Chen, Y.Z., Xing, G.Q., Peng, X.F. (2007). The space curve mesh equation and its kinematics experiment. 12th IFToMM World Congress, Besançon. [13] Chen, Y.Z., Xiang, X.Y., Luo, L. (2009). A corrected equation of space curve meshing. Mechanism and Machine Theory, vol. 44, no. 7, p. 1348-1359, D0I:10.1016/j. mechmachtheory.2008.11.001. [14] Ding, J., Chen, Y.Z., Lv, Y.L. (2012). Design of Space-Curve Meshing-Wheels with Unequal Tine Radii. Strojniški vestnik - Journal of Mechanical Engineering, vol. 58, no. 11, p. 633641, D0I:10.5545/sv-jme.2012.493. [15] Chen, Y.Z., Lv, Y.L., Ding, J., Chen, Z. (2013). Fundamental design equations for space curve meshing skew gear mechanism. Mechanism and Machine Theory, vol. 70, p. 175188, D0I:10.1016/j.mechmachtheory.2013.07.004. [16] Chen, Y.Z., Luo, L., Hu, Q. (2009). The contact ratio of a space-curve meshing-wheel transmission mechanism. Journal of Mechanical Design, vol. 131, no. 7, p. 074501, D0I:10.1115/1.3116343. [17] Lv, Y.L., Chen, Y.Z., Cui, X.Y. (2015). A contact ratio and interference-proof conditions for a skew line gear mechanism. Transactions of the Canadian Society for Mechanical Engineering, vol. 39, no. 3, p. 647-656, D0l:10.1139/ tcsme-2015-0051. [18] Chen, Y. Z., Yao, L. (2016). Design formulae for a concave convex arc line gear mechanism. Mechanical Sciences, vol. 7 , no. 2, p. 209-218, D0I:10.5194/ms-7-209-2016. [19] Wu, X.T. (2009). Principle of Gearing. Jiaotong University Press, Xi'an. (in Chinese) [20] Li, T., Pan, C., Gao, F., Wang, X. (2012). Research on sliding ratios of conjugate surfaces of two degrees of freedom meshing transmission of spherical gear pair. Journal of Mechanical Design, vol. 134 , no. 9, p. 091002, D0I:10.1115/1.4006322. [21] Narayanan, G. (2016). Effect of sliding friction on spline surface failure under misaligned condition in aero engines. International Journal of Structural Integrity, vol. 7, no. 5, D0I:10.1108/IJSI-07-2015-0024. [22] Qian, S., Guo, D., Liu, S., Lu, X. (2011). Experimental investigation of lubrication failure of polyalphaolefin oil film at high slide/roll ratios. Tribology Letters, vol. 44, no. 2, p. 107115, D0I:10.1007/s11249-011-9820-8. [23] Ciulli, E., Stadler, K., Draexl, T. (2009). The influence of the slide-to-roll ratio on the friction coefficient and film thickness of EHD point contacts under steady state and transient conditions. Tribology International, vol. 42, no. 4, p. 526-534, D0I:10.1016/j.triboint.2008.04.005. [24] Ding, J., Chen, Y. Z., Lv, Y. L., Song, C. (2014). Positionparameter selection criterion for a helix-curve meshing-wheel mechanism based on sliding rates. Strojniški vestnik - Journal of Mechanical Engineering, vol. 60, no. 9, p. 561-570, D0I:10.5545/sv-jme.2013.1574. Strojniški vestnik - Journal of Mechanical Engineering 64(2018)6, 373-382 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2017.5120 Original Scientific Paper Received for review: 2017-11-30 Received revised form: 2018-03-18 Accepted for publication: 2018-04-16 Fatigue Design of Ferritic-Pearlitic Nodular Cast Iron Components with Surface Discontinuities Marton Groza* - Karoly Varadi Budapest University of Technology and Economics, Department of Machine and Product Design, Hungary Surface and subsurface discontinuities are one of the most important factors affecting the fatigue life of structural cast components. Their location, shape and size vary from component to component, most of them are completely harmless, but one critical defect can lead to inservice failure. Allowable surface discontinuity size finite element (FE) results can be a practical engineering tool for casting design, process planning, and the quality inspection process. Different methods based on the continuum and fracture mechanics applicable on the multiaxial high-cycle fatigue (HCF) and fatigue limit prediction for components with surface discontinuities are compared with experimental results on ISO1083/JS/500-7 nodular cast iron (NCI). Results also confirm, that the fatigue properties of the analysed material in standards truly represent low-end material strength. A design methodology is presented based on the Defect Stress Gradient approach for the display of an allowable surface discontinuity size FE-result for complex components under proportional loading conditions in HCF. Keywords: high-cycle fatigue, finite element analysis, fracture mechanics, surface defects, multiaxial fatigue Highlights • Fatigue test results for the ISO1083/JS/500-7 NCI material grade with and without surface discontinuities are compared and analysed. • Different methods are applied to describe the reduction of fatigue strength with the equivalent surface discontinuity size. • Allowable surface discontinuity size FE-results are plotted based on the DSG approach with the consideration of the mean-stress dependent scatter in the fatigue strength. • The proposed method to obtain allowable surface discontinuity size results is directly applicable in industrial applications. 0 INTRODUCTION The numerous advantageous characteristics of nodular cast iron (NCI) components, such as their steel-like mechanical properties and economical manufacturability, are linked with a higher level of scatter in material behaviour compared to steels. Engineering materials are classified based on their minimum mo no tonic tensile properties; therefore, a given standardized NCI grade can contain materials with 100 MPa difference in tensile and fatigue strength. Recommended values for the fatigue strength in standards and design guidelines do not always correspond to the lower-end strength material with the prescribed minimum monotonic properties. Apart from the large-scale scatter originating from differences in the microstructure and wall-thickness, internal and surface discontinuities also heavily influence the fatigue life of structural cast components. Surface and subsurface discontinuities have the most detrimental effect on the fatigue strength from these factors; their position, size, and morphology vary from component to component. A critical defect can lead to in-service failure, but most of the discontinuities do not affect the expected service life. A fatigue assessment-based quality inspection methodology is, therefore, necessary to balance in-service safety and the economical production of cast components. Allowable surface discontinuity size finite element (FE) result fields are a practical engineering tool during casting design, process planning, and the quality inspection process. The current study employs different methods to obtain the allowable surface discontinuity size as an FE-result field. Improving the accuracy of the fatigue assessment calculations through the consideration of the non-idealistic properties of components is currently a heavily researched area. Nadot et al. [1] investigated the effect of casting defects of the fatigue limit of NCI material. From the standpoint of the defect-induced fatigue process, each specimen is considered individual during testing. For the assessment of individual fracture events, the "step-by-step" method proposed by Bellows et al. [2] has proven to be essential. Cheng et al. [3] the high frequency cutoff of machined surface topography is defined, which is related to a material parameter a0. Besides, the proposed analytical expressions are applied to predict the stress concentration factors (SCFs derived analytical formulas for the modelling of the surface topography in the fatigue assessment process based on the theory of critical distance. Rotella et al. [4] plotted a defect size map based on FE-calculations for cast A357-T6 under multiaxial loading using the defect stress gradient (DSG) criterion. *Corr. Author's Address: BME, Department of Machine and Product Design, Mflegyetem rkp. 3, Budapest 1111, Hungary groza.marton@gt3.bme.hu 373 The comparison of experimental and simulation results has shown that the DSG criterion provides a good estimation of the crack initiation sites. Vincent et al. [5] have shown that the defect size that impacts the fatigue limit is largely dependent on the grain size of the material. Schönbauer et al. [6] investigated small artificial surface defects on the torsional and tensile fatigue behaviour of precipitation-hardened chromium-nickel-copper stainless steel. They concluded, that fatigue behaviour is governed by the maximum principal stress even under torsion and emphasized the importance of the crack initiation phase for defects with p > 50 ^m notch root radius. Mukherjee et al. [7] analysed the relationship between fatigue and microstructural properties in NCI material with cooling rate calculations, microstructural characterization by 2D SEM and 3D X-ray tomography and fatigue testing on sand and metal moulded samples. Casting defects can be considered to be notches with complex morphology; most of the advancements in the area of notched fatigue research are, therefore, applicable with some modifications. The effect of notches on the crack initiation life can be approached from the standpoint of fracture mechanics with the analysis of short-crack propagation. Ostash et al. [8] applied a fracture mechanics approach for macrocrack initiation life computation on notched mild-steel and aluminium specimens with promising results. A conventional stress-strain based approach was applied in the work of Gates and Fatemi [9] for multiaxial notched fatigue analysis with the utilization of the theory of critical distances and the Fatemi-Socie critical plane damage parameter. New approaches are currently being developed for a unified description of the total fatigue life; Liu et al. [10] proposed a method based on the theory of damage mechanics to describe damage evolution in a notched specimen. The current study compares different approaches to simulate the Kitagawa diagram (fatigue strength vs surface discontinuity size) for ISO1083/JS/500-7 NCI material. These methods are used to plot equivalent allowable surface discontinuity size FE-results. Experimental and literature fatigue data are compared to quantify the scatter in the defect-free fatigue properties. Design values of the surface discontinuity size-dependent fatigue strength are derived to ensure safety during the quality inspection process. 1 EXPERIMENTAL RESULTS FOR NCI MATERIAL As part of a large-scale research program with the aim of improving the fatigue design of NCI components, load-controlled fatigue tests have been conducted on specimens with different microstructures and surface properties. The material chosen for this analysis is a ferritic-pearlitic NCI, designated as ISO1083/JS/500-7 according to the international standard ISO 1083 [11]. 1.1 Results with Artificial Surface Defects and Casting Skin One hundred (100) fatigue test specimens with a rectangular cross-section were created with a complex artificial surface defect introduced by means of a special casting core. Fig. 1 shows the specimen geometry and computed tomography (CT) image from the longitudinal cross-section. The equivalent size of the artificial complex defect is 2110 ^m according to the Varea parameter from Murakami [12]. A detailed view of the surface defect in the specimens can be seen in Fig. 2. The chemical composition shown in Table 1 is consistent with the expectations to achieve a casting material for the IS01083/JS/500-7 grade. In Table 1, Sc refers to the degree of saturation, and has been calculated as Sc=%C / (4.23 - 0.3 • (%Si+%P)). An evaluation of the obtained microstructural features in accordance to ISO 945-1 [13] is summarized in Table 2. The metallographic inspections were conducted on a cross-section of the clamping area from a tensile specimen. B8*40 (DIN 50125 [14]) tensile specimens have been machined to standard ISO 1083 Y1-shaped test samples with 265 mm length. The tensile tests were carried out at room temperature according to EN ISO 6892-1 [15]. The average values of 12 tensile tests are a proof strength of Rp0.2=348 ± 21 MPa, an ultimate tensile strength of and an elongation at fracture of Rm=583 ± 43 MPa and an elongation at fracture of A = 14 % ± 2 %. Table 1. Average chemical composition of the test specimens (average of five measurements) %C %Si %Mn %P %S %Cr %Cu %Mg Sc 3.62 2.33 0.19 0.013 0.005 - 0.36 0.04 1.03 Table 2. Average metallographic properties of the test specimens (average of five measurements) Graphite Nodules Nodularity Form III Form IV Form V Form VI Pearlite cont. [%] [1/mm2] [%] [%] [%] [%] [%] [%] 10.7 ± 0.14 657±207 88.8 ± 1 0.84 ± 0.15 0.24 ± 0.1 28 ± 2.6 70.92 ± 2.6 46.6 ± 3.8 6 Fig. 1. Fatigue test specimen with casting skin and complex artificial surface defect, CT image of the cross-section of the gauge length 100,000 No. of load cycles Fig. 3. Design and 50 % failure probability S-N curves for defective ISO1083/JS/500-7 with the test results for R01 pulsating tensile load Fig. 2. Geometry of the testing section of the fatigue test specimen with casting skin and complex artificial surface defect 100,000 No. of load cycles Fig. 4. Design and 50 % failure probability S-N curves for defective IS01083/JS/500-7 with the test results for R0 5 pulsating tensile load Fatigue tests were carried out on a RUMUL MAGNODYN fatigue machine at room temperature at R01 and R05 pulsating tension loads with fixed amplitude with a testing frequency of 40 Hz. Cracks initiated exclusively from the large artificial defect. Fig. 3 and 4 show the fitted and the design (97.5 % failure probability and 95 % significance level considering the lognormal distribution of fatigue life) S-N curves with the experimental data for R01 and R05 pulsating tension load. 1.2 Results for Low and High-Strength Material Fatigue tests were conducted on low- and high-strength ISO1083/JS/500-7 material, to estimate the possible scatter in the high-cycle fatigue properties of smooth specimens. The specimens have been machined from NCI components from the production line, the distinction "low" and "high" strength refers to different batches of the same component. The tensile properties of the studied materials were determined on B8x40 (DIN 50125 [14]) tensile specimens according to the standard EN ISO 6892-1 [15] and are summarized in Table 3. The geometries of the fatigue test specimens are displayed in Fig. 5. The fatigue tests were carried out on an INSTRON 8874 servo-hydraulic fatigue machine with a fixed-load ratio of R0.05 and a fixed-load maximum with a testing frequency of 30 Hz. In the case of the low-strength material, fatigue testing was carried out at three different stress levels to obtain the full S-N curve. In the case of the high-strength cast iron, the step-by-step method [2] was applied with 20 MPa steps in the stress maximum to determine the fatigue strength at 106 cycles. Fracture surfaces were analysed, and the initiation sites were located in the near-surface region without the influence of notable subsurface or internal defects. Results are summarized in Table 4. It is interesting to note that the previous loading steps below the fatigue limit do not have a noticeable impact in the fatigue strength; therefore, the application of the step-by-step approach in the fatigue testing process is justified. Fig. 6 shows the obtained test results with the fitted S-N curve and the S-N curve based on the FKM guideline [16]. Table 3. Summary of tensile monotonic and fatigue strength at cycles for the ISQ1083/JS/5Q0-7 cast iron material at different R-ratios Designation Rp0.2 [MPa] Rm [MPa] A [%] or.1 [MPa] &R0 [MPa] 1 ^R0.05 [MPa] OR0.1 [MPa] tfR0.36 [MPa ] ^0.5 [MPa] FKM [16] 320 500 7 170 135 - - - - low-strength 328 491 15 (194) 135 - - - high-strength 378 649 9.7 (277) 188 - - - Yamabe and Kobayashi [17] 363 555 8 275 - - - - Rabb [18] 338 625 - 260 - 160 120 - artificial defect 348 583 14 (119) - 88 78 Table 4. Fatigue test results for the high-strength ISQ1083/ JS/500-7 material No. No. of cycles [MPa] <05 [MPa] 1,000,000 325 1,000,000 345 1 1,000,000 365 184 1,000,000 385 153,644 405 1,000,000 385 2 1,000,000 405 193 80,476 425 3 855,122 405 190 4 1,000,000 385 185 Fig. 5. Specimen geometry for a) high-strength material and b) low-strength material testing 350 300 250 200 150 100 50 0 10,000 x low-strength test data ----low-strength curve o high-strength test data ...................FKM X X ~xrx---x. "yxx. ® XT^'*--*___x .......x~. ~ ------- 1,000,000 100,000 No. of load cycles Fig. 6. Comparison of experimental results and the FKM calculation guideline [16] at R0Q5 pulsating tensile load for low- and high-strength ISQ1083/JS/500-7 2 FATIGUE LIFE OF IS01083/JS/500-7 NCI Table 3 provides a quick view of monotonic tensile strength and tensile fatigue strength at 106 cycles of the investigated material. If no tests were conducted at R-1 tension, an estimated value of the fully-reversed fatigue strength obtained by the Goodman approach is shown in brackets for an easier comparison of the different results. Fig. 7 shows the effect of mean stress on the stress amplitude corresponding to fracture at 106 cycles. The modified Goodman approach is used due to its simplicity and slight conservatism for the simulation of the mean stress effect. For the description of the mean test results the mean value of ctR*- and the mean value of the yield strength was used to fit the modified Goodman curve. 100 200 300 400 Mean stress [MPa] Fig. 7. Effect of mean stress on the fatigue limit of ISQ1083/JS/500-7, experimental results and the modified Goodman curve at 106 cycles The fatigue data included in the FKM guideline [16] proves to be a good estimate of the low-strength ISO1083/JS/500-7 fatigue properties. Based on FKM [16] data and the modified Goodman approach, a design curve is proposed for the defect-free material by using the values for CTm and Rp02 from the guideline. 3 FATIGUE LIFE OF MATERIAL WITH SURFACE DISCONTINUITIES The current study employs four concurrent approaches to simulate the relationship between equivalent surface discontinuity size and the fatigue limit, which is the Kitagawa curve in other terms. 3.1 The Defect Stress Gradient Method 3.2 The Murakami Approach The Murakami approach [20] provides a good estimate of the fatigue limit for a large number of different defective materials. It has undergone several smaller developments over the years, which has led to robust applicability. The fatigue limit ow^ can be expressed for a given equivalent defect size -Jarea and load ratio R: The DSG methodology from Vincent et al. [19] in previous work [Leopold G, Nadot Y. J ASTM Int 7:2010], to represent the effect of a defect in the fatigue criterion by means of a stress gradient term. This general methodology called Defect Stress Gradient (DSG describes the effect of surface discontinuities on the high-cycle fatigue strength under multiaxial loading conditions. The local hotspot stress is modified with a defect and material-type dependent stress gradient function. In the current analysis, the approach is applied with the Crossland equivalent stress, characterizing the multiaxial fatigue strength of the defect-free material, which limits the application to proportional loading conditions. The Crossland parameters have been identified on the fatigue limits for three different loading types based on the experimental results of Nadot et al. [1] on defect-free NCI, and are aCr = 1.13 and fiCr=255 MPa. The av parameter describes the effect of the defect on the multiaxial fatigue strength for spherical surface defects and is av = 209 ^m. The equivalent DSG stress aDSG predicts crack initiation at the value fiCr, and is calculated with the equation: OY max 7Cr _max _0 oCr -oCr va ißo (1) equivalent stress at the tip of the defect, aCr is the stress at the same location without the stress concentration effect of the defect. The Crossland equivalent stress is calculated with the following equation for a representative load cycle: °Cr = J a +aCr°h,max. (2) With the simplification, that the ratio a™ /a0Cr is equal to the theoretical elastic stress concentration factor for hemispherical surface defects in an infinite body K?ph = 2.06, the allowable hemispherical ' I-DSG surface defect size \JareaalLsph can be expressed according to the following equation: yja O I all.sph ( -1) (3) o = F ,(HV +120 )((area^j 1 (1 - * ) 2 (4) where Fioc is the location factor and is 1.43 for a surface defect; a is a material parameter and is given as 0.371 + HV- 10-4 for NCI material. The calculation process can be used to express -Murakami area allsc with short crack considerations: V area aiis Ftoc (Hv +120) 0-4 2 a i f (5) HV = 200 hardness was considered for the investigated NCI material as the main input parameter of the Murakami approach. A differentiation is given for Varea > 1000 ^m, where the equivalent crack is considered to be microstructurally long. In this regime, the threshold stress intensity factor (SIF) range for crack propagation is a constant value, estimated as: AK'h = 3.3 -10-3 • (HV +120) -V10003 (1 - * ) 2 • (6) This leads to the following equation for the equivalent defect size 4areaalUc considerations: V —Murakami 1 area all.ic =- AK' v YAaiff y with long crack y . (7) 7 Kspha0 - B "-t uCr Per 3.3 Fatigue Limit Prediction with Linear Elastic Fracture Mechanics Considering complete equivalence between surface defects and cracks of equivalent size, linear elastic fracture mechanics (LEFM) can be applied to model the state of the crack propagation threshold, in other Murakami = a terms the fatigue limit. The NASGRO crack growth model by Forman and Mettu [21] describes the phenomenon of crack closure with a crack-opening function from Newman [22]. The NASGRO model for the NCI material of interest has been calibrated on the experimental results of Bjorkblad [23]. The threshold stress intensity range in the NASGRO model is calculated as follows: AKfl =AK o a + a. 1 - f (8) (1 - A )(1 - R) where a = Varea is the crack length, Y = 2/n is the shape factor at the bottom of a hemispherical surface crack embedded in an infinite body under tension. The material parameters considered in the calculations are summarized in Table 5. The allowable equivalent i-LEFM defect size ylareaatt according to the LEFM approach is given by: 4 area aU = — n AK, Y Act i f (9) Table 5. Material properties for the NASGRO model for ISO1083/ JS/500-7, based on experimental data from [23] Symbol Value Unit ao 0.038 mm a 2.5 - Smax / °0 0.3 - AKo 7.5 MPaVm Cth+ 1 - 3.4 Fatigue Limit Prediction with The Effective Threshold The equivalence of surface defects and cracks with the same size considered in the LEFM approach is, however, not fully correct. Long fatigue cracks have plastic wake along the crack path, whereas a crack initiated from a defect is a short crack with different parameters to consider in a crack propagation analysis. From a theoretical standpoint, cracks initiated from defects should be modelled with the equivalent size of the defect with the application of the threshold for microstructural crack initiation AKtheff. The effective threshold SIF range is independent of the load ratio. The equivalent allowable defect size can be estimated with the following equation: I-iKtkf 1 ( ^ s! area all 1 th.eff Y Act I.eff (10) AK thef = 3. 75 MPa>/m based on the where experimental results of Nadot et al. [24]. 3.5 Comparison with Experimental Results Figs. 8 to 10 show the comparison between the simulated curves with the different methods and the experimental results. The results obtained on specimens with complex artificial defects are displayed with the mean fatigue strength at cycles with a scatter band of three times the standard deviation. Some of the experimental data used for comparison originate from the work of Nadot et al. [1]. The fracture mechanics method based on the effective threshold SIF range (AKthef) leads to highly conservative results compared to specimens cycled until complete fracture. This approach models the start of microcrack propagation, hence it neglects parts of both the crack initiation and propagation phases. The fracture mechanics method utilizing the threshold SIF range (AKth) can lead to non-conservative results in cases of small defects, but overall describes the Kitagawa relationship well at different load ratios. This approach considers complete equivalence between surface defects and fatigue cracks, which can be a non-conservative assumption since cracks initiating from defects do not have the same closure levels as fully formed fatigue cracks of the same size do. The method of Murakami has a good general description of the Kitagawa-curve for ferritic-pearlitic NCI. This approach, however, should be applied carefully, since the utilized HV-AKth eJf correlation does not necessarily stand for the high-strength ferritic NCI grades. The DSG approach also models the multiaxial fatigue behavior of the defect-free and defective material in the HCF range as well, and therefore more general than the other analysed methods do. The curves recommended for design have been derived from the curves obtained with the DSG approach with the application of a safety factor estimated from the scatter of the defect-free results for this material grade. The mean stress-dependent safety factor is essentially the distance of the mean experimental and design curves on the Haigh diagram in Fig. 7. Its value is 1.64 for R_b 1.58 for R01, and R05 due to the smaller level of scatter in strength under higher mean stresses. 4 ALLOWABLE DISCONTINUITY SIZE MAPS The derived formulas for the allowable surface discontinuity size can be used to plot corresponding FE-results for a general complex case. Within the DSG approach the Crossland equivalent stress is employed to describe multiaxial loading effects. For the Murakami and the fracture mechanics-based a approaches, the maximum principal stress range can be used to describe multiaxial stress states, since that corresponds to the crack opening mode in fracture mechanics. It has to be noted, that these criteria are only applicable to proportional loading conditions. A simple analysis example is employed to demonstrate the allowable discontinuity size FE-results obtainable by the methods introduced in Chapter 3. Fig. 11 shows the analysis setup for a rectangular 3D body under pure bending. The R0 cyclic load was set up with the aim of inducing bending stresses leading to first yield according to the standard proof strength of 320 MPa at the maximum load. E = 169 GPa and u = 0.275 were applied to set up the linear elastic material model for the FE-solution. Fig. 8. Kitagawa curves for R_j tension-compression for NCI with artificial surface defects, with additional experimental results from [1] The simulations were conducted with ANSYS Academic Research 17.2 software. Fig. 12 shows the allowable discontinuity size maps obtained by the different approaches. The Fig. 9. Kitagawa curves for R01 tension for NCI with artificial surface defects, with additional experimental results from [1] Fig. 10. Kitagawa curves for NCI under R0.5 tension loading with artificial surface defects Fig. 11. a) FEA model setup and b) the normal stress [MPa] result field at the maximal load Fig. 12. Allowable discontinuity size FE-result fields with a) Murakami's, b) effective threshold, c) DSG and d) LEFM approaches differences of the compared approaches follow the trend discussed in Chapter 3.5. Fig. 13 displays the design values of the allowable surface discontinuity size derived from the results obtained with the DSG approach with the application of the mean stress-dependent safety factor. Illustrations of the corresponding discontinuity sizes are linked with the different colour bands to provide some level of physical connection. From the illustration of the proposed methodology, it can be concluded that even large surface defects can be tolerable from the standpoint of safe in-service life, if their position on the component is carefully monitored. In the current case displayed on Fig. 13, 1 mm to 1.5 mm surface defects can be allowed on more than two thirds of the component surface, which is marked by blue on the defect size map. (The legend indicates 1.5 mm to 2 mm for this colour band, but the distribution inside one band is unknown. For this reason, the allowed defect size range is shifted to a lower class.) The area-ratio of the different colour bands is governed by the nature of the Kitagawa diagram, which shows a decreasing sensitivity (especially the DSG approach) to the increase of the defect size in the investigated range. The critical area, where no discontinuities are allowed, must be generally quite narrow, otherwise the component would not even comply with the traditional fatigue assessment. Fig. 13. Allowable discontinuity size FE-result field for design with the DSG approach Allowable surface defect size maps have been plotted for a simple analysis example: for the comparison of the methods in their planned field of application. Surface discontinuities of 1000 ^m can reduce the fatigue strength of NCI by a factor of 2.5. (hemispherical defects in cylindrical test specimens under R_j tension-compression [1]). The fracture mechanics method based on the effective threshold SIF range (AKtfief) leads to highly conservative results compared to specimens cycled till complete fracture. The fracture mechanics method utilizing the threshold SIF range (AKth) can lead to non-conservative results in case of small defects, but overall describes the Kitagawa relationship well at different load ratios. The method of Murakami has a good general description of the Kitagawa-curve for ferritic-pearlitic NCI. This approach however should be applied carefully, since the utilized HV-AKth correlation does not necessarily stand for the high-strength ferritic NCI grades. The DSG approach also models the multiaxial fatigue behaviour of the defect free and defective material in the HCF range and is, therefore, more general than the other analysed methods. The fatigue properties recommended by the FKM Guideline [16] are in good agreement with the tested low-strength ISO1083/JS/500-7 material; they truly represent the lower end mechanical properties for this grade. 6 ACKNOWLEDGEMENTS The recent study was realized within the Knorr-Bremse Scholarship Program supported by the Knorr-Bremse Rail Systems Budapest. 7 NOMENCLATURES 5 CONCLUSIONS Fatigue tests have been conducted to analyse the scatter of fatigue properties for the ISO1083/JS/500-7 NCI material grade. A 40 % difference in the fatigue strength (aR00 05 ) is possible between polished samples machined from 500-7 NCI components. The step-by-step [2] approach has proven to give good results for the analysed material; the loading history below the fatigue limit had a negligible effect on the fatigue strength at load cycles. Different methods applicable to the assessment of surface discontinuities have been compared in the current paper. a material parameter in the DSG criterion, [M-m] a0 material constant depending on the microstructure, [mm] A component of the crack-opening function by •J area Newman, [-] equivalent defect size from Murakami, [MPa] Cth+ fitting parameter describing the .-dependence of AKth, [-] E modulus of elasticity, [GPa] f crack-opening function by Newman, [-] Floe location factor, [MPa] Hv Vickers hardness, [MPa] J2 a R R p0.2 R ^m Smax / °0 a aCr Par AKo AK th.eff amx _106 ^Cr max _arnp eff V [1] [2] [3] amplitude of the square root of the second invariant of the deviatoric stress tensor in a load cycle, [MPa] fatigue load ratio, [-] monotonic proof strength at 0.2 % elongation, [MPa] ultimate tensile strength, [MPa] ratio of the maximum applied stress to the flow stress in the NASGRO model, [-] plane stress/plain strain constraint, [-] material parameter in the Crossland criterion, [-] material parameter in the Crossland criterion, [MPa] the value of AKth at R0, [MPaVm] crack propagation threshold stress intensity factor range, [MPaVm] effective value of the threshold stress intensity factor, [MPaVm] applied cyclic stress max., [MPa] fatigue strength in stress amplitude at 106 cycles, [MPa] fatigue strength in stress amplitude at 106 cycles and load ratio "X", [MPa] Crossland equivalent stress, [MPa] max. of the hydrostatic stress in a load cycle, [MPa] fatigue limit for a given load-ratio, [MPa] the effective maximum principal stress amplitude, [MPa] the effective maximum principal stress range, [MPa] Poisson's ratio, [-] 8 REFERENCES Nadot, Y., Mendez, J., Ranganathan, N., (2004). Influence of casting defects on the fatigue limit of nodular cast iron. International Journal of Fatigue, vol. 26, no. 3, p. 311-319, DOI:10.1016/S0142-1123(03)00141-5. Bellows, R.S., Muju, S., and Nicholas, T. (1999). Validation of the step test method for generating haigh diagrams for Ti -6Al - 4V. International Journal of Fatigue, vol. 21, no. 7, p. 687-697, DOI:10.1016/S0142-1123(99)00032-8. Cheng, Z., Liao, R., Lu, W., Wang, D. (2017). Fatigue notch factors prediction of rough specimen by the theory of critical distance. International Journal of Fatigue, vol. 104, p. 195205, DOI:10.1016/j.ijfatigue.2017.07.004. Rotella, A., Nadot, Y., Augustin, R., Piellard, M., LHeritier, S. (2017). Defect size map for cast A357-T6 component under multiaxial fatigue loading using the defect stress gradient (DSG) criterion. Engineering Fracture Mechanics, vol. 174, p. 227-242, DOI:10.1016/j.engfracmech.2016.12.008. [5] Vincent, M., Nadot, Y., Nadot-Martin, C., Dragon, A. (2016). Interaction between a surface defect and grain size under high cycle fatigue loading: Experimental approach for Armco Iron. International Journal of Fatigue, vol. 87, p. 81-90, DOI:10.1016/j.ijfatigue.2016.01.013. [6] Schonbauer, B.., Yanase, K., Endo, M. (2017). Influences of small defects on torsional fatigue limit of 17-4PH stainless steel. International Journal of Fatigue, vol. 100, part 2, p. 540548, DOI:10.1016/j.ijfatigue.2016.12.021. [7] Mukherjee, K., Faster, S., Hansen, N., Dahl, A.B., Gundlach, C., Frandsen, J.O., Sturlason, A. (2017). Graphite nodules in fatigue-tested cast iron characterized in 2D and 3D. Materials Characterization, vol. 129, p. 169-178, DOI:10.1016/j. matchar.2017.04.024. [8] Ostash, O.P., Chepil, R.V., Vira, V.V. (2017). On the assessment of fatigue life of notched structural components. International Journal of Fatigue, vol. 105, p. 305-311, DOI:10.1016/j. ijfatigue.2017.09.007. [9] Gates, N., Fatemi, A. (2016). Notch deformation and stress gradient effects in multiaxial fatigue. Theoretical and Applied Fracture Mechanics, vol. 84, p. 3-25, DOI:10.1016/j. tafmec.2016.02.005. [10] Liu, J., Zhang, R., Wei, Y., Lang, Sh. (2017). A new method for estimating fatigue life of notched specimen. Theoretical and Applied Fracture Mechanics, DOI:10.1016/j. tafmec.2017.07.017. [11] ISO 1083:2004(E). Spheroidal Graphite Cast Irons -Classification. International Organization for Standardization, Geneva, p. 31. [12] Murakami, Y. (2002). Metal Fatigue: Effects of Small Defects and Nonmetallic Inclusions, Elsevier Science, Oxford. [13] ISO 945-1:2017. Microstructure of Cast Irons - Part 1: Graphite Classification by Visual Analysis, International Organization for Standardization, Geneva. [14] DIN 50125:2016-12. Testing of metallic materials - Tensile test pieces, DIN, Berlin. [15] ISO 6892-1:2016. Metallic Materials - Tensile Testing - Part 1: Method of Test at Room Temperature, International Organization for Standardization, Geneva. [16] FKM (2012). Analytical Strength Assesment of Components, VDMA Verlag, Frankfurt am Main. [17] Yamabe, J., Kobayashi, M. (2006). Influence of casting surfaces on fatigue strength of ductile cast iron. Fatigue & Fracture of Engineering Materials & Structures, vol. 29, no. 6, p. 403-415, DOI:10.1111/j.1460-2695.2006.01017.x. [18] Rabb, B.R. (2004). Influence of Occasional Underloads on Fatigue. European Congress on Computational Methods in Applied Science and Engineering, Jyvaskyla, p. 229-250. [19] Vincent, M., Nadot-Martin, C., Nadot, Y., Dragon, A. (2014). Fatigue from defect under multiaxial loading: defect stress gradient (DSG) approach using ellipsoidal equivalent inclusion method. International Journal of Fatigue, vol. 59, p. 176-187, DOI:10.1016/j.ijfatigue.2013.08.027. [20] Murakami, Y. (2002). Mechanism of Fatigue in the Absence of Defects and Inclusions. Metal Fatigue: Effects of Small Defects and Nonmetallic Inclusions. Elsevier Science Oxford, DOI:10.1016/B978-008044064-4/50001-3. [21] Forman, G.R., Mettu, S.R. (1990). Behavior of surface and corner cracks subjected to tensile and bending loads in Ti-6A1-4V Alloy, NASA Technical Memorandum, NASA, Houston. [22] Newman, J.C. Jr. (1984). A crack opening stress equation for fatigue crack growth. International Journal of Fracture, vol. 24, no. 4, p. R131-R135, DOI:10.1007/BF00020751. [23] Bjorkblad, A. (2005). Conventional vs. closure free crack growth in nodular iron. Samuelsson, J., Marquis, G. (eds.). Competent Design by Castings Conference: Improvements in a Nordic Project, p. 273-286. [24] Nadot, Y., Mendez, J., Ranganathan, N., Beranger, A.S. (1999). Fatigue life assessment of nodular cast iron containing casting defects. Fatigue & Fracture of Engineering Materials & Structures, vol. 22, no. 4, p. 289-300, DOI:10.1046/j.1460-2695.1999.00162.x. Strojniški vestnik - Journal of Mechanical Engineering 64(2018)6, 383-392 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2017.5123 Original Scientific Paper Received for review: 2017-12-01 Received revised form: 2018-03-29 Accepted for publication: 2018-04-26 Pneumatic Muscle Actuated Wrist Rehabilitation Equipment Based on the Fin Ray Principle Ioana Petre - Andrea Deaconescu - Flavius Sarbu - Tudor Deaconescu* Transylvania University of Brasov, Romania Rehabilitation of the hand entails the deployment of equipment that, by motion, accelerates the regaining of its specific functions and allow the swift social reintegration of patients at a minimum cost. In this context, the paper presents and discusses a pneumatic muscle actuated piece of equipment that allows the simultaneous rehabilitation of the radiocarpal, metacarpophalangeal and interphalangeal joints. The main novelty of this equipment consists in the solution selected for the support and mobilization of the hand; the construction put forward is based on a couple of bar mechanisms that form a Fin Ray type structure, specific to fish fins. The paper presents the constructive solution, the structure of the actuation system and the experimental results that highlight both the benefits and disadvantages of the proposed equipment. Of particular interest is the compliance of the system that contributes to its adaptability to the particular pain threshold of the patient. Keywords: compliance, fin ray effect, pneumatic muscles, wrist rehabilitation equipment Highlights • New bio-inspired wrist rehabilitation equipment is proposed. • A fin ray type structure is designed to support and mobilize the palm. • The structural diagram of the rehabilitation equipment, its construction, and its pneumatic control diagram are presented. • Attaining the angular motion limits of the hand joints and the variation of the rehabilitation equipment compliance are both demonstrated. 0 INTRODUCTION The wrist is an extremely complex joint, and its integrity is essential for everyday activity. The basic skeleton of the palm and the fingers includes 27 bones, eight of which (the carpal bones) form the most complex joint of the human body - the wrist. The actual hand consists of the five metacarpal bones followed by the 14 phalanges of the fingers. While the wrist ensures a wide range of motions in flexion, extension, circumduction, radial and ulnar deviation, its complexity renders it vulnerable to a series of affections likely to limit its functionality. Lesions of the wrist can be the cause of problems in the long term, ranging from articular stiffness, pain, numbness to partial or total disability. Rehabilitation by mobilization of the wrist has to commence as soon as physically possible, and the risk of relapse has been diminished significantly. The mobilization of the joints aims at correcting the retractions of the soft tissue, improvement of blood circulation, recovery of the gripping function, increase of muscle force and strength. Wrist rehabilitation can be achieved by a procedure known as continuous passive motion that entails setting the injured joint into motion by mechanized means, without straining the patient's muscles. This is possible by means of specially designed equipment, capable of applying the optimum rehabilitation motions to the wrist. As its primary control parameters, passive mobilization has the applied force, the stroke, the velocity, acceleration, duration and frequency of the displacement; all these quantities need to be adapted to the patient's clinical state and the set target. This means that equipment for passive mobilization has to allow the adjustment between certain limits of all the mentioned parameters. Equipment based on continuous passive motion currently available on the marketplace provides in addition to mobilizing the damaged joint also an assessment of muscle performance, thus allowing patient state monitoring during the entire period of recovery. Of the various manoeuvres possible as part of the mobilization exercises, slow progressive motion at varying amplitudes is most often applied. The force applied by the physical therapist at the maximum possible amplitude is typically dosed depending on the onset of pain, but also according to the therapist's experience with patients displaying either very high or very low pain thresholds. The necessity of rehabilitation of the human body joints by mobilization was documented as early as nearly 5000 years ago. Chinese documents, 4700 years old, held information about therapeutic postures and motions aimed at pain relief [1]. It was *Corr. Author's Address: Transilvania University, 29 Eroilor Bd., Brasov, Romania, tdeacon@unitbv.ro 383 in ancient Greece where the fundamentals of kinetic therapy were truly laid by Aristoteles, considered its founder. Later it was Hippocrates who noticed that motion contributes to regaining muscle force, while immobilization leads to atrophy. The utilization of wrist rehabilitation equipment has been a usual practice for centuries. One of the oldest achievements are the mechanical devices of Zander based on the pulley and weight mechanism built in 1857 [2]. In 1919, the Otto Bock Company put forward a series of rehabilitation devices designed for all joints of the upper and lower limbs [3]. At present several manufacturers of such equipment are present on the marketplace, like Kinetek, producer of Kinetec Maestra Hand & Wrist CPM [4], Remington Medical marketing the WaveFlex Hand [5] or Kinex with its Kinex Hand [6]. Of interest in the recovery of the wrist is also the system with three degrees of freedom that can be attached to the MIT-MANUS robot [7]. The limits of motion of this equipment are flexion/extension 60°/60°, abduction/adduction 30°/45°, pronation/ supination 70°/70°. Paper [8] presents a novel lightweight hand exoskeleton robot, which is called an advanced service robots (ASR) glove. The study shows that the proposed system can flex and extend fingers in their range of motion and grasp objects efficiently. The vast majority of rehabilitation equipment for continuous passive motion available on the marketplace at present is actuated by electric motors [9]. An analysis of the actuation systems of wrist rehabilitation equipment has revealed that about 72 % are electrically driven, 20 % pneumatically, 5 % hydraulically, and the remaining 3 % is covered by other energy sources [10]. At present, given the diversification, miniaturization and significant improvement of pneumatic components, the deployment of compressed air increasingly becomes a solution to be considered. Thus, in recent years. pneumatically actuated wrist rehabilitation equipment has been the object of several patents, a number of which are presented in [11] to [13]. Another example of pneumatically driven rehabilitation equipment is the parallel manipulator suitable for the complex motion of the human wrist joint [14]. A study of wrist rehabilitation equipment currently available on the marketplace has revealed that mobilization is achieved separately for each joint. It could also be observed that most equipment is actuated electrically, and pneumatic drives are neglected. Consequently, the authors of this paper put forward innovative equipment that mobilizes the wrist and finger joints at the same time, the actuation being ensured by a pneumatic muscle. Another novel aspect of the proposed solution is the bio-inspired construction of the palm support, based on the fin-ray effect, specific to fish fins. The benefit of the solution put forward consists in the fact that the complex motion performed by each of the hand segments, even though not all are affected, helps to improve their muscle and articular tone. Thus, in addition to the prophylaxis of joint pain achieved by the proposed equipment, it also helps to maintain and/or increase the muscle excitability of the healthy segments. The first part of the paper presents the biomechanics of the hand and fingers, highlighting the limits of the motions that have to be reached by the rehabilitation equipment. Furthermore, the actuation motor is presented, with a focus on a specific property of pneumatic muscles known as compliance. The second part describes the construction of the rehabilitation equipment and its actuation, followed by the third part in which the experimental results are presented. The paper concludes with a section of conclusions. 1 DESIGN REQUIREMENTS The development of rehabilitation equipment for the hand joints entails a thorough knowledge of the motions of the hand, of the forces and couples developed by it. It is further essential to identify the limits of these motions so that the new equipment can achieve them. With respect to this, the anatomic principles of the motions carried out by the joints of the hand are presented further on. 1.1 Biomechanics of the Joints of the Hand The possible motions of the joints of the hand are incredibly numerous and complex. A rehabilitation device capable of conducting all these motions would require extremely complicated kinematics and would be economically inefficient. For this reason, rehabilitation equipment available on the marketplace at present is typically designed for a single joint, for either the wrist or the finger joints. The proposed rehabilitation system discussed in this paper achieves mobilization by the continuous passive motion of the radiocarpal joint and at the same time of the metacarpophalangeal and interphalangeal joints. The proposed equipment can reproduce in the case of all mentioned joints the motions of flexion-extension (Fig. 1). Fig. 1. Flexion/extension of the wrist and joints of the fingers The flexion of the wrist or the palmar flexion is obtained by rotating the palm in the volar direction. Zero position is obtained when the forearm is in flexion at 90° and in pronation. In this case, the maximum rotation angle is of 80°. The extension of the wrist or dorsiflexion is obtained by rotating the palm in a dorsal direction, the maximum rotation angle being of 70°. In the case of the fingers, zero position is considered for the extended palm and fingers. For the metacarpophalangeal joints, the angle of flexion increases from finger II to finger V from 90° to 100°. The extension of the fingers varies from one individual to another and ranges from 0° to 90°, the latter value describing hyperextension of the fingers. Proximal flexion is the movement of the middle phalanx towards the inner palm, towards the proximal phalanx, until the limit of the motion is reached. The amplitude of this motion is of 100° to 120°. This flexion has small values for fingers I and II and greater ones for fingers IV and V. Distal flexion is the motion of the distal phalanx towards the inner palm, towards the middle phalanx. To obtain the maximum value, the proximal joint has to be in flexion. The maximum flexion value of this joint does not exceed 90°. The extension is possible only in the distal joints, its maximum amplitude being of 20°, and is encountered only in certain individuals. A larger extension angle causes the so-called hyperextension. Fig. 2 illustrates the main joints of the hands, of interest for the present study. Metacarpal Distal V interplialangeal joint i Proximal — ijCvj», ''""Uw lnterphalaiigeal >0» joint Metacarpophalangeal" joint Fig. 2. Joints of hand 1.2 Requirements of the Novel Rehabilitation Equipment The study conducted on hand rehabilitation equipment available on the marketplace has revealed specific characteristics of motions carried out by these, the values of which are given in Table 1. Table 1. Limits of the motions carried out by rehabilitation equipment of the wrist and fingers Joint Motion Rotation angle [°] Rotation speed [°/min] Wrist Flexion 0 to 80 Extension 0 to 70 Metacarpo- Flexion 0 to 90 Joints phalangeal Extension 0 to 90 of the Proximal Flexion 0 to 120 55 to 440 fingers Distal Flexion 0 to 90 Extension 0 to 90 In addition to the possibility of conducting the recovery motions within the limits resulting from studying the biomechanics of the joints of the hand, such rehabilitation equipment needs to satisfy requirements including: • compact construction and reduced weight such as to render the equipment portable to the patient (in the rehabilitation clinic or even at the patient's home). The portability of the entire rehabilitation system is assured by its construction, including the actual equipment (of dimensions 500 mm x 300 mm x 250 mm and approximately 3 kg mass) and the compressor (of dimensions 310 mm x 150 mm x 370 mm and 2.5 kg mass); • reduced activation time, entailing only the duration of placing the equipment on a table and connecting it to the electrical source; • utilization for its construction of robust, easy-maintenance components, typically manufactured from aluminium; • easy control and user-friendly interface. This facility requires only minimum operation proficiency of the user, given the easy command of the rehabilitation motions; • adaptive behaviour to the pain threshold of the person undergoing rehabilitation treatment. The latter requirement is achieved by using so-called adjustable compliance actuators (ACAs) for driving the equipment. A motor of this kind is selected given its benefits, such as its capacity of minimizing mechanical shocks, safe interaction with people or the ability to store and release energy in/from its passive-type elastic elements. The utilization of adjustable compliance motors can ensure the adaptability of the rehabilitation equipment to a given concrete working situation determined by the patient's pain threshold, which in some cases may be different than initially envisaged. If the recovery motion is opposed by a force generated, for example, by the increased rigidity of the joint to be rehabilitated, and implicitly by the pain felt by the patient, the pneumatic muscle will absorb (due to air compressibility) the induced shock, without forcing the movement. 2 CONSTRUCTION OF THE REHABILITATION EQUIPMENT 2.1 Constructive Principle Starting from the observation that maximum flexion of the wrist is obtained when flexion of the fingers is simultaneously induced, and the maximum extension of the wrist is achieved with the fingers being extended, the idea underlying the development of the novel rehabilitation equipment was bio-inspired, based on the Fin Ray effect, specific to fish fins. A pneumatic muscle was used as the motion generator element. Fish can move with the help of their fins, which ensure stability and propulsion. Similarities were observed between the movement of the palm and the fingers and the movement of fish fins, as the flexion-extension of the palm is similar to the oscillatory movement of the tail fin of fish (Fig. 3). In both cases, the bone structure is set into motion by antagonistic pairs of muscles located on each side of the skeleton. There are, nevertheless, differences between the movements of the palm and of the fishtail that need consideration when designing the construction of the palm rehabilitation equipment. One of these is the fact that the fishtails move symmetrically in relation to its relaxed resting position, while in the case of the palm and of the fingers the angles of flexion and extension, respectively, have different values; therefore, the movements of the fingers are not symmetrical in relation to the resting position. Another difference is connected to the movement amplitude of the fingers that is greater than that of the fishtail tip [15]. Fig. 3. Movements of the fishtail and of the palm Fish fins are elastic structures built around thin bone fibres forming an assembly that generates the Fin Ray effect. A mechanical Fin Ray effect-generating structure is obtained by several series-linked bar mechanisms. For the mechanism to operate similarly to a fin, the component mechanisms need to be arranged in a pyramid structure. These considerations underpinned the conceived rehabilitation equipment. Fig. 4 shows its structure next to that of a hand [16]. Fig. 4. Structural diagram of the developed equipment The assembly designed to support and mobilize the palm consists of two series-linked bar mechanisms forming a fin ray-type structure. The length of lever AB is given by the distance between the wrist and the metacarpophalangeal joints, while distance BE represents the length between the metacarpophalangeal and the proximal joints. Each of joints B and C that link the two bar mechanisms includes a torsion spring with a 270° angle between its arms; the role of these springs is to cancel a degree of mobility of the entire assembly. For this system, the maximum flexion of the palm is obtained for extended fingers, and maximum palm extension is achieved when flexion of the fingers is induced. Fig. 5 shows the positions of the hand on the mechanism for the maximum flexion and maximum extension of the wrist, respectively. Forearm H Metacarpals Proximal phalanges Intermediate phalanges Distal phalanges Fig. 5. Limit positions of the mechanism The first bar mechanism (ABCD) ensures the rotation of the palm. The counter-clockwise rotation of lever AB causes the extension of the palm, namely its lifting. At the same time, the second bar mechanism (BCFE) responsible for the movement of the finger joints will rotate clockwise. Thus, the counter-clockwise rotation of the palm will cause it to make a fist, while the clockwise rotation will open it. It can be noticed that this equipment provides (as an innovation) the simultaneous motion of the palm and finger joints, ensuring mobilization of several muscle sets at the same time. The degree of mobility of this mechanism, not considering the two torsion springs, is given by Eq. (1): k M = he • n e - r) = 3 • 6 - 8-(3 -1) = 2, (1) i=1 where ye is the dimension of the actual kinematic space, n the number of mobile elements, k the number of couples, fe the dimension of the space of relative velocities for couple k. Introducing the two springs linking element AB to BE and DC to CF, respectively, cancels a degree of mobility of the mechanism, as the movement of element BE becomes dependent on the movement of element AB. For bar mechanism BCFE to function correctly, the torsion springs mounted in couples B and C have to be twisted, thus tensioned over the entire working interval of this mechanism. The equipment for the rehabilitation of the hand joints includes two Fin Ray type mechanisms positioned in parallel and linked rigidly one to another, in order to create the necessary width for the positioning of the hand. The rehabilitation equipment is set into motion by joint A, by means of a pinion-rack mechanism with a pneumatic muscle acting as a linear motor that drives the rack to conduct a linear reciprocating motion. The pneumatic muscle is fed compressed air at its fixed end, while its mobile end is linked to the rack. Fig. 6 shows a photographic image if this rehabilitation equipment. Fig. 6. Construction of the rehabilitation equipment The pneumatic muscle contracts upon being fed compressed air and thus provides the required motor force. In this respect, in order to achieve the whole angle required for rotating the wrist in flexion/ extension, the inferior position of the mechanism is ensured by the muscle in relaxed state (not charged with compressed air), while its superior position can be obtained by the contraction of the muscle (at maximum feed pressure). A kinematic analysis of the rehabilitation equipment was achieved by means of the Mechanism modules of Pro/Engineer software. This analysis yielded the theoretical limits of the movements achievable by the rehabilitation equipment, featured in Table 2. It needs to be pointed out that the proposed rehabilitation equipment ensures the achievement of the movement limits required for the flexion and extension of the hand joints. Obtaining the different values of the flexion/extension angles is achieved by controlling the stroke of the pneumatic muscles, i.e. by feeding it compressed air at various pressures. Table 2. Limits of the movements achievable by the rehabilitation system Wrist [°] Metacarpophalangeal joints [°] Flexion Extension Flexion Extension Literature 80 70 90 to 100 0 to 0 Equipment 80 70 90 9.67 Fig. 7 shows the pneumatic actuation diagram of the rehabilitation equipment. The compressed air reaches the pneumatic muscle by means of a proportional pressure regulator (MPPES-3-1/4-6-010), controlled by a reference module MPZ-1-24DC-SGH-6-SW (all made by Festo, Germany). By means of rotational potentiometers, the reference module can generate up to six different values of the reference voltage, which are transmitted in the form of signals to the proportional regulator. If none of these reference values is used, the signal transmitted to the pressure regulator is a voltage adjustable via an external potentiometer. The continuous adjustment of the air pressure allows for modifying the compliance of the entire system. 2.2 Actuation of the Rehabilitation Equipment 3 EXPERIMENTAL RESULTS The dimensions of the pneumatic muscle utilized for the construction of the rehabilitation equipment were obtained upon conducting kinematic and dynamic calculations that yielded the magnitude of the required strokes and forces. Following these calculations, a DMSP-10-300N pneumatic muscle made by Festo was selected, ensuring a maximum displacement of 60 mm and the development of a required minimum force of 35 N. The first category of measurements was aimed at establishing a correlation between the values of the extension and flexion angles of the wrist achievable by the double bar mechanism and the feed pressure of the pneumatic muscle. For this, an E6A2-CS5C Omron rotary encoder was mounted on the axis of joint A. The horizontal line through joint A is considered the basis for measuring the angle of rotation. The angles for that support AB is above the horizontal axis are considered positive values (extension), while the negative values correspond to the position of support AB beneath the horizontal (flexion) (Fig. 8). Fig. 7. Control of the rehabilitation equipment by means of a proportional pressure regulator Fig. 8. Convention of signs for the flexion and extension angles The measured data points were obtained by continuously charging the pneumatic muscle with air at a pressure increasing from 0 bar to 6 bar. Upon reaching this level, the rotation angle was also measured for reducing the pressure from 6 bar to 0 bar. Table 3 features the measured values, and Fig. 9 shows the corresponding graph. The obtained experimental results revealed certain deviations of the limit angle magnitudes from those obtained theoretically (the data in Table 2). Thus, in the case of maximum flexion the deviation is of 1.8° at the beginning of muscle contraction, while in the case of maximum extension the deviation is of 2.7°. These deviations are explained by the effect of the two torsion springs that prevent reaching the ends of the stroke. In the absence of these springs, the mechanism achieves the limits found in the theoretical study. Another conclusion concerns the occurrence of hysteresis, caused mainly by the pneumatic muscle but also by the torsion springs. Table 3. Variation of the flexion/extension angles versus the feeding pressure of the pneumatic muscle Pressure [bar] 0 1 2 3 4 5 6 Rotation angle [°] -78.2 -64.7 -49.8 -25.0 7.1 35.0 67.3 Pressure [bar] Rotation angle [°] 67.3 57.2 45.7 21.3 -18.0 -49.0 -74. Wrist flexion/extension / -5G a ngle m " / ,J / f * —2(H W Exten sion / / ✓ / / t j ' w- t i Pressure [bar] -20 r t Flexion /' / / ' > / f / -70 ✓ ✓ J Fig. 9. The flexion/extension angle of the wrist versus pressure The main constructive element of Festo pneumatic muscles is a flexible tube made from chloroprene covered by a sealed envelope made from inelastic aramid fibres displayed in diamond patterns forming a 3D-mesh. The hysteresis of the pneumatic muscles is caused by the deformation of the flexible tube, but also by the internal friction between each aramid fibre and the elastic material enveloping it. The hysteresis of the pneumatic muscles increases the non-linearity of the rehabilitation equipment and consequently the complexity of the required control system. Further experimental research was aimed at establishing the link between the magnitude of the flexion/extension angle of the metacarpophalangeal joints (the angle between the extension of lever AB and BE) and the feed pressure of the pneumatic muscle. This angle was measured with a second rotary encoder, similar to the first one, placed in the axis of joint B. In this case the basis considered for the measurement was the axis of lever AB; the extension angles were considered positive and the flexion angles negative. The obtained results (Table 4 and Fig. 10) reveal the real maximum limits of the flexion and extension angles achievable on this equipment (maximum flexion of -87.6° and maximum extension of 8.3°). Table 4. Variation of the flexion/extension angles of the metacarpophalangeal joints versus the feed pressure of the pneumatic muscle Pressure [bar] 0 1 2 3 4 5 6 Rotation angle [°] 8.3 -13.7 -26.6 -39.4 -50.9 -68 -87.6 Pressure [bar] 6 5 4 3 2 1 0 Rotation angle [°] -87.6 -78.5 -68.2 -47.9 -30.2 -16 7.1 10; Metacarpophalangeal ^Sexion'extension angle [°] Extensic n Pres sure [bar] : ; Flextoi -10 -20 _ . -30 \ N\, \ ^ . v |Q \ X x V X 1 It * \ X Vv -60 ( \ \ 1 v \ V -70 V V > < X- ' \ X \ -90 Fig. 10. The flexion/extension angle of the metacarpophalangeal joint versus pressure The deviations in relation to the theoretical values shown in Table 2 are determined by the presence of the torsion springs, namely the same causes that prevented the attaining of the theoretical limits of the wrist rotation angle. In this case, the occurrence of hysteresis can also be noticed, caused by the behaviour of the pneumatic muscle. Another set of experiments was aimed at proving the compliance of the rehabilitation equipment. When 6 5 4 3 2 0 the dependency between the force generated by a motor and the displacement caused by it is of a nonlinear type, the stiffness of the actuator and implicitly its compliance are variable. This is the case of the pneumatic muscle that actuates the rehabilitation equipment. In order to determine the compliance of the rehabilitation equipment, the forces must be measured, developed by the pneumatic muscles in order to attain the various angular positions of the palm support. Measurements were conducted by means of a force sensor supplied by Festo, with the following characteristics: measuring range: 0 kN to 1 kN; supply voltage: 24 V DC; output voltage: 0 V to 10 V. The data provided by the sensor were processed by means of dedicated software FluidLab®-P V1.0, developed for the collection and processing of pneumatic system data. Considering a maximum force of 15 N given by the weight of the hand and of the mechanical elements of the mechanism that need to be rotated, and knowing the dimensions of the components of the actuation mechanism, the variation of the force required for the actuation of the mechanism is described by Eq. (2): l 35 F = G — cosa = 15---cosa = 35 • cosa, (2) r 15 where the notations are: G is weight of the mobile assembly, l distance from joint A to the application point of force G, r radius of the pinion pitch circle and a the flexion/extension angle of the wrist. Fig. 11 shows the graph corresponding to this equation. Fig. 11. Dependency of the required force to be generated by the pneumatic muscle vs the flexion/extension angle of the wrist From Eq. (2) and Fig. 11 result the necessary values of the strokes and forces to be developed by the pneumatic muscle to ensure various angular positions of the palm support (Table 5). Table 6 features the measurement results of the forces developed by the pneumatic muscle, and Fig. 12 shows a comparison of the necessary and the achieved forces of the rehabilitation equipment. Table 5. Required force of the pneumatic muscle a [°] 80 50 20 -10 -40 -70 Stroke s [mm] 0 5 10 15 20 25 F [N] 6.1 22.5 32.89 34.46 26.82 11.99 Table 6. The force developed by the pneumatic muscle Pressure [bar] 0 1 2 3 4 5 6 Stroke of the rack [mm] 0 4.1 8.3 12.5 16.6 20.8 25 Force developed by the muscle [N] 0 26.9 64.8 106.3 147.5 185.8 219.9 Force [N] 20 /fi** Stroke [mi 10 1 5 2 0 2 5 Fig. 12. Dependency of the force developed by the pneumatic muscle on its stroke From the graph, it can be observed that at the beginning of the lifting motion of the hand (counterclockwise rotation), the force developed by the muscle is smaller than required. This entails a delay of the motion until sufficient air pressure is reached such as to overcome this inconvenience. Nevertheless, this time lag of the rotation in relation to the starting time of feeding the pneumatic muscle with compressed air does not influence the efficiency of the rehabilitation exercise. The regression function corresponding to the curve that describes the evolution of the force developed by the pneumatic muscle is: F = 0.0233 • .5 • s - 3.969. (3) The correlation coefficient attached to this equation is of 0.998, which indicates that the proposed function describes the studied phenomenon with excellent precision. The stiffness of the analysed system (k) is calculated by Eq. (4), and its inverse, compliance (C), by Eq. (5): dF k =--= -0.0466 • s - 8.5, ds C = k=- 1 (4) (5) -0.0466 • s - 8.5 where 5 denotes the magnitude of the axial contraction of the pneumatic muscle as it is fed compressed air. Figs. 13 and 14 show the variation of the rehabilitation equipment stiffness and its compliance, respectively, as the pneumatic muscle shortens axially upon being fed compressed air. —6 k [N/mm] Stroke [mm -S 10 i 5 2 0 2(5 -to Fig. 13. Stiffness of the rehabilitation equipment 0.05 C [mm/N ] Stroke [mm] -0.05 -0.1 10 1 5 20 2 5 -0,15 Fig. 14. Compliance of the rehabilitation equipment The rehabilitation equipment proposed in this paper is characterized as visible in Figs. 13 and 14 by linearly decreasing stiffness and increasing compliance, as the stroke of the pneumatic muscle advances and the pressure grows. Consequences of such evolution of stiffness and compliance are a greater response time of the system to load variations, and evidently lesser precision, an aspect however of no significance to the efficiency of the rehabilitation exercises of the joints of the hand. On the other handin cases of rehabilitation exercises of the joints of the hand according to their degree of mobility, a compliant system like the one presented in this paper provides the benefit of adaptive behaviour to the concrete situation, allowing mobilization without causing pain to the patient. 4 CONCLUSIONS This paper presents a piece of pneumatic muscle actuated wrist rehabilitation equipment. Its novelty consists in the assembly developed for supporting and mobilizing the hand, based on two bar mechanisms that form a fin ray-type structure, specific to fish fins. A further innovative element is the utilization of a pneumatic muscle as the actuation motor of the proposed system. The conducted theoretical and experimental research has shown that the rehabilitation equipment ensures flexion and extension of the radiocarpal, metacarpophalangeal and interphalangeal joints that attain the angular limits of a healthy hand. Mobilization by means of a pneumatic muscle benefits from the compliance of the rehabilitation equipment, which contributes to the adaptability of the studied system to the individual concrete pain threshold of the patient. The disadvantage of using pneumatic muscles as actuators is the occurrence of hysteresis that diminishes the positioning accuracy of the mechanical assembly designed for the mobilization of the joints. Taking into consideration, however, that, in the case of rehabilitation of the hand joints, patient comfort takes precedence over positioning precision, the proposed and discussed equipment appears as a viable alternative to rehabilitation equipment available on the marketplace at present. Upon completion of the research related to the prototype, the authors of the paper intend to move to clinical tests in view of improving the actuation system of the equipment. 5 REFERENCES [1] Sbenghe, T. (2008). Kinesiology - Movement Science. Medical Publishing House, Bucharest. [2] Motet, D. (2010). Physical Therapy Encyclopedia. vol. 2, Semne Publishing House, Bucharest. [3] Otto bock Group (2017). from http://www.ottobock.com/en/, accessed on 2017-08-15. [4] Maestra (2017). Kinetec, from http://kinetec.fr/en/kinetec-selection/cpm-continuous-passive-motion/attelle-kinetec-maestra-detail.html, accessed on 2017-08-15. [5] Waveflex Hand CPM (2017). Remington Medical, from http:// www.remingtonmedical.com/product/detail/A1, accessed on 2017-08-15. [6] Kinex Hand (2017). Kinex CPM Devices, from https://www. kinexmedical.com/cpm-devices, accessed on 2017-08-15. [7] Krebs, H.I., Volpe, B.T., Williams, D., Celestino, J., Charles, S.K., Lynch, D., Hogan, N. (2007). Robot-aided neurorehabilitation: a robot for wrist rehabilitation. IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 15, no. 3, p. 327335, D0l:10.1109/TNSRE.2007.903899. [8] Hadi, A., Alipour, K., Kazeminasab, S., Elahinia, M. (2017). ASR glove: A wearable glove for hand assistance and rehabilitation using shape memory alloys. Journal of Intelligent Material Systems and Structures, vol. 29, no. 8, p. 1575-1585, D0I:10.1177/1045389X17742729. [9] Maciejasz, P., Eschweiler, J., Gerlach-Hahn, K., Jansen-Troy, A., Leonhardt, S. (2014). A survey on robotic devices for upper limb rehabilitation. Journal of NeuroEngineering and Rehabilitation, vol. 11, no.3, p. 1-29, D0I:10.1186/1743-0003-11-3. [10] Gopura, R.A.R.C., Bandara, D.S.V., Kiguchi, K., Mann, G.K.I. (2016). Developments in hardware systems of active upper-limb exoskeleton robots: A review. Robotics and Autonomous Systems, vol. 75, p. 203-220, D0l:10.1016/j. robot.2015.10.001. [11] Dongming, Y., Binrui, W., Yinglian, J. (2011). Pneumatic Muscle Flexible Elbow Joint Device with Buffer Spring and Flexible Shaft Sleeve. Patent CN202071080, China Institute of Metrology, Beijing. [12] Canjun, Y., Jiafan, Z., Jie, Z., Ying, C. (2007). Flexible Exoskeleton Elbow Joint Based on Pneumatic Power. Patent CN200984250, China Institute of Metrology, Beijing.. [13] Zhang, J., Liu, Y., Jia, Y., Shen, Z. (2012). Device for Hand Rehabilitation. Patent W02012120393, World Intellectual Property Organization, Geneva. [14] Takaiwa, M., Noritsugu, T. (2005). Development of wrist rehabilitation equipment using pneumatic parallel manipulator. Proceedings of IEEE International Conference on Robotics and Automation, p. 2302-2307, D0I:10.1109/ R0B0T.2005.1570456. [15] Deaconescu, T., Deaconescu, A. (2016). Applications of pneumatic muscles developed within the Festo national fluid actuation and automation training centre in Bra$ov. Recent Review, vol. 17, no. 4(50), p. 540-547. [16] Filip, 0., Deaconescu, T. (2017). Mathematical modelling of a Fin Ray type mechanism, used in the case of the wrist rehabilitation equipment. The 4th International Conference on Computing and Solutions in Manufacturing Engineering, vol. 94, D0I:10.1051/matecconf/20179401006. Strojniški vestnik - Journal of Mechanical Engineering 64(2018)6, 393-400 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2017.5189 Original Scientific Paper Received for review: 2017-12-30 Received revised form: 2018-02-19 Accepted for publication: 2018-04-16 New Hybrid System of Machine Learning and Statistical Pattern Recognition for a 3D Visibility Network Matej Babič1* - Karolj Skala2 - Dookhitram Kumar3 - Roman Šturm4 1 Jožef Stefan Institute, Ljubljana, Slovenia 2 Ruder Boškovic Institute, Croatia 3 University of Technology, Mauritius 4 University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Intelligent systems are an excellent tool to use for solving complex problems in the field of industrial applications. We use the mathematical method of fractal geometry and network theory when laser-hardening techniques are applied. The microstructure of the robot-laser-hardened specimens is very complex; however, we can present it by using a 3D visibility network. We convert the scanning electron microscope (SEM) images of the microstructure to a 3D graph and calculate the density of the visibility network of these 3D networks. We have analyzed the topographical properties of the hardened specimens by using the algorithm for the construction of a visibility network in a 3D space. We develop a new hybrid system of machine learning for predicting carbide content of the hardened specimens by using multiple regression, neural networks, and a genetic algorithm. We find the statistical significance of the relationship between attributes of the hardened specimens, the topological properties of visibility graphs, and carbide content of the hardened specimens._ Keywords: fractal geometry, hybrid system, laser hardened specimens, visibility network, statistical pattern recognition Highlights • We have calculated the statistical properties of the data of the parameters of the hardened specimens. • We have described the carbide content of the hardened specimens using the topological properties. • We have presented a new intelligent hybrid system model to predict the carbide content of the hardened specimens. • Our new method has many applications in pattern recognition, computer graphics, computational geometry, and so on. 0 INTRODUCTION Robot laser hardening [RLH] [1] is a heat treatment similar to inductive or conventional flame hardening. We can analyze the microstructures of RLH specimens using 3D visibility networks (graphs). A visibility network [2] is a graph of visible areas, which presents a set of nodes and obstacles in the Euclidean plane or space. Fractal geometry [3] was developed by Mandelbrot, who built the Mandelbrot set y=z2 + c with the help of a computer. He provided a new approach in the scientific discipline as he set out and designed a new way of thinking about structures and shapes. The Hurst parameter H [4] is the correlation between random steps X1 and X2, which is followed by the time-to-time difference At. Hurst parameter occurs in many areas of applied mathematics, including fractals and chaos theory, and is used in many fields ranging from biophysics to network computers. The parameter was originally developed in hydrology. However, modern techniques to estimate the Hurst parameter H come from fractal mathematics. The fractal dimension has been used to measure the roughness of sea coasts. The relationship between the fractal dimension D and the Hurst parameter H is given by the equation D = 2 - H for 2D objects and D=3 - H for 3D objects. We developed a new method to estimate the Hurst parameter H of a 3D object. In this paper, we introduce a new hybrid intelligent system to predict the carbide content of RLH specimens from the topological property of the density of the visibility network and fractal dimension. 1 PREPARATION OF MATERIAL SPECIMENS The study was done on the standard tool steel labeled as DIN 1.7225 [5]. The tool steel was surface hardened by laser at different speeds and different powers. We use a robot laser cell RV60-40 (Reis Robotics Company). The maximum power of the robot-laser cell is 3000 W. We hardened specimens with an output power of 1500 W. Therefore, we modified the speed parameter v £ [2, 5] mm/s and the temperature parameter T £ [800, 2000] °C. Each sample was prepared by etching and polishing (IMT, Institute of Metals and Technology Ljubljana, Slovenia) for a microscope evaluation (IJS, Jozef Stefan Institute). Fig. 1 was made by field emission scanning electron microscopy, JMS-7600F, JEOL. We wanted to know whether the microstructure of the RLH fractal patterns found a structure from which the Hurst parameter H could be estimated. Fig. 2 presents a 3D graph of the *Corr. Author's Address: Jožef Stefan Institute, Ljubljana, Slovenia, , babicster@gmail.com 393 microstructure of the RLH specimen. It is converted from Fig. 1 by using Fig. 1 color depth. Fig. 1. SEM picture of microstructure of the RLH specimen Fig. 2. 3D graph of microstructure of the RLH specimen 2 DESCRIPTION OF METHOD We used fractal geometry and a visibility network to determine the complex microstructure of the RLH specimens. Fig. 3 presents the random vertices of the 3D graph. For better visual presentation, we use 5x5 vertices of the 3D graph. 3D graph presents microstructure of the RLH specimens. Babic et al. [6] present a solution for constructing a visibility network in a 3D space. Fig. 4 presents the results of the problem of constructing a visibility network in a 3D space. We use the statistical topological property of the density of the visibility Fig. 3. Vertices in a 3D graph Fig. 4. Visibility network of vertices in a 3D graph network for pattern recognition from SEM images of the RLH specimens. The density q was calculated for each visibility network by Eq. (1): q = - 2m n x (n -1) (1) where m is the number of edges, n is the number of vertices in the visibility network. We present a method of estimating the Hurst exponent H for 3D objects [7]. First, we use the program Image J to find all the coordinates (x, y, z) of the SEM picture. Secondly, we estimate the Hurst exponent H by using the z-coordinates, which present a long continuous graph (Fig. 5). Also, all points (xz, y0, zz) present a space component on a 2D graph for all points (xz, zz). All points (xi,_y1,zJ) present a second space component on a 2D graph for all points (xj, zj). We made a space component for all yj, Vz. Then we combined all these space components into one space component. For this long space component, we can estimate the Hurst exponent H. We use the fractal dimension for pattern recognition from SEM images. 70 60 50 40 30 20 10 ■mil pi HiiiRiiMii 11 I III ¡1 Uli m i h mm.....mm i Iii11 I M. 'i'IL'II II,1 lU'Llll Inn n„i i mil h lull MIL. JIIILlUlULULIIIIIIil I 11 III nI I I I II 11IilRnR I I 'I II il imiij inn Ullillilililllllliii mini ujm tu jiiiu Ii if Iii,IIw Fig. 5. Long continuous graph We use the method of a visibility network and fractal geometry for statistical pattern recognition. To model the results, we use intelligent system methods, that is, multiple regression, neural network (NN), and genetic programming (GP). NNs [8] have the ability to solve a variety of problems. The sophistication of NNs is primarily due to their ability to imitate the principle of the functioning of the biological brain, which means that they solve problems similarly to humans. We used the Neuralyst program to create a model with NNs. Neuralyst is a software tool used within Excel. It has the ability to model NNs. We used a multi-tasking neural system with backpropagation and no back links. We had the option of setting different attributes. Table 1 presents the attributes of the NNS. Table 1. Attributes of the NN Learning speed [-] 0.6 Inertial coefficient [-] 0.5 Test mass tolerance [-] 0.02 Tolerance of the learning set [-] 0.03 Number of layers [-] 4 Fig. 6 presents a general multi-layer NN system. Table 2. Attributes of the GP Size of the population of organisms 500 Maximum number of generations 100 Reproduction probability 0.4 Crossover probability 0.6 Maximum permissible depth in the creation of the population 10 Maximum permissible depth after the operation of crossover of two organisms Smallest permissible depth of organisms in generating new organisms Tournament size used for selection of organisms Genetic programming [9] is similar to genetic algorithms and differs only in terms of the presentation method. Individual component in genetic algorithms is presented by a sequence of numbers, and the individual component in genetic programming is presented by a computer program. GP automatic writing of programs according to the nature of natural selection (evolution). At the beginning, we have some randomly written programs, which represent the initial Fig. 6. General multi-layer NN system 6 2 7 population. Then by crossing and selection, we get the next generation. Table 2 presents the attributes of the genetic programming. We used the genetic operations of reproduction and crossover. Fig. 7 presents an example of an organism in genetic programming. (22 -(^-))*(17+cos(Y)) Fig. 7. Organism in GP Multiple regression (MR) [10] is designed to investigate linear causes the relationship between a single dependent variable and one or more independent variables. With it, we determine the statistical feature in the power of connection and we predict the values of the dependent variable. The impact of each of the independent variables is estimated to be independent of the interactions between independent variables. Hybrid evolutionary computation [11] is a generic, flexible, robust, and versatile method for solving complex global optimization problems and could be used in practical applications. We present a new intelligent hybrid systems model in Fig. 8. Data information is presented in Table 3. Table. 3. Attributes of the hardened specimens S x1 x2 x3 x4 x5 x6 Y S1 1000 2 1.91 2.30 0.191 85 39 S2 1000 3 1.96 2.26 0.224 85 45 S3 1000 4 1.95 2.26 0.210 85 43 S4 1000 5 1.94 2.34 0.235 85 41 S5 1400 2 1.92 2.22 0.246 85 36 S6 1400 3 1.98 2.39 0.228 85 49 S7 1400 4 1.95 2.25 0.201 85 45 S8 1400 5 1.98 2.29 0.215 85 48 S9 1000 2 1.97 2.18 0.247 39 46 S10 1000 3 1.86 2.18 0.232 45 32 S11 1000 4 1.98 2.41 0.219 43 45 S12 1000 5 1.94 2.21 0.241 41 42 S13 1400 2 1.98 2.26 0.225 36 28 S14 1400 3 1.58 2.27 0.238 49 19 S15 1400 4 1.97 2.43 0.208 45 41 S16 1400 5 1.81 2.29 0.197 48 38 S17 800 0 1.97 2.23 0.289 85 47 S18 1400 0 1.98 2.24 0.277 85 52 S19 2000 0 1.97 2.26 0.245 85 50 S20 950 0 1.96 2.28 0.217 85 66 S21 850 0 1.95 2.32 0.212 85 80 S22 0 0 1.91 2.30 0.195 85 39 3 RESULTS AND DISCUSSION The attributes of the hardened specimens influence on the carbide content (Table 3). The specimens are labeled as S1 to S22. Attribute x1 represents the temperature [°C] and x2 represents the speed of RLH [mm/s]. Attributes x3, x4, and x5 represent the keys for pattern recognition. Parameter x3 represents the complexity in 2D, x4 represents the complexity in 3D, x5 represents the density of the visibility networks in a 3D space, and x6 represents the carbides in specimens. The last attribute Y is the measured surface carbide Fig. 8. Intelligent system model and new hybrid intelligent system content of the RLH specimens. Specimen S22 presents the material before the RLH process. In Table 4, we denote measurement (M) data (D) by MD, prediction data obtained by MR, prediction data obtained by neural network (NN), prediction data obtained by GP, and prediction data obtained by the hybrid system by H. Specimen S14 has a minimal carbide content after hardening, that is, 19 %. Table 4 present the statistical properties of the experimental and predicted patterns. Table 5 present topological properties of 3D visibility network. The measured and predicted carbide content of the RLH specimens is presented in the Fig. 9. The MR model is presented by Eq. (1), GP model is presented by Eq. (2), and hybrid model is presented by Eq. (3). We calculated precision of the GP model, NN and of the MR model by calculating average of absolute difference between measured and predicted data divided by measured data. The GP model has 85.57 % precision, the NN has 90.57 % precision, MR has 80.87 % precision and hybrid system has 62.98 % precision. Table 5 presents the topological properties of the 3D visibility network. Measured and predicted parts of carbides of the LHR specimens depend on attributes X], X2, X3, X4, X5 and x6 are presented in Fig. 9. We use the statistical topological property of the visibility networks to describe the carbide content in the microstructures. Image analysis of the SEM images of the RLH specimens is an interesting approach. With MR, GP and NN, we predict the carbide content in the microstructures. Finally, we present a new hybrid system of intelligent systems. For measured and predicted parts of carbides of the LHR specimens data, we calculated Kendall correlation coefficient. The best results for prediction give us NN, because the Kendall correlation coefficient (0.021) is most close to experimental data (0.131). Table 3 presents the topological properties of the 3D visibility network and attributes of RLH specimens. In this way, we can see how the attributes of speed and temperature influence the topological structures of visibility graphs in 3D space. Table 6 presents the statistical properties of the topological properties of the extreme number, number of edges, and triadic census type 16 to 300 of the 3D visibility network for RLH specimens. Firstly, we calculated the basic statistical properties of the mean, standard deviation, standard error, median, geometric mean, and harmonic mean of the topological properties of visibility graphs in 3D space of RLH specimens. We found significant positive relationships between the kurtosis, Fisher's G2, the coefficient of variation, the coefficient of dispersion, Table. 4. Measured and predicted data S MD MR NN GP H S1 39 55.28 38.98 38.7 34.77 S2 45 50.93 44.53 45.0 37.32 S3 43 47.84 42.25 45.0 40.61 S4 41 45.89 42.27 41.4 44.78 S5 36 42.41 36.98 40.8 44.18 S6 49 49.69 48.32 44.3 43.32 S7 45 40.97 44.86 46.7 46.45 S8 48 40.00 48.06 44.6 48.86 S9 46 42.74 45.97 46.2 40.85 S10 32 35.98 32.08 18.8 47.98 S11 45 48.01 44.98 40.0 44.99 S12 42 33.04 42.26 42.0 55.08 S13 28 39.91 46.48 27.7 47.05 S14 19 19.88 18.99 22.8 67.35 S15 41 42.15 41.54 45.2 51.32 S16 38 27.23 37.08 38.0 60.50 S17 47 59.78 47.86 47.4 30.50 S18 52 50.58 50.11 49.7 39.58 S19 50 42.10 51.01 50.1 48.54 S20 66 62.98 66.65 47.5 28.93 S21 80 65.82 79.25 46.4 27.35 S22 39 79.59 84.60 84.9 58.25 Table. 5. Topological properties of 3D visibility network S Extreme number Number of edges Triadic census type 16 to 300 S1 120823 3500351 4865624 S2 125787 3308776 4191425 S3 123943 3335861 4267175 S4 124833 3355735 4353872 S5 124626 3314397 4212248 S6 131540 3190001 3796016 S7 126962 3311163 4196282 S8 130799 3173601 3741603 S9 123393 3355056 4256560 S10 126395 3386391 4483986 S11 124296 3315948 4207031 S12 123829 3355735 4353872 S13 128143 3451450 4862060 S14 122500 3685175 5877473 S15 120818 3338595 4199754 S16 116812 3733624 5848517 S17 133031 3178192 3774789 S18 130974 3182544 3819193 S19 131043 3170121 3746658 S20 95090 4151533 7284078 S21 106916 5653616 1764141 S22 86871 5735036 1466536 Model of multiple regression: Y = -99.0509 - 0.0178 x x1 - 3.25717 x x2 + 46.65489 x x3 + 38.06208 x x4 - 63.3748 x x5 + 0.165097 x x6. (2) Model of genetic programming Y = 0.129654 x (-x6 - x22 - 2 + 2 + X X^ Xr :(-x 2 + x22) / \ X2 Xi + X2 X (X2 + X2 X X3 ) + X6 (3 X (x2 + x3) + X5 )X (3 + 2X X3 X (2 + x3) + x6)) (x2 + x3 + X6) X X 7.71283 + 2 x x2 - x6 — + - v X2 - X3 X V ^ y (3) X 6 X 6 2 X 6 4 6 X6 + — v 3 y X 6 + X 4 4 X6 + Hybrid model Y2 + Y (A + B) + AB - C — 0, A — 0.129654 ' ( —2 + —2 ) - X3 X3 + X2 B — 99.0509 + 0.0178x1 + 3.25717x2 - 46.65489x3 + 63.3748x5 -0.165097x6, CC — X2 X3 —6 , X1 + X2 (x2 + X2 X3 ) + X6 2 + -6 X2 — X3 (x3 (x2 + x3 ) + x6 )(x3 + 2x3 x (x2 + x3 ) + x6)) f \ X6 —6 + — (x3 X2 + X6 , f x6V —6 + ~ 7.71283 + 2x2 - x6 + 2 6 +- X2 X3 X2 X1 X 6 X 2 X x 6 + 90 80 70 60 50 40 30 20 10 0 A> $ & A* ^ 4> & & 4> 4> 4> > js o> o- o> o> o> o> o> o> ♦ Experimental data y = -0.1417X + 45.766 T = 0,131 ■ Prediction with MR y = -0.4789X + 51.998 T = - 0,08 A Prediction with NN y = -0.5965X + 53.91 T = 0,021 x Prediction with GP y = -0.4445X + 48.439 T = - 0.008 X Prediction with H y = -0.1457X + 46.61 T = -0,11 Density of the visibility networks in a 3D space of RLH specimens Fig. 9. Measured and predicted carbide content of the LHR specimens depend on attributes xh x2, x3, x4, x5 and x6 and the topological properties of visibility graphs in 3D space of RLH specimens. Table. 6. Statistical properties of topological properties of 3D visibility network SP Extreme number Number of edges Triadic census type 16 to 300 Mean 121792 3599223 4253132 Standard deviation 11531.35 714391.3 1201353 Standard error 2458.49 152308.8 256129.4 Median 124461 3346826 4209640 Geometric mean 121194 3546497 4055020 Harmonic mean 120509 3505535 3790200 Kurtosis 5.90 7.454659 4.719475 Fisher's G2 4.03 5.993685 2.517122 Coefficient of variation 0.094 0.1984849 0.2824632 Coefficient of dispersion 0.05 0.1024648 0.1701942 4 CONCLUSIONS The paper presents a new method of constructing visibility networks in a 3D space, a new method of describing the complexity of 3D space, and a new hyper-hybrid system of machine learning for the use in mechanical engineering to predict the topographical properties of materials. The paper presents a method of using visibility graphs in 3D space and fractal geometry to analyze the complexity of RLH specimens. Analyzing the complexity of RLH surfaces is a very hard problem. This new method has many applications in pattern recognition, computer graphics, computational geometry, and so on. The main findings are: 1. We use the method of network theory and fractal geometry to analyse the microstructure. 2. For prediction of the carbide content of hardened specimens, we use intelligent system methods, namely a neural network, multiple regression, and a genetic algorithm. The best results for prediction give us neural network. 3. We present the new hybrid spiral sequences. 4. The paper introduces a new method of machine learning in metallurgy. 5. We find the statistical significance of the relationship between attributes of the hardened specimens and the experimental and predicted pattern data. 6. The paper compares three methods, namely multiple regression, neural network and genetic programming, with a hybrid system of intelligent systems. 5 REFERENCES [1] Petrovič, S., D., Šturm, R. (2014). Fine-structured morphology of a silicon steel sheet after laser surface alloying of Sb powder. Strojniški vestnik - Journal of Mechanical Engineering, vol. 60, no. 1, p. 5-11, DOI:10.5545/sv-jme.2013.1347. [2] Ghosh, S.K. (1997). On recognizing and characterizing visibility graphs of simple polygons. Discrete & Computational Geometry, vol. 17, no. 2, p. 143-162, DOI:10.1007/ BF02770871. ISSN 0179-5376. [3] Mandelbrot, B.B. (1983). The Fractal Geometry of Nature. W.H. Freeman & Co., San Francisco, DOI:10.1119/1.13295. [4] Stoev, S., Pipiras, V., Taqqu, M.S. (2002). Estimation of the self-similarity parameter in linear fractional stable motion. Signal Processing, vol. 82, no. 12, p. 1873-1901, DOI:10.1016/ S0165-1684(02)00317-1. [5] Alta steel (2017). from http://www.altaspecialsteel.com/DIN-1-7225-AISI-4140-38.html, accessed on 2017-12-30. [6] Babič, M., Hluchy, L., Krammer, P., Matovič, B., Kumar, R., Kovač, P. (2017). New method for constructing a visibility graph-network in 3D space and new hybrid system of modeling. Journal of Computing and Informatics, vol. 36, no. 5, p. 1107-1126, DOI:10.4149/cai.2017.5.1107. [7] Babic, M., Kokol, P., Guid, N., Panjan, P. (2014). A new method for estimating the Hurst exponent H for 3D objects. Materials and Technology, vol. 48, no. 2, p. 203-208. [8] Graves, A., Schmidhuber, J. (2009). Offline Handwriting Recognition with Multidimensional Recurrent Neural Networks, in Bengio, Y., Schuurmans, D.,; Lafferty, J., Williams, C.K.I., Culotta, A. (eds.). Advances in Neural Information Processing Systems, vol. 22, p. 545-552. [9] Koza, J.R. (1992). Genetic Programming: On the Programming of Computers by Means of Natural Selection, MIT Press. [10] Armstrong, J.S. (2012). Illusions in regression analysis. International Journal of Forecasting, vol. 28, no. 3, p. 689694, DOI:10.1016/j.ijforecast.2012.02.001. [11] Ravi, V., Naveen, N., Pandey, M. (2013). Hybrid classification and regression models via particle swarm optimization auto associative neural network based nonlinear PCA. International Journal of Hybrid Intelligent Systems, vol. 10, p. 137-149, DOI:10.3233/HIS-130173. Strojniški vestnik - Journal of Mechanical Engineering 64(2018)6, 401-411 © 2018 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2017.5085 Original Scientific Paper Received for review: 2017-11-14 Received revised form: 2018-02-01 Accepted for publication: 2018-02-14 A Parametric Thermal Analysis of Triangular Fins for Improved Heat Transfer in Forced Convection Saroj Yadav - Krishna Murari Pandey* National Institute of Technology Silchar, India Thermal management systems in electronic devices require a reduction in size due to improve the overall performance of the system. The aim is to improve the heat transfer with reduction in weight of the system. Fins are the extended surfaces that ease the heat transfer process by increasing the wetted surface area. The thermal diffusion in a fin is always affected by parameters like the size, the shape, the material, the relative arrangement, orientation and position of the fins, the working fluid and its velocity, etc. Moreover, an extended surface may affect the pressure gradient in the flow domain. In this article, a three dimensional (3D) system of aluminium fin system has been numerically modelled. Simulations are done for conjugate heat transfer problem with fins of triangular shape. The thermal analysis is performed for various input parameters viz., arrangements, orientation and the number of fins. The effect of these parameters are analysed based on Nusselt number, convective heat transfer coefficient, coefficient of friction and coefficient of pressure. Keywords: conjugate heat transfer, extended surfaces, finite element method, heat transfer co-efficient, Nusselt number, triangular fins Highlights • The extended surface improves the heat transfer rate in any system by improving the thermal interaction between the solid surface and the fluid. • Heat transfer rate in a fin is always affected by the size, the shape, the material, the relative arrangement, the orientation and position of the fins, the working fluid and its velocity, etc. • Presence of an extended surface affects the pressure gradient in the flow domain. • Variation of Nusselt number, frictional coefficient, pressure coefficient and heat transfer co-efficient are affected by the arrangement, orientation and the number of fins. 0 INTRODUCTION Any system is under working condition generates heat. Sometimes, the amount of heat that it generates reaches the threshold limit of the system. Therefore, in order to avoid the damaging of the components of the system due to burning or overheating, it is very important to remove the heat in an effective way. Thus, heat transfer in enhanced rate is one of the important topics of thermal engineering. Enhancing the heat transfer rate with the controlled system dimension in the equipment like electronics, heat exchangers, aeronautical, internal combustion engine, etc. have always been a challenging task. To conquer this, researchers are continuously giving their effort in designing or reconstructing a new device with high thermal performance. These devices are intended for cost-effective and energy efficient transfer of heat in the compact areas without any thermodynamic loss. In this concern, an array of fins is one of the possible ways to increase the heat transfer rate from the base surfaces. This method can increase the heat transfer area that comes in contact with the fluid. Fin geometries and fin arrays can create turbulence in the fluid flowing, this further increase the heat transfer coefficient (h) between systems and surrounding and thus improve the performance. However, it gives rise to the pressure drop, which is a critical condition in most high performance applications. Consequently, it can be said that an optimized fin geometry and fin array combination is a conciliation of performance, pressure drop, weight, and size. The fin materials usually have a high thermal conductivity, thus it conducts the heat from the wall in high rate. Mostly, fins are used to enhance convective heat transfer by creating the temperature difference between the object and the environment. This in turn increases the convection heat transfer coefficient. Pin fin heat sinks are one of the competent heat exchanging devices used in many electric cooling appliances including central processing unit (CPU), transformer and thyristor. Thus most of the high thermal performance application based fins are made of copper or aluminium. Aluminium is always a preferred material for fins used in electronic system for cooling applications, due to its higher thermal conductivity and light weight. Researcher proposed many techniques through experimental, analytical and numerical analysis for the performance augmentation of extended surfaces, specifically pin fin heat exchanger. Pin-fin arrays induced turbulence in the flow field which further *Corr. Author's Address: National Institute of Technology Silchar, Department of Mechanical Engineering, Silchar, Assam, India, kmpandey2001@yahoo.com 401 helpful in enhancement of heat transfer. Axtmann et al. [1] focused their studies on thermally inactive pin-fins and presented the heat transfer results in terms of Nusselt number (Nu) and also discussed about other individualities like pressure drop and thermal performance parameter of investigated configurations. §ara [2], experimentally investigated the heat transfer, pressure and performing individualities for the array of staggered square pin fins attached on a flat surface in a rectangular duct and compared it with those for the inline arrangement. Jeng and Tzeng [3] also explore the heat transfer and the pressure drop in square pin-fins with variable inter-fin pitches for understandings of pin-fin arrays and systematic experiments has been performed, consequently the optimal inter-fin pitches are provided. Metzger et al. [4] presented variation in stream wise heat transfer, overall array heat transfer, and overall flow friction behavior for large aspect ratio ducts, that contained uniformly spaced staggered arrays of circular pin fins. Agarwal et al. [5] explored the effectiveness of pin fins for heat transfer in a channel cooled by air and agitated by a translational oscillating plate. Park et al. [6] designed a staggered pin-fin radial heat sink. The system was optimized for cooling of light-emitting diode (LED) device. They developed a numerical model to simulate various pin-fin array heat sinks. The results were verified experimentally. Huang et al. [7] estimated the optimal diameters for perforated pin fin array based on the desired temperature difference between base plate averaged temperature and ambient temperature and system pressure drop. The heat transfer and friction factor characteristics of the perforated rectangular fins were determined by Sahin and Demir [8]. The performances of three dimensional (3D) conjugate thermo-fluid analysis of micro pin-fins were presented by Abdoli et al. [9]. They considered a conventional circular shape, a hydrofoil shape, a modified hydrofoil shape and a symmetric convex lens shape of a fin. Joo and Kim [10] presented an analytical comparison of the thermal performances of optimized plate-fin and optimized pin-fin heat sinks under fixed volume condition. McNeil et al. [11] measured the heat-transfer coefficient and pressure drop measurements for a heat sink comprising micro pin-fins. In the following, a study has been performed on triangular fins with different geometrical parameters. The solver and the developed model are validated with the experimental data and with empirical relations available in the literature. Further, the variation of Nu, h, coefficient of friction (Cf) and coefficient of pressure (Cp) are studied by varying the arrangement, orientation and the number of fins. 1 FORMULATIONS AND GEOMETRY Thermal analyses of extended surfaces of different shapes are performed in this article. 1.1 Experimental Set-Up Experiments are carried out to understand the performance of a set of aluminium pipe bundle arranged in a staggered manner. Fig. 1a shows the experimental set-up used in the current study. The set-up consists of a vertical square duct with an axial flow fan placed just below the outlet temperature sensor. Air enters the duct from the bottom and leaves it from the top. The staggered arrangement of the pin fin system with the base plate as shown in Fig. 1b. The fins with the base plate are connected with the heater unit to supply heat from the bottom of the plate (Fig. 1b). The removable set-up (Fig. 1b) of fins with heating unit is placed inside the duct in the test section. The placement of the fins in the test section provides normal direction of flow to the incoming fluid (Fig. 1a). Fixed temperature sensors are provided b) a) Aluminum fin Base plate Fig. 1. The experimental setup of a) heat convection apparatus, and b) pipe bundle heating element at the inlet and the outlet of the duct to measure the temperatures of incoming and outgoing air, respectively. A portable thermocouple is also provided to measure the temperature at different locations of the fins. The display and control unit helps in controlling the speed of the fan and the amount of supplied heat to the heating unit. 1.3 Governing Equations Thermal behavior of any system is determined by the participated modes of heat transfer, and the thermo physical properties of the system. 1.3.1 Experimental Analysis 1.2 Computational Model Consideration is given to a 3D numerical model of aluminum pipe bundle placed in an air duct of dimension L*H*L (Fig. 2a). The fins of diameter D and height l, arranged in a staggered formation are placed over a base plate of dimension Lb^Lb^Hb (Fig. 2b). The arrangement is located inside the air duct at a distance (H - Lb) / 2 from the inlet. With a longitudinal and transverse pitch of SL and ST, the system is analyzed for a forced convection scenario. The considered numerical model is validated prior to the analysis of different cases. Air is driven inside the duct at an isothermal room temperature Tf and velocity V in the positive y direction. Numerically, the outlet is considered to be at atmospheric pressure (Po) and at outflow condition (n -VT = oj. The walls of the air duct are at no slip conditions (|V| = oj and at an isothermal room temperature Tf (Fig. 2a). The base plate of the fin system is supplied with a total uniform heat of Qs. Base plate a) Base Plate Inlet vS^ b) Fig. 2. Schematic diagram of the a) computational domain, and b) fin geometry In this study, the experimental analysis has been performed considering the force convection over the fin surfaces. The system is analysed based on the performance parameters like h and averaged Nu. The facility helps in measurement of various thermo physical properties of the system like system dimensions, the velocity and the temperature at various points. These quantities forms the basis for calculation of h and Nu in the system. For the current system (Fig. 1a), at various inlet flow velocities, with known values of temperature at inlet, outlet and over the fins, the h is calculated from the thermal energy balance. In ideal scenario, neglecting any heat loss from the back of the heater and from the duct walls (Fig. 1a), the rate of heat convection from the surface of the fin will be equal to the rate of heat carried by the incoming air. Mathematically, it can be expressed as, hA 'fin (( - Tn ) = mcp (Tout - Tn ), (1) where the Afin is the wetted surface area of the fin, the Tfin is the average surface temperature and m(= pAductV) is the mass flow rate inside the duct. Considering diameter D of the fin as the characteristic length, the average Nu is calculated using h L / kf. 1.3.2 Uncertainty in Experimental Data Every experiment needs analysis of uncertainty in the measured quantities. The uncertainty in the dependent derived quantities are due to various errors that appear in the measured quantities. In the present case, the uncertainty in the measurement of h and Nu are calculated using root mean square combination Eq. (2). In general, for any quantity A, it is expressed as E*=iE -—E , dx, (2) where E, is the uncertainty in the generic variable x. In the calculations of uncertainty of h, the generic variables are various areas, Aduct and Afin, velocity V, and temperatures Tin, Tout and Tfn. Similarly, for Nu, it will be h and D. 2 =1 1.3.3 Numerical Analysis The numerical study of the fin geometry requires combination of multiple physics of heat transfer and fluid flow to capture variation of temperature, velocity and pressure field in the domain (Fig. 2a). Considering forced convection over the present system, the conjugate heat transfer energy equation is expressed as, PCp if"+ V ■yT\ = V021 , J m \ -Ô.0166 1-0.0352 = -0.0414 15 ▼-0.0414 Fig. 11. Contour profile of y-component of velocity v [m/s] in staggered arrangement of fins in a plane at x = 0, with e = a), 0°, b) 30°, c) 45°, d) 60°, e) 90° c) %, 32.5 %, 44.8 % and 42.7 %, respectively, more frictional losses in the flow compared to the case with 0°. On the other hand, the same set of fins at different angles produces 50.4 %, 74.8 %, 90.3 % and 83 % more pressure losses in the domain compared to the fin system parallel to the flow, i.e., 0°. 2.3 Number of Fins With the understanding of heat transfer process in different arrangements of the fins, next the study has been extended to observe the effect of number of fins in the system. Four different cases with staggered arrangement of fins are considered with 2-1-2 (5 fins), 3-2-3 (8 fins), 4-3-4 (11 fins) and 5-4-5 (14 fins) staggered arrangement. The fins are place at angle 6 = 0° with the flow direction. A schematic diagram of the same is shown in Fig. 12. To maintain uniformity, the volume of all the fins and the aspect ratio (altitude: base) of individual fin is maintained the same as the earlier cases i.e. 2.5. Table 2 shows the results obtained from the numerical analysis for same boundary and thermal load conditions. It has been observed that as the numbers of fins are increased, the wetted surface area of the fins also increases, although the volume is remaining constant. However, the characteristic length of the fins decreases with increase in the number of fins. This affects the value of Nu and h. As the number of fins increases the Nu and the h decreases. Fig. 13 show variation of Cf and Cp for various configurations. As the number of fins increases, the wetted surface area of the domain increases. The increment in the surface area leads to more amount of shear stress experienced by the fin system. Hence, as the number of fins increases, the value of Cf also increases (Fig. 13a). However, the variation of Cp depicts a different scenario. As one moves from staggered 2-1-2 configuration to 3-2-3 configurations, the wetted surface area increases along with differential pressure. However, the relative rate of increment of pressure drop is found more than the increment in the wetted surface area in case of 2-1- a) b) c) Fig. 12. Schematic diagram of staggered arrangement of fins placed parallel to the flow direction with; a) 2-1-2; b) 3-2-3; c) 4-3-4; d) 5-4-5 configuration 2 configuration compared to 3-2-3 configurations. Hence, an increment in Cp is observed. As one move from 3-2-3 to 4-3-4 configuration, the ratio of pressure drop with respect to dynamic pressure is found to reduce and hence, the value of Cp decreases further. Table 2. Comparison of values of Nu and h for various numbers of fins placed parallel to the flow direction with staggered arrangements Cases Number of fins Nu h [W/(m2K)] Wetted surface area [mm2] Characteristic length [mm] 2-1-2 5 1.48 44.22 469.66 1.95 3-2-3 8 1.87 59.24 482.69 1.90 4-3-4 11 1.26 41.97 493.45 1.85 5-4-5 14 1.11 39.06 502.73 1.82 0.16 0.14 0.12 0.1 (J" 0.08 0.06 0.04 0.02 0 0.14369 0.12985 0.122 0.1073 a) 2 1.9 1.8 1.7 1.6 Y/E. By considering that e = y/R Y(c), 0 < y < c Y(y)n- c< y < h ' (25) The bending moment at any stage of the plastic deformation is [• h M = 2b oydy. Jo (26) By substituting Eq. (25) into Eq. (26) and integrating we obtain M = [3(Re)n - (1 -n)(R)2]. (27) M 2 + nL R Re The above equation applies for a non hardening material when n = 0 obtaining D ( M_ Re ^ J Me R I M < Me (3 - 2MM)-1/2- M > Me . (28) Fig. 7. The plastic deformation model. Considering the second partial derivative of Eq. (22) respect to x we obtain 1 R d2^ d x2. (29) The Fig. 7 illustrates a frontal view of the flange as a cantiliver beam with a terminal load, where l is the total length of the flange and 9 is the inclination angle. This assumption involves the following relations Fe = a F = l M Me l (30) where Fe is the load at the elastic/plastic boundary, a is the distance to the elastic/plastic boundary and Me is the bending moment at the elastic/plastic boundary. By inserting Eq. (29) and Eq. (30) into the Eq. (28) we obtain d2v_ J (3 - 2a)-1/2- x > a dx2 | x, x