Strojniški vestnik
Journal of Mechanical Engineering
no. 5
year 2021
volume
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: Demat d.o.o., printed in 250 copies Founders and Publishers
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
University of Maribor, Faculty of Mechanical Engineering, Slovenia
Association of Mechanical Engineers of Slovenia Chamber of Commerce and Industry of Slovenia, Metal Processing Industry Association
President of Publishing Council
Mitjan Kalin
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
Vice-President of Publishing Council
Bojan Dolšak
University of Maribor, Faculty of Mechanical Engineering, Slovenia
Cover:
The cover image shows a digital image correlation (DIC) measurement system, which can acquire heterogeneous strain fields (shown on the bottom). By using an inverse identification procedure, constitutive model parameters can be determined by comparing the measured strain fields with the calculated ones. However, prior to the procedure, proper spatial alignment of both data needs to be ensured.
Image courtesy:
University of Ljubljana,
Faculty of Mechanical Engineering,
Laboratory for Numerical Modelling and
Simulation, Slovenia
ISSN 0039-2480, ISSN 2536-2948 (online)
International Editorial Board
Kamil Arslan, Karabuk University, Turkey
Hafiz Muhammad Ali, King Fahd U. of Petroleum & Minerals, Saudi Arabia
Josep M. Bergada, Politechnical University of Catalonia, Spain
Anton Bergant, Litostroj Power, Slovenia
Miha Boltežar, University of Ljubljana, Slovenia
Filippo Cianetti, University of Perugia, Italy
Janez Diaci, University of Ljubljana, Slovenia
Anselmo Eduardo Diniz, State University of Campinas, Brazil
Igor Emri, University of Ljubljana, Slovenia
Imre Felde, Obuda University, Faculty of Informatics, Hungary
Imre Horvath, Delft University of Technology, The Netherlands
Aleš Hribernik, University of Maribor, Slovenia
Soichi Ibaraki, Kyoto University, Department of Micro Eng., Japan
Julius Kaplunov, Brunel University, West London, UK
Iyas Khader, Fraunhofer Institute for Mechanics of Materials, Germany
Jernej Klemenc, University of Ljubljana, Slovenia
Milan Kljajin, J.J. Strossmayer University of Osijek, Croatia
Peter Krajnik, Chalmers University of Technology, Sweden
Janez Kušar, University of Ljubljana, Slovenia
Gorazd Lojen, University of Maribor, Slovenia
Darko Lovrec, University of Maribor, Slovenia
Thomas Lubben, University of Bremen, Germany
George K. Nikas, KADMOS Engineering, UK
Tomaž Pepelnjak, University of Ljubljana, Slovenia
Vladimir Popovič, University of Belgrade, Serbia
Franci Pušavec, University of Ljubljana, Slovenia
Mohammad Reza Safaei, Florida International University, USA
Marco Sortino, University of Udine, Italy
Branko Vasič, University of Belgrade, Serbia
Arkady Voloshin, Lehigh University, Bethlehem, USA
General information
Strojniški vestnik - Journal of Mechanical Engineering is published in 11 issues per year (July and August is a double issue).
Institutional prices include print & online access: institutional subscription price and foreign subscription €100,00 (the price of a single issue is €10,00); general public subscription and student subscription €50,00 (the price of a single issue is €5,00). Prices are exclusive of tax. Delivery is included in the price. The recipient is responsible for paying any import duties or taxes. Legal title passes to the customer on dispatch by our distributor. Single issues from current and recent volumes are available at the current single-issue price. To order the journal, please complete the form on our website. For submissions, subscriptions and all other information please visit: http:// www.sv-jme.eu.
You can advertise on the inner and outer side of the back cover of the journal. The authors of the published papers are invited to send photos or pictures with short explanation for cover content.
We would like to thank the reviewers who have taken part in the peer-review process.
© 2021 Strojniški vestnik - Journal of Mechanical Engineering. All rights reserved. SV-JME is indexed / abstracted in: SCI-Expanded, Compendex, Inspec, ProQuest-CSA, SCOPUS, TEMA. The list of the remaining bases, in which SV-JME is indexed, is available on the website.
The journal is subsidized by Slovenian Research Agency.
Strojniški vestnik - Journal of Mechanical Engineering is available on
https://www.sv-jme.eu.
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5
Contents
Contents
Strojniški vestnik - Journal of Mechanical Engineering volume 67, (2021), number 5 Ljubljana, May 2021 ISSN 0039-2480
Published monthly
Papers
Andraž Maček, Janez Urevc, Miroslav Halilovič: Flat Specimen Shape Recognition Based on Full-
Field Optical Measurements and Registration Using Mapping Error Minimization Method	203
Abdulbasit Mohammed, Hirpa G. Lemu, Belete Sirahbizu: Determining Optimum Rotary Blade Design
for Wind-Powered Water-Pumping Systems for Local Selected Sites	214
Mingxing Han, Yinshui Liu, Yitao Liao, Shucai Wang: Investigation on the Modeling and Dynamic
Characteristics of a Novel Hydraulic Proportional Valve Driven by a Voice Coil Motor	223
Mohsen Soori, Mohammed Asmael: Virtual Minimization of Residual Stress and Deflection Error in the
Five-Axis Milling of Turbine Blades	235
Mario Santos, Jaime Santos: Ultrasonic Scattering Attenuation in Nodular Cast Iron: Experimental and
Simulation Studies	245
Ivan Okorn, Marko Nagode, Jernej Klemenc: Operating Performance of External Non-Involute Spur
and Helical Gears: A Review	256
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 203-213 © 2021 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2021.7111	Original Scientific Paper
Received for review: 2021-01-26
Received revised form: 2021-04-14
Accepted for publication: 2021-04-23
Flat Specimen Shape Recognition Based on Full-Field Optical Measurements and Registration Using Mapping Error Minimization Method
Andraž Maček - Janez Urevc - Miroslav Halilovič*
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
In the paper, an alignment methodology of finite element and full-field measurement data of planar specimens is presented. The alignment procedure represents an essential part of modern material response characterisation using heterogeneous strain-field specimens. The methodology addresses both the specimen recognition from a measurement's image and the alignment procedure and is designed to be applied on a single measurement system. This is essential for its practical application because both processes, shape recognition and alignment, must be performed only after the specimen is fully prepared for the digital image correlation (DIC) measurements (white background and black speckles) and placed into a testing machine. The specimen can be observed with a single camera or with a multicamera system. The robustness of the alignment method is presented on a treatment of a specimen with a metamaterial-like structure and compared with the well-known iterative closest point (ICP) algorithm. The performance of the methodology is also demonstrated on a real DIC application.
Keywords: full-field measurements, digital image correlation (DIC), specimen shape recognition, surface registration, iterative closest point (ICP)
Highlights
•	A methodology is proposed for aligning FEA data and the full-field measurement data of planar specimens.
•	The methodology offers both the shape recognition and alignment process to be performed on a single full-field measurement system.
•	The main advantage of the proposed methodology is its robustness.
•	The performance of the methodology is presented on synthetic as well as real DIC data.
•	The measured specimen can be observed with a single camera or a multi-camera system.
0 INTRODUCTION
Advanced constitutive models allow precise adjustment of the material mechanical response to specific loading conditions. Their flexibility is a result of a large number of free parameters. However, as the number of parameters increases, information gathered from standard experiments becomes insufficient or a huge amount of different experiments is required [1] and [2]. An alternative to the approach is taking into account the full-field kinematic information. Such data are normally acquired through digital images from where the displacement fields can be calculated using the digital image correlation (DIC) techniques. The approach gives rise to the development of complex specimen shapes [3], making it possible to identify material parameters through a single experiment [4]. One of the challenges, still not sufficiently resolved and addressed in this work, is the alignment of planar specimens between the modelling (numerical) data and experimental DIC data.
Full-field measurements offer a huge amount of information on the specimen surface (several 10'000), especially compared to classical extensometers.
However, their dependency on the sought material parameters generally cannot be explicitly determined, and therefore direct material identification cannot be executed. Researchers resort to the use of inverse identification techniques, where the finite element model updating (FEMU) [5] represents one of the established methods with extensive research work performed on the full-field measurements [3] (e.g. DIC [6] response calculation or by employing different optimization methods (e.g. genetic algorithm [7] or simulated annealing [8]). The method is based on an iterative comparison between the measured and calculated specimen's responses. The optimal values of material parameters are then determined by minimizing their discrepancy.
However, before finite element model (FEM) results and DIC data can be compared, coordinate systems of both sets of data need to be aligned. The latest research shows that such alignment is essential for reliable identification of material parameters in the case of specimens with non-homogeneous strain field where sharp strain gradients occur. Fehervary et al. [9] examined material parameter fitting results of planar tests when the sample orientation was unknown
*Corr. Author's Address: Faculty of Mechanical Engineering, University of Ljubljana, Aškerčeva 6, 1000 Ljubljana, Slovenia, miroslav.halilovic@fs.uni-lj.si
203
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 203-213
or deviated. The study showed that after a certain threshold of misalignment, reliable parameters can no longer be found. The authors also concluded that the level of threshold seemed to be material dependent. Lava et al. [10] studied two different methods used to compute full-field error maps between experimental DIC data and FEA data to, among others, investigate the effects of processing parameters, model form error (such as incorrect boundary conditions) and mesh alignment. The authors showed that even a small misalignment can have a surprisingly large effect on the strain error maps and exposed the necessity to develop robust and precise methods of alignment. Similarly, Ruybalid et al. [4] performed various virtual test cases to assess the performance of FEMU and integrated digital image correlation (IDIC) method when subjected to different error sources, among others, misalignment of the specimen. Both methods are shown to be sensitive to misalignment. Namely, the increase of the misalignment further increases the error on the identified material parameters.
It seems interesting that despite the importance of properly addressing the alignment of DIC and FEM data and considering the numerous publications addressing material characterization using planar specimens, not much information can be found on how researchers align both data sets in their works. However, as pointed out in a recent article from Polyga (a developer of 3D scanners with more than 10 years of experience), scanning flat objects can be particularly difficult even for experienced 3D scanning technician [11]. The process of aligning DIC in FEM data set generally consists of two steps, specimen's edge detection and alignment (point set registration). It is important to note that the procedure needs to be performed on the same optical system as DIC measurements so that the DIC and FEM data can then be compared, which adds to the complexity of the problem.
One of the best-known methods for point sets registration is the iterative closest point (ICP) algorithm introduced by Besl and McKey [12]. The algorithm consists of the closest point search and a minimization of the matching error, applied iteratively to the two surfaces to be matched [13]. Many variants of the algorithm have been introduced since, affecting the algorithm at different stages [14] and [15], e.g. the selection of points, matching, rate of convergence, etc. However, despite the widespread of the ICP method, its convergence in the general case of initial misalignment cannot be guaranteed [16]. Also, optical full-field displacement measuring methods most
commonly cannot accurately measure the specimen shape.
To address the problem of aligning full-field measurements with FEM data, Bruno et al. [17] utilized a linear transformation matrix to map the location of the calculation point into the measurement picture pixel position. The projection matrix was determined using user input coordinates of the calculation point and the pixel positions of three arbitrary points. The spatial position of all three points must be known a priori and all of them must be observed by the measurement. A similar approach was employed by Silva [18] who transformed the model coordinate system by recognizing that some reference points observed in the image have known numerical coordinates. Both calibration procedures are easy to implement but with the state of the art specimen designs, which exhibit smooth shapes and no dominant features [19] they become increasingly difficult to be used and prone to errors.
In the present work, a robust methodology is presented for aligning full-field measurements (e.g. DIC) and modelling (e.g. FEM) data of flat specimens. The procedure addresses both the contour recognition of a specimen as well as registration. Namely, for a successful DIC measurement, the specimen must be covered with a speckle pattern consisting of a white background and black speckles. In practice, this represents a problem because the distinction between the specimen and the background can quickly be lost. We determine the specimen's geometry by taking two consecutive images with changing the background illumination. From their comparison, the background can be subtracted, leaving just the specimen's geometry. The result is a black-white image of the measured specimen (white) and the background (black). To perform the alignment, the CAD model is also projected onto the image as black and white pixels. By minimizing the intersection between both images, i.e. the measured specimen and numerical model, both data sets are aligned.
In the following, the methodology is presented first in the case of 2D and then 3D measurements. The performance of the method is demonstrated on two examples. In the first example, a synthetic one, we demonstrated the method on a problem of aligning a specimen with a metamaterial-like structure. In the example, we also analyse the effect of misalignment and compare the performance of the method with the ICP. In the second example, the method is demonstrated on real DIC data where its ability to handle the presence of noise and mismatches in the geometry is also presented. Discussion of the
204
Macek, A. - Urevc, J. - Halilovic, M.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 203-213
results and the performance of the method is given in the Discussion section, which is followed by the conclusion.
1 METHODS
The method will be presented first in the case of 2D DIC measurements, where the specimen is observed with a single camera, located normal to its front surface (Fig. 1a). Afterwards, the method will be generalized for the case of 3D DIC measurement. Although in the latter case, additional coordinate transformations are needed to take into account the specimen's perspective (Fig. 1b), the alignment procedure is in both cases the same.
a)
b)
&

Fig. 1. Camera arrangement for a) 2D, and b) 3D DIC measurement
1.1 2D measuring System
Specimen recognition. Let us denote the numbering of pixels in the acquired (DIC) image in a form of a set:
P = {(u, v )e N x N |u < n and v < m|,
(1)
with the size of the image being n x m pixels.
To determine the set of pixels which represent the specimen surface, two images of the same specimen are needed. The first image is a regular one used for the measurement itself. The background on the acquired image is usually dark and cannot be clearly distinguished from the specimen. For the second image, we brighten up the background using an illumination panel, as shown schematically in Fig. 2 by the two left-most images (corresponding to specimen recognition).
The region of the specimen, which we wish to determine, is defined by a set Sm c P, which is a set of pixels, whose difference in grayscale values between both images is below a threshold value s:
Sm = {(u, v )e P||g (u, v)-g2 (u, v )|<e}, (2)
where gj(u,v) and g2(u,v) represent grayscale values of the (u, v)th pixel in the individual image, respectively.
Alignment. First, positions of pixels (u, v) of the specimen's shape (set Sm ) need to be transformed into a physical location associated with specimen dimensions (X, Y), see Fig. 2. This is performed with a calibration procedure, which can be written in a form of an arbitrary function c
c : N2 ^ R2; (u,v(X, Y).
(3)
By assuming that the position of the camera is normal to the specimen surface and neglecting the optical distortions [20], the calibration only scales (u, v) as
[ = c (u, v) = k <1 I, J	[n - u J
(4)
with k being the scaling factor.
The relation between coordinates (X, Y) and the modelling space coordinates (x,y) can be written in a form of a mapping function f
f :]

2; (X,Y(x,y),
(5)
where the mapping, in general, carries out translational and rotational rigid body transformation only. For a planar case, we can write the mapping as
= f (X, Y ) =
cos a - sin a sin a cos a
|Ax lAy
(6)
with a, Ax and Ay being the angle of rotation, translation in x and y direction, respectively. All three parameters (a, Ax, Ay) are unknowns and must be determined in such a manner that the modelling specimen, being in (x, y) space, coincides with the measured one, being in (X, Y) space.
So far, from the calibration and mapping procedure, one can determine each pixel location in the modelling space as (x, y) = f ◦ c) (u, v). If we denote the specimen region in the (x, y) space as (see Fig. 2), the set of pixels that are located inside the region denoted as a set Ss c P, follows as
={(«,v)e P| {f oc){u,v)

(7)
By comparing both images, the mismatch between the recognized measured and modelled specimen shape is determined as the symmetric difference between the two pixel sets Sm and Ss
E = S„ aS..
(8)
Thus, E contains all the pixels of both sets Sm and Ss except for the ones in their intersection. The number of those pixels, i.e. the cardinality of the error set |E|, represents the measure of the mismatch level.
X
Flat Specimen Shape Recognition Based on Full-Field Optical Measurements and Registration Using Mapping Error Minimization Method
205
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 203-213
CAD model (used for FEM)
Simulation shape
( ßs	Y)-
	
J
Alignment
V
i/Inverse map (a. Ax, Ay)
Methods' flowchart
Acquired £\ DIC data (Ax, Ay, a)-
Images
J_

Eg. (4) I 1 Eg. (6)


Eg. (2)
Eg. (7)

Eg. (8)
JK
I Eg (9) I I Eg. (11) I (Ax, Ay, a)
Iterations
Overlaid shapes
TJy\
Simulation shape
Pixel grid
' Iterations
Simulation shape
x <¡="
Pixel grid
Pixel grid
Fig. 2. The proposed alignment method
If the shapes are perfectly aligned, there will be no mismatched pixels, thus the cardinality of the error set would be zero.
The unknown mapping parameters (a, Ax, Ay) are determined as a solution to the following unconstrained optimization problem. Let us define the cost function CF (a, Ax, Ay) as
CF (a, Ax, Ay ) = |e| which needs to be minimized
min CF (a, Ax, Ay).
(9)
(10)
The result of the iterative optimization procedure (e.g. steepest descent method, Levenberg-Marquard algorithm) is a set of optimal values of mapping parameters ja, Ax, Ay J that minimize the cardinality of the error set
a, Ax, Ay} = arg min CF (a, Ax, Ay). (11)
J	{a, Ax ,Ay]
Once the optimal values of parameters are obtained, the measured results can be mapped into the modelling space via Eq. (6) as
■ = f (X ,Y ) =
cos a - sin a sin a cos a
jAx[
I* I
-. (12)
206
Macek, A. - Urevc, J. - Halilovic, M.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 203-213
1.2 3D Measuring System
The 3D measuring system consists of several digital cameras, each of them observes the specimen from its direction, as presented in Fig. 1b. The result of a 3D DIC measurement is a point cloud given in a measurement coordinate system (X, Y, Z), Fig. 3, which in general is not aligned with the specimen. In the case of planar specimens, a best-fit plane is normally constructed over the point cloud. The coefficients of the best-fit plane are determined using linear regression. Consequently, the location of the measured points in 3D (X, Y, Z) can be expressed with the 2D location on the best-fit plane (X, Y). The reduction to 2D is essential to enable the use of equations derived in the previous subsection.
Calibration procedure c, is performed for each camera with the purpose to map the pixels of the acquired images into the (X, Y) coordinate system. Each camera acquires its image with the corresponding pixel set P,, i = {1, 2, ..., TV} (analogously to Eq. (1)) with N being the number of cameras.
calculates kuYi, kvYi and nYi by considering the known pixels' positions (u, v) and the corresponding position on the Y axis of the best-fit plane.
The mapping f of measured points on the best-fit plane (X, Y) to the modelling space (x,y) is performed as presented in the 2D case via Eq. (6). Sets Smi, SSii and E, which again represent sets of pixels of the measured specimen, modelling specimen and the mismatch between the two, respectively, are constructed for each camera individually as already presented.
The level of the mismatch considering all cameras is expressed as the sum of all individual levels of mismatches squared
CF (, Ax, Ay ) = E;|
(14)
with CF (a, Ax, Ay) being the cost function. As in the case, . optimal values of mapping variables are determined by minimizing CF.
2D
ja, Ax, Ay J
Measured data are then mapped onto the modelling space via Eq. (12).
Measured pomts^v.'/VyN
I
Best-fit plane Coordinate System

Y
Measurement
X Coordinate System
Fig. 3. Best-fit plane
The calibration procedure for each camera ci is approximated using a linear transformation
IX
= ci (u, v) =
kuY,i kvY
(u, V)e P,
(13)
where k* and n* represent the constants of a two-dimensional linear function. Constants are obtained for each camera individually by employing two linear regressions. First regression calculates kuXi, kvXJ and nXJ by considering the known pixels' positions (u, v) and the corresponding position on the X axis of the best-fit plane. The second regression analogously
2 EXPERIMENTAL INVESTIGATION
The performance of the method is presented on two examples. For a practical purpose, both examples are presented in the following along with the corresponding results. Discussion of results is then provided in Section 3. Minimization of the cost function (Eq. (10)) was in both cases performed using the steepest descent method. Derivatives of the cost function were approximated using the finite difference method.
In the first example, the synthetic one, we demonstrated the method on a problem of aligning a specimen with a metamaterial-like structure. The purpose of the experiment is to present the advantage of the proposed alignment algorithm in comparison with the well-known ICP algorithm. The effect of misalignment on the performance of methods is also studied.
In the second example, the entire methodology is presented on real DIC data (i.e. shape recognition and alignment), where the ability of the method to handle the presence of noise and mismatches in the geometry is also tested. The example presents a heterogeneous strain field specimen proposed for the calibration of plastic anisotropy.
2
n
n
Flat Specimen Shape Recognition Based on Full-Field Optical Measurements and Registration Using Mapping Error Minimization Method
207
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 203-213
2.1 Synthetic Experiment
The geometry of the specimen is shown in Fig. 4. The shape consists of cut-out slots exhibiting an even pattern, mimicking the shape of metamaterial [21].
The specimen in Fig. 4 also represents the reference (modelling) geometry. The measured geometry is obtained by introducing rigid translation and rotation to the reference one [10].
0 0 0 0 _ . o o o oo
0 0 0 0 0 0 CD CD CD CD CD
0 0 0 0 0 0 CD CD CD CD CD
Fig. 4. Investigated synthetic specimen (dimensions in mm)
The effect of initial rotation and translation was examined to test the robustness of the method, mimicking multiple measurement setups. The translation was assumed equal in the horizontal x and vertical y direction (for the sake of brevity, to reduce the number of cases). Analysed cases include all possible combinations of rotations for angle a, ranging from 0° to 90° by a step 15°, and translations Ax=Ay, ranging from 0 mm to 2 mm with a step of 0.5 mm. An example of the initial misalignment (a = 2°, Ax=Ay=0.5 mm), presents Fig. 5a. The corresponding error image, i.e. pixels from the error set E, Eq. (8), is presented in Fig. 5b in white.
_
BEro
dgqQt'n
ojtctto)
Fig. 5. Example of the analysed case: a) initial misalignment between simulation and measurement, b) the corresponding error image E, Eq. (8)
Besides the proposed method, the ICP was also employed in this example for comparison. Although many different versions of the ICP exist, one of the most basic closed-form versions of the algorithm was
assumed. Data used for the ICP are those presented in Fig. 5a. The closest point search was performed using the grid closest point and no false matches rejection was assumed.
Results. The performance of both methods is demonstrated by presenting the number of iterations that were needed to align each initial misalignment.
Results of both methods are presented in Fig. 6. Initial rotational misalignment a is displayed on the vertical axis and the translational Ax, Ay on the horizontal axis. Colours represent the number of iterations used for a successful alignment. Cases where the alignment was not achieved (exceeding 100 iterations or converged to an inappropriate position) are displayed in white.
ICP
New method
	Used
iterations	
|	80
	60
	40
	20
1	■0
C	Not foui
0.5 1. 1.5 2.
0.5 1. 1.5 2.
Ax,Ay (mm)	Ax,Ay (mm)
Fig. 6. Performance and convergence of ICP and the proposed alignment method
The proposed method managed to align the specimen in almost all the analysed cases of the initial misalignment. On the other hand, the ICP failed to converge or converged to a wrong minimum at initial misalignment rotations greater than 45°. With the increase of the initial rotation, the number of iterations considerably increased for the ICP, reaching over 80 iterations. The proposed method needed approx. 15 iterations on average for a successful alignment, irrespective of the initial rotation.
2.2 Application to Real DIC Data
Practical application of the method is demonstrated on a specimen presented in [22] (Fig. 7), designed to induce a biaxial strain-stress state. Due to the specimen's diverse geometry and no distinct features (such as sharp corners or round holes), the problem also represents a good alignment test problem.
Measurements were performed with a DIC measuring system Q-400 Dantec Dynamics GmbH, (Ulm, Germany). The measuring setup is presented in Fig. 8. We utilized 4 digital cameras, which in pairs acquire images of the front and
208
Macek, A. - Urevc, J. - Halilovic, M.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 203-213
backside of the specimen, representing a case of a 3D DIC measurement. Such setup enables biplane measurement, where we simultaneously measure the front and back specimen surface in a common coordinate system. The calibration procedure was performed using a special two-sided calibration target. The purpose of such a set-up is enhanced characterization of material mechanical behaviour [22].
t = 0.68
70
Fig. 7. Investigated specimen (dimensions in mm) [22]
Fig. 8. Measuring setup utilizing the digital image correlation
Table 1. Details of the measuring system
Company
DANTEC Dynamics GmbH, (Ulm, Germany)
Model
Q-400
Cameras
Manta G-507 (4 pieces)
Image resolution
2464 x 2056 pixel
Objective focal distance
35 mm
Field of view
approx. 65 mm x 40 mm
Patterning technique
matt white spray paint base coat _with black speckles_
Pattern feature size
approx. 3 pixel
DIC software
DANTEC Dynamics, Istra 4D (ver. 4.6)
Facet size
19 pixel
Grid spacing
12 pixel
Spatial smoothing
none
Temporal smoothing
Number of acquired data points
20400
The DIC method was carried out with Istra 4D software. Technical details of the measuring system and the adopted DIC settings are summarized in Table 1.
The shape of the measured specimen was determined as presented in section Methods. Two images of the specimen were acquired by changing the background lightening as presented in Figs. 8 and 9. The first image was acquired using the typical lighting setup, presented on the left-hand side of Fig. 10a. For the second image, we illuminated the specimen from the backside, making the background on the acquired image brighter. The acquired image is presented on the right-hand side of Fig. 10a. The shape of the specimen is obtained from both images using Eq. (2).
Fig. 9. Image acquisition using a backlight
Fig. 10. Acquisition of specimen geometry: a) left-hand side: specimen image using typical lighting (Fig. 8) and right-hand side: using the backlight (Fig. 9), b) the obtained specimen geometry
none
Flat Specimen Shape Recognition Based on Full-Field Optical Measurements and Registration Using Mapping Error Minimization Method
328
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 203-213
Results. The determined shape of the measured specimen is presented in Fig. 10b. From the image, one could easily further determine the specimen contour, which most alignment methods use for the alignment procedure. In the proposed method, however, the entire shape is assumed, presented in Fig. 10b in white.
The measured and modelling specimen are together presented in Fig. 11a, given in the measured specimen's coordinate system. The measured shape is presented with the experimentally obtained image and the modelling shape with yellow dots. The corresponding error between both shapes is presented in the form of an error image, Fig. 11b. White pixels correspond to the error set E, obtained using Eq. (8), representing the regions of misalignment.
The initial value of the cost function CF, corresponding to the square number of white pixels in Fig. 11b, is approx. 1012 px2, Eq. (9). The value was then significantly reduced by the optimization procedure in only 5 iterations and the procedure successfully converged in 28 iterations. The final value of the cost function was reduced to approx.
( Model, specimen shape x-direction —^direction
Fig. 11. The initial position of the measured and modelled specimen (the latter shown with yellow dots) in the measuring coord. system, a) presentation of both shapes, b) the corresponding error image of the method
Fig. 12. Alignment of the specimen: a) camera arrangement, b) the initial misalignment and c) the optimal alignment 210	Macek, A. - Urevc, J. - Halilovic, M.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 203-213
106 px2 with the optimal values of the mapping parameters determined to be: a = 2.108°, Ax = 4.907 mm, Ay = 3.091 mm.
The final result of alignment is shown graphically in Fig. 12. The corresponding Fig. 12a schematically presents the experimental set-up. The initial misalignment of the measured and modelling specimen shape is presented in Fig. 12b, shown in the measuring coord. system. The final alignment of specimens is presented in Fig. 12c. Although it is difficult to quantify the level of the alignment, practically no difference between the aligned shapes can be observed.
Handling of the artifacts. To test the ability of the method in handling geometrical mismatches (due to machining or manufacturing tolerances) we purposely misplaced the holes in the measured specimen in the first case, see the upper row in Fig. 13, Case A, where the artifacts are marked with red on the ideal shape shown in white. In the second case, Case B, the roundings on the outer contour were purposely modified to some extent on the measured specimen as if they were not properly machined. The initial position of the modelled specimen is shown with magenta dots. As can be seen in the lower row of the figure, the alignment was performed successfully in both cases.
Fig. 13. Handling of the artifacts; upper row: the assumed artifacts, marked with red on the ideal specimen shape (white), initial modelled specimen shape position, marked with magenta, and lower row: the alignment results; Cases A and B analysed the mismatches in the measured geometry and Case C analysed the presence of noise (outliers)
The ability of the method to handle noise is analysed in Case C, where noise (outliers) are added along the entire contour of the measured specimen. As can be seen in the lower row of the figure, the method managed to successfully perform the alignment. We have also analysed the assumed cases with the ICP method. In Cases A and B, ICP performed the alignment as the proposed method whereas in Case C, the ICP method diverged.
3 DISCUSSION
The purpose of this study was to present a new methodology for aligning full-field measurements. More specifically, the method deals with the alignment of planar specimens, which is a field becoming increasingly popular due to its potential when combining full-field measurements with material characterization. As presented in [22], by using the full-field measurement techniques and monitoring the specimen loads, it is possible to characterize the anisotropic material behaviour from a single heterogeneous strain field specimen by employing inverse identification techniques (eq. finite element model updating, virtual fields method). However, as pointed out in the introduction, the outcome of the material characterization process crucially depends on the accuracy of the alignment between the measured and modelling data.
In the first example, the method was applied to a synthetic experiment. The shape of the specimen resembled a structure of a metamaterial. The example was chosen due to its periodic structure which can cause difficulties in the aligning procedure. We analyse the initial misalignment effect and the performance of the method compared with the ICP algorithm. The results demonstrated the robustness of the proposed method, which managed to successfully perform the alignment even at 90° of the initial rotation. This was not the case for ICP, which managed to perform the alignment up to 45° of the initial rotation only. We need to point out, however, that only the most basic version of the closed-form ICP was employed in this work. There are numerous modifications of the method [23] that could perform better. On the other hand, the performance of the proposed method could also be enhanced, such as by using pixels weights or by modifying the optimization algorithm. By employing such modifications, it is possible to successfully perform the alignment for the entire region of initial misalignments in Fig. 6 (data omitted). However, as is the case with the ICP, the convergence of the proposed
Flat Specimen Shape Recognition Based on Full-Field Optical Measurements and Registration Using Mapping Error Minimization Method
211
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 203-213
alignment algorithm cannot be guaranteed in a general case.
In the majority of cases on which we tested alignment methods, the performance of the ICP and the proposed method was comparable. In general, the advantage of the ICP is its computation efficiency. The method uses singular value decomposition (SVD) for the computation of alignment parameters whereas in the proposed approach an optimization approach is assumed. However, the benefit of the method is its robustness. As demonstrated with the first example, the initial mismatch did not considerably affect the number of iterations whereas they significantly increase in the case of the ICP. This property of the proposed approach comes from using the entire specimen region for the error estimation whereas in the ICP only the contour of a specimen is assumed.
The accuracy of the proposed algorithm on the one hand depends on the quality of the measuring image and on the other on the accuracy of the calibration. The mapping between the observed and modelled geometry is a composite of two mappings. First mapping c1 of images, provided by the cameras, onto the unified plane takes place, subsequently followed by the mapping f of these points onto the space of the modelled specimen. It is important to note that the alignment algorithm addresses only the latter because the former is defined during a calibration procedure of cameras, needed for any full-field measurement system.
In the second example, the specimen recognition and the alignment process were both demonstrated. It can be seen from the acquired specimen's images (Fig. 11) that a simple threshold effect for the contour recognition could not be employed. Due to the speckle pattern, there are similar values of brightness in the background as well as on the specimen. For the presented approach of shape recognition, however, this was not an issue. The alignment of the specimen was in the example performed by using the proposed approach. However, ICP or any other method could be employed as well. From analysing the example (data omitted), similar performance was obtained for both the ICP and the proposed method, with ICP being computationally more efficient. Methods were also analysed for dealing with image artifacts, such as noise and geometrical mismatches due to machining or manufacturing tolerances. In the case of handling the noise, the proposed method turned out to be superior, otherwise, it is known that the ICP algorithm is susceptible to such artifacts [24].
In comparison to alignment procedures published in the field of material characterisation [17] and
[18], the benefit of the proposed approach is that it needs no user input for the shape recognition or for the alignment procedure and it is fairly easy to be implemented. Although we cannot claim that the presented approach for solving the shape recognition and the alignment problem is novel since in both areas there are numerous publications, we found no studies in the field of material characterisation that would address both problems and join them in a form of a methodology.
4 CONCLUSIONS
An integral part of advanced material characterisation by using full-field measurement techniques is the alignment between the FEA data and the experimental DIC data. Despite numerous publications in the field of treating planar specimens, studies that address the subject are scarce.
In the paper, a methodology is presented that enables the alignment of data from a single measuring system. In practice, the alignment approach needs to be addressed once the specimen is placed in the measuring system, more precisely, it needs to be addressed by the measuring system used for the DIC itself. However, because DIC does not recognize the pattern near an object's contour, the location of the contour is not exactly known. With the presented methodology we address both problems, specimen shape recognition, and the alignment procedure. The practical application of the method is presented on two examples, which demonstrate the robustness of the method, its comparison with the ICP algorithm, and its application on real DIC data.
Although the results show that the methodology manages to successfully perform the alignment, it is difficult to actually quantify its accuracy (from both the shape recognition and the alignment process). Such an analysis can be performed by means of computing full-field error maps [10]. For the present methodology, this remains to be performed.
5 ACKNOWLEDGEMENTS
The authors acknowledge the financial support from the Slovenian Research Agency (research core funding No. P2-0263).
6 REFERENCES
[1] Moayyedian, F., Kadkhodayan, M. (2015). Combination of modified Yld2000-2d and Yld2000-2d in anisotropic pressure dependent sheet metals. Latin American Journal of Solids
212
Macek, A. - Urevc, J. - Halilovic, M.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 203-213
and Structures, vol. 12, no. 1, p. 92-114, D0I:10.1590/1679-78251372.
[2]	Starman, B., Vrh, M., Halilovič, M., Štok, B. (2014). Advanced modelling of sheet metal forming considering anisotropy and young's modulus evolution. Strojniški Vestnik - Journal of Mechanical Engineering, vol. 60, no. 2, p. 84-92, D0I:10.5545/sv-jme.2013.1349.
[3]	Grediac, M., Hild, F., Pineau, A. (eds.) (2013). Full-Field Measurements and Identification in Solid Mechanics, John Wiley & Sons, Inc, Hoboken, D0I:10.1002/9781118578469.
[4]	Ruybalid, A.P., Hoefnagels, J.P.M., van der Sluis, O., Geers, M.G.D. (2016). Comparison of the identification performance of conventional FEM updating and integrated DIC. International Journal for Numerical Methods in Engineering, vol. 106, no. 4, p. 298-320, D0I:10.1002/nme.5127.
[5]	Kavanagh, K.T., Clough, R.W. (1971). Finite element applications in the characterization of elastic solids. International Journal of Solids and Structures, vol. 7, no. 1, p. 11-23, D0I:10.1016/0020-7683(71)90015-1.
[6]	Wang, D., Mottershead, J.E. (2016). Measurement precision and spatial resolution with Kriging digital image correlation. Strojniški vestnik - Journal of Mechanical Engineering, vol. 62, no. 7-8, p. 419-429, D0I:10.5545/sv-jme.2016.3736.
[7]	Thengade, A., Dondal, R. (2012). Genetic Algorithm - Survey Paper. IJCA Proceedings on National Conference on Recent Trends in Computing, p. 6.
[8]	Henderson, D., Jacobson, S.H., Johnson, A.W. (2003). The theory and practice of simulated annealing, Glover, F., Kochenberger, G.A. (eds.). Handbook of Metaheuristics (vol. 57), Kluwer Academic Publishers, Boston, p. 287-319, D0I:10.1007/0-306-48056-5_10.
[9]	Fehervary, H., Vastmans, J., Vander Sloten, J., Famaey, N. (2018). How important is sample alignment in planar biaxial testing of anisotropic soft biological tissues? A finite element study. Journal of the Mechanical Behavior of Biomedical Materials, vol. 88, p. 201-216, D0I:10.1016/j. jmbbm.2018.06.024.
[10]	Lava, P., Jones, E.M.C., Wittevrongel, L., Pierron, F. (2020). Validation of finite element models using full field experimental data: Levelling finite element analysis data through a digital image correlation engine. Strain, vol. 56, no. 4, art. ID 12350, D0I:10.1111/str.12350.
[11]	Polyga (2019). 3D Scanning flat objects, from https://www. polyga.com/blog/3d-scanning-flat-objects/, accessed on 2021-01-25.
[12]	Besl, P.J., McKay, D.N. (1992). A method for registration of 3-D shapes. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 14, no. 2, p. 239-256, D0I:10.1109/34.121791.
[13]	Schutz, C., Jost, T., Hugli, H. (1998). Multi-feature matching algorithm for free-form 3D surface registration. Proceedings. Fourteenth International Conference on Pattern Recognition (Cat. No.98EX170), Vol. 2. p. 982-984, D0l:10.1109/ ICPR.1998.711852.
[14]	Rusinkiewicz, S., Levoy, M. (2001). Efficient variants of the ICP algorithm. Proceedings Third International Conference on 3-D Digital Imaging and Modeling, p. 145-152, D0I:10.1109/ IM.2001.924423.
[15]	Bing, J., Vemuri, B. C. (2011). Robust point set registration using gaussian mixture models. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 33, no. 8, p. 16331645, D0I:10.1109/TPAMI.2010.223.
[16]	Eggert, D.W., Lorusso, A., Fisher, R.B. (1997). Estimating 3-D rigid body transformations: a comparison of four major algorithms. Machine Vision and Applications, vol. 9, no. 5-6, p. 272-290, D0I:10.1007/s001380050048.
[17]	Bruno, L., Furgiuele, F.M., Pagnotta, L., Poggialini, A. (2002). A full-field approach for the elastic characterization of anisotropic materials. Optics and Lasers in Engineering, vol. 37, no. 4, p. 417-431, D0I:10.1016/S0143-8166(01)00120-8.
[18]	Silva, G.H.C. (2009). Identification of Material Properties Using Finite Elements and Full-Field Measurements with Focus on the Characterization of Deterministic Experimental Errors. Ecole des Mines de Saint Etienne, Saint Etienne.
[19]	Souto, N., Andrade-Campos, A., Thuillier, S. (2016). A numerical methodology to design heterogeneous mechanical tests. International Journal of Mechanical Sciences, vol. 107, p. 13, D0I:10.1016/j.ijmecsci.2016.01.021.
[20]	Yoneyama, S., Kikuta, H., Kitagawa, A., Kitamura, K. (2006). Lens distortion correction for digital image correlation by measuring rigid body displacement. Optical Engineering, vol. 45, no. 2, p. 9, D0I:10.1117/1.2168411.
[21]	Mirzaali, M.J., Janbez, S., Strano, M., Zadpoor, A.A. (2018). Shape-matching soft mechanical metamaterials. Scientific Reports, vol. 8, art. ID 965, D0I:10.1038/s41598-018-19381-3.
[22]	Macek, A., Starman, B., Mole, N., Halilovic, M. (2020). Calibration of Advanced yield criteria using uniaxial and heterogeneous tensile test data. Metals, vol. 10, no. 4, art. ID 542, D0I:10.3390/met10040542.
[23]	Liu, Y. (2004). Improving ICP with easy implementation for free-form surface matching. Pattern Recognition, vol. 37, no. 2, p. 211-226, D0I:10.1016/S0031-3203(03)00239-5.
[24]	Gao, W., Tedrake, R. (2019). FilterReg: Robust and Efficient Probabilistic Point-Set Registration using Gaussian Filter and Twist Parameterization. 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), p. 1108711096, D0I:10.1109/CVPR.2019.01135.
Flat Specimen Shape Recognition Based on Full-Field Optical Measurements and Registration Using Mapping Error Minimization Method
213
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 214-222 © 2021 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2021.7104	Original Scientific Paper
Received for review: 2021-01-22
Received revised form: 2021-04-04
Accepted for publication: 2021-05-14
Determining Optimum Rotary Blade Design for Wind-Powered Water-Pumping Systems for Local Selected Sites
Abdulbasit Mohammed1 - Hirpa G. Lemu2* - Belete Sirahbizu1
1 Addis Ababa Science and Technology University, Ethiopia 2 University of Stavanger, Norway
The design of a windmill rotor is critical for harnessing wind energy. In this work, a study is conducted to optimize the design and performance of a rotor blade that is suitable for low wind conditions. The windmills' rotor blades are aerodynamically designed based on the SG6043 airfoil and wind speed data at local selected sites. The aerodynamic profile of the rotor blade that can provide a maximum power coefficient, which is the relation between real rotor performance and the available wind energy on a given reference area, was calculated. Different parameters, such as blade shapes, chord distributions, tip speed ratio, geometries set angles, etc., were used to optimize the blade design with the objective of extracting maximum wind power for a water pumping system. Windmill rotor of 10.74 m, 7.34 m, and 6.34 m diameter with three blades were obtained for the selected sites at Abomsa, Metehara, and Ziway in south-east Ethiopia. During the rotary blades performance optimization, blade element momentum (BEM) theory and solving iteration by MATLAB® coding were used. Keywords: boundary element method, rotary sizing, wind power, hydraulic power, power coefficient, water pumping system
Highlights
•	Optimum rotary blade for wind power at low wind speed has been designed based on SG603 airfoil and analysed for three selected local sites.
•	The blade design is optimized with the objective of extracting maximum wind power for each selected site for a water-pumping system.
•	The relation between real rotor performance and the available wind energy on the selected given reference areas is calculated to find the optimum aerodynamic profile of the blades.
•	Blade element momentum theory has been used and iterated in MATLAB to identify the parameters for the performance optimization.
0 INTRODUCTION
Currently, significant global progress to benefit from the need for renewable energy by converting the energy of the wind to a mechanical form by designing and constructing diverse forms of windmills exists [1]. The main interest in wind energy is that the conversion process does not release carbon emissions to the environment, and it leads to less consumption of resources. In order to mitigate the side effects of fossil fuels, the use of renewable energy sources has become extremely important. In this respect, wind energy is considered one of the most promising renewable energy sources, which forces of nature can continuously generate, and the rotor blade design is decisive in optimally harvesting wind energy.
Studies on wind turbine rotor blade design have been continuously conducted based on the performance of aerodynamic analysis. Supreeth et al. [2] used a heuristic approach to study rotary blade design for a fixed pitch and small-scale horizontal axis wind turbine. The sectional chord and twist angle distribution for the idealized, optimized, and linearized blades were analytically determined. Prasad and Virupaxi [3] designed the rotary blade of windmills
for irrigational purposes in Bagur of Hosadurga district (Karnataka, India). Various parameters were optimized to achieve maximum power coefficient. Lopez-Lopez et al. [4] studied the effect of the blade slenderness and the wind speed on the dynamics and instability of wind turbine blade under large deflection. They used a simple cantilever beam model and employed a Galerkin approach validated with experimental results. Lee and Shin [5] also improved an optimization framework for blade aerodynamic design under realistic conditions.
The design parameters, such as angle of twist, chord distribution, maximum chord and the objective function relation, were obtained. Shubham et al. [6] designed and optimized the micro wind turbine blades at a rated wind speed of 8.4 m/s using QBlade, with which the SG6043 airfoil was selected. and simulation was done on a 1.2 m long blade. Chan [7] reported the design of a 1.2 kW wind turbine blade for rural application for Mandalay Hill of Mandalay City (Myanmar) at a rated wind velocity of 7 m/s. The blade airfoil, angle of pitches and distribution of sectional chord along the blade span were determined. Deepak et al. [8] designed and optimized the blade of a small wind turbine by considering various factors,
214
*Corr. Author's Address: University of Stavanger, N-4036, Stavanger, Norway, Hirpa.g.lemu@uis.no
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 214-222
such as tip loss, hub loss, drag coefficient, and wake. By using blade element momentum (BEM) methods, power performance was simulated. Muhsen et al. [9] optimized the rotor blade parameters that affect blade design for small turbines with low performance due to low Reynolds number. In the optimization processes, the performance of the final designed blade performance was investigated by employing the blade element momentum theory. Norona and Krishna [10] designed a 300 W micro-horizontal axis wind turbine using BEM theory on SG6040 and SG6041 airfoils. The BEM procedure was codified in MATLAB software, and a 3 m diameter wind turbine was analysed with QBlade. In general, the literature study indicates that the efficient conversion of the wind energy directly to a mechanical form requires closer study of the local wind data and designing rotor blades that are tailormade for the wind condition.
This study focused on designing optimum rotary blades for a wind-powered water-pumping system for local selected sites at Abomsa, Metehara and Ziway (Ethiopia). This design of the blades is based on the design of an aerodynamic type of SG6043 airfoil reported by Giguere and Selig [11] and Osei et al. [12] and analysis of collected wind speed data of the sites [13] from the National Meteorological Agency (NMA). The blade parameters, including chord distribution length, angle of twist, tip speed ratio, power coefficient of the blade, and attack angle, were determined. These blades' parameters were optimized to attain maximum power coefficients by using BEM theory procedures, and different iterations were solved by MATLAB coding.
1 DESIGN METHODS AND PROCEDURES
The general procedures followed to design the rotor blades of the windmills for water pumping application are summarized in Fig. 1.
Data		Rotor
collection		siting
Blade design
No. of blades ■ Blade profile - Blade geometry
Power computation
Fig. 1. Design procedure of rotary blade of wind turbines
The weather data were collected from selected local cites: Abomsa, Metehara and Ziyay, all from the eastern part of the Oromia regional state (Ethiopia). The next step is developing new designs and configurations of the rotor blades (rotor blade sizing) that can efficiently transform the energy in the wind to direct mechanical energy to drive water pumps. Fig. 2 describes the design steps that must be taken
to determine the optimum rotary size of windmills to be used at those particular sites and for a particular purpose.
Fig. 2. Illustration of the design steps for sizing of windmill rotor for water-pumping system
Daily water requirement (QD) depends on the total number of beneficiaries (NB) of the site and the daily water consumptions per capita (qpp), expressed as Eq. (1).
Qd = Nb • qp
(1)
Hydraulic power requirement (Phyd) parameter is a function of the daily water requirement (QD) and total pumping head (H) as given by Eq. (2) [14] to [16], where the total pumping head includes the static water head (HSt) and friction losses (Hf), considering head losses in the pipe, velocity head losses, and other minor losses.
P,d = 0.113452 • Qd • H,
(2)
where H = HSt + Hfr. The availability of wind power potential on a monthly basis is essential to evaluate as given in Eq. (3) [14] to [18].
P = 10 . v
wind 2 üir
(3)
Reference area (Ar) is related to the rotor area needed to capture sufficient power from the wind as given in Eq. (4) [14] to [18].
P „
(4)
Ar =
hyd
P
Determining Optimum Rotary Blade Design for Wind-Powered Water-Pumping Systems for Local Selected Sites
215
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 214-222
1.1 Blade Design Procedures
1.2 Rotor Blade Performance
After the rotary blade size of horizontal axis windmills
was determined, the following steps were considered
in rotary blade design.
a.	Determining the minimum Cd/ C¡ ratio: This is important to determine the value of design lift and angle of attack corresponding to a minimum value of Cd/ C¡ ratio.
b.	Determining design lift coefficient and design angle of attack (ad): From airfoil data corresponding to minimum Cd/ C¡ ratio, the value of lift coefficient and angle of attack are found.
c.	Choosing design tip speed ratio (A) and the number of blades (B): These parameters can be chosen from tip speed ratio (TSR) vs the number of blade table depending on the form of operation.
d.	Determining geometries of optimum blades: When rotor blade size is determined, then it is divided into 10 to 20 equal segments or elements; for each unit element, the blade geometries such as local tip speed ratio A(rt), relative wind angle p(r,), chord distribution c(rl), pitch angle 6(r), twist angle ^(r¿), and local solidity a(r) can be computed as shown in Fig. 3.
e.	Linearization of the twist angle and blade chord length: Using the optimum blade shapes as a guide, a blade shape that promises to be a good approximation is selected. Since the chords as well as the blade twist vary in a non-linear manner along the blade, it is possible to linearize the chords and the twist angles using Eq. (5).
r + a~
- (r ) = a
) = a3 (R - r ),
(5)
where ah a2 and a3 are constant coefficients obtained from the chosen chord distributions and twist angles.
Fig. 3. Computation procedures to determine blade geometries
After the chords and twist angles are linearized, the geometry of the rotary blade has a practical shape. The performance study of the designed rotor blade is conducted by using BEM theory [19]. Obtaining axial a(r,)j, and angular induction factors d(ri)j of the designed blade is the base of rotor blade performances. It is possible to determine the inflow angle using the axial and angular induction factors, as follows (Eq. (6)).
) =tan-1
1 - « (r )
i
(i+á (r )J )(r)
(6)
To describe the lifting power reduction by airflow around the blade tip, the tip loss correction factor F(r) for each /'th section is calculated based on Prandtl's method [20]:

2
— I cos n
wir
(7)
The axial and tangential induction factors for each sectional /'th element at iteration jth are then calculated using Eq. (8). If the axial induction factor is greater than 0.4, the Froude momentum principle is no longer accurate. Thus, Bangga [21] suggested correction of the axial induction factor, if a(r) > ac, using Eq. (9).
i	v1
4F(ri )sin2^ (ri)
a (ri )j = á (ri )j =
a„
i (Ci c°w(r )+cd sin^ (r ))
4F (r )cosy (r )siny (r ) ■ (r )(cisiny(r ) - cd i may (r ))
+1
V1
-1
> (8)
l(r¡ ) = - {2 + Z(r )-
}/(z(rl )(1 - 2ac ) + 2)2 + 4(z(r )a- -l) j, (9)
where ac is commonly about 0.2 and Z(r) is defined as:
* (r )=■
4F (r¡ )sin2 y/(r¡ )
(r ){ci ,i cos¥ (r )+Cd ,i sw (r ))
a
(10)
Finally, the local thrust coefficient CT(r) and power coefficient CP(r) can be computed from the following relations.
B
2
R
e
216
Mohammed, A. - Lemu, H.G. - Sirahbizu, B.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 214-222
C (r ) = a(rt )(1 - a (rt )
(Ci jcow(r ) + Cd,i )) F (r )smV(r ) :
Cp (r )— F (r )a3 (r )â (r )J
(1 - « (r )j
C
1 -
C
cot^(r ).
(ii)
2 ANALYSIS OF RESULTS
The rotary blades have been designed according to data collected by NMA of Ethiopia from three selected sites (Abomsa, Metehara, and Ziway, Ethiopia). The data was collected for one year (October 2018 to September 2019) at 10 m height. Table 1 shows important data determined from the three selected
sites, and total water demand and pumping required can be calculated by estimating the daily water consumption per capita to be 20 l/person in Ethiopia [14] and [15]. The monthly mean wind speeds were extrapolated at 20 m heights (hub heights), and monthly specific wind power potentials for the three selected sites were evaluated (Table 2).
Table 1. Important parameters for three selected sites
Parameters	Abomsa	Metehara	Ziway
No. of beneficiaries [-]	4086	4191	3867
Total demand [l]	81720	83820	77340
Pumping required [m3/day]	81.72	83.82	77.34
Total head [m]	87	79	74
Tip speed ratio, (X) [-]	3	3.5	4
Number of blades, (B) [-]	3	3	3
Based on the number of beneficiaries and total water demands, hydraulic power can be determined,
Table 2. Available mean wind speeds and specific wind power for three selected sites
			Abomsa			Metehara			Ziway	
Year	Month	Vm at 10 [m/s]	Vm at 20 [m/s]	p 1 wind [W]	Vm at 10 [m/s]	Vm at 20 [m/s]	p wind [W]	Vm at 10 [m/s]	Vm at 20 [m/s]	p wind [W]
	Oct	2.82	3.42	24.59	2.93	3.56	27.58	3.82	4.64	61.12
2018	Nov	2.81	3.41	24.33	2.72	3.30	22.06	3.92	4.76	66.04
	Dec	2.99	3.63	29.31	2.63	3.19	19.95	3.56	4.32	49.47
	Jan	3.82	4.64	61.12	2.63	3.19	19.95	3.67	4.46	54.20
	Feb	3.63	4.41	52.44	2.92	3.55	27.30	3.05	3.70	31.11
	Mar	3.27	3.97	38.34	2.53	3.07	17.76	2.91	3.53	27.02
	Apr	2.33	2.83	13.87	2.65	3.22	20.40	2.66	3.23	20.64
2019	May	2.76	3.35	23.05	2.75	3.34	22.80	3.01	3.65	29.90
	Jun	2.58	3.13	18.83	3.39	4.12	42.71	4.14	5.03	77.80
	July	2.42	2.94	15.54	3.36	4.08	41.59	3.20	3.89	35.93
	Aug	2.01	2.44	8.90	3.23	3.92	36.95	3.64	4.42	52.88
	Sept	2.16	2.62	11.05	2.75	3.34	22.80	2.72	3.30	22.06
Table 3. Hydraulic power, reference area and rotor diameter for three selected sites										
Year			Abomsa			Metehara			Ziway	
	Month "	Phyd	Ar	Dr	Phyd	Ar	Dr	Phyd	Ar	Dr
	Oct	806.3	32.8	6.5	750.9	27.2	5.9	649.1	10.6	3.7
2018	Nov	806.3	33.1	6.5	750.9	34.0	6.6	649.1	9.8	3.5
	Dec	806.3	27.5	5.9	750.9	37.7	6.9	649.1	13.1	4.1
	Jan	806.3	13.2	4.1	750.9	37.7	6.9	649.1	12.0	3.9
	Feb	806.3	15.4	4.4	750.9	27.5	5.9	649.1	20.9	5.2
	Mar	806.3	21.0	5.2	750.9	42.3	7.3	649.1	24.0	5.5
	Apr	806.3	58.1	8.6	750.9	36.8	6.9	649.1	31.5	6.3
2019	May	806.3	35.0	6.7	750.9	32.9	6.5	649.1	21.7	5.3
	Jun	806.3	42.8	7.4	750.9	17.6	4.7	649.1	8.3	3.3
	July	806.3	51.9	8.1	750.9	18.1	4.8	649.1	18.1	4.8
	Aug	806.3	90.6	10.7	750.9	20.3	5.1	649.1	12.3	4.0
	Sept	806.3	73.0	9.6	750.9	32.9	6.5	649.1	29.4	6.1
Determining Optimum Rotary Blade Design for Wind-Powered Water-Pumping Systems for Local Selected Sites	217
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 214-222
a)
Abomsa Re — 0.23 Metehara Re = 0.25 Ziway üe- 0.29
10 20 Attack angle, a [deg]
Abomsa Re = 0.23 Mclehara Re = 0.25 Ziway Re = 0.29
0	10 20 30
kj	Attack angle, a [dcg]
Fig. 4. Lift coefficient a) C, and b) C,/Cd ofSG6043 with Reynolds numbers of three sites
Fig. 5. Flow diagram of performing rotor blades design
which is constant through the year, the reference area and rotor diameter for the three sites were calculated and summarized in Table 3. The rotor size is chosen on the basis of "design month", which is the month in which the water demand is highest in relation to the wind power resources, i.e., the month when the system will be most heavily loaded [14] and [16].
Therefore, according to the data given in Table 3, August, March and April are the design months for Abomsa, Metehara and Ziway, respectively.
For three rotor blade design, Giguere and Selig [11] and Osei [12] recommended airfoil type SG6043 for low Reynolds numbers because this type of airfoil operates in regional low wind speed. The designed C{
Table 4. Geometric properties of optimum and linearized rotor blades for Abomsa site
rt/R		Optimal blade				Linearized blade		
	W(rt)	m	¿(r)	c(rt)	W(rt)	m	¿(r)	c(rd
0.07	52.46	47.46	40.17	0.39	24.83	19.83	20.79	0.77
0.13	45.47	40.47	33.18	0.59	23.34	18.34	19.30	0.73
0.20	39.36	34.36	27.07	0.68	21.86	16.86	17.82	0.70
0.27	34.23	29.23	21.94	0.69	20.37	15.37	16.33	0.67
0.33	30.00	25.00	17.71	0.66	18.89	13.89	14.85	0.64
0.40	26.54	21.54	14.25	0.63	17.40	12.40	13.36	0.61
0.47	23.69	18.69	11.40	0.59	15.92	10.92	11.88	0.58
0.53	21.34	16.34	9.05	0.54	14.43	9.43	10.39	0.54
0.60	19.37	14.37	7.08	0.51	12.95	7.95	8.91	0.51
0.67	17.71	12.71	5.42	0.47	11.46	6.46	7.42	0.48
0.73	16.30	11.30	4.01	0.44	9.98	4.98	5.94	0.45
0.80	15.08	10.08	2.79	0.41	8.49	3.49	4.45	0.42
0.87	14.03	9.03	1.74	0.38	7.01	2.01	2.97	0.39
0.93	13.10	8.10	0.81	0.36	5.52	0.52	1.48	0.35
1.00	12.29	7.29	0.00	0.34	4.04	-0.96	0.00	0.32
218	Mohammed, A. - Lemu, H.G. - Sirahbizu, B.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 214-222
a)
c)
b)
1.0 2.0 3.0 Rotor radius, R [m]
Fig. 6. Local chord and twist angle comparison of optimal and linearized blades for three sites, a) Abomsa, b) Metehara, and c) Ziway
Table 5. Geometric properties of optimum and linearized rotor blade for Metehara site
r / R		Optimal blade				Linearized blade		
	W(rt)	m	¿(r)	c(rd	W(rt)	m	¿(r)	c(rd
0.07	51.24	46.24	40.61	0.43	23.04	18.04	18.98	0.73
0.13	43.32	38.32	32.69	0.63	21.68	16.68	17.62	0.70
0.20	36.67	31.67	26.04	0.69	20.33	15.33	16.27	0.66
0.27	31.32	26.32	20.69	0.67	18.97	13.97	14.91	0.63
0.33	27.07	22.07	16.44	0.63	17.62	12.62	13.56	0.60
0.40	23.69	18.69	13.06	0.59	16.26	11.26	12.20	0.57
0.47	20.98	15.98	10.35	0.54	14.91	9.91	10.85	0.54
0.53	18.79	13.79	8.16	0.49	13.55	8.55	9.49	0.50
0.60	16.98	11.98	6.35	0.45	12.19	7.19	8.13	0.47
0.67	15.47	10.47	4.84	0.42	10.84	5.84	6.78	0.44
0.73	14.19	9.19	3.56	0.39	9.48	4.48	5.42	0.41
0.80	13.10	8.10	2.47	0.36	8.13	3.13	4.07	0.38
0.87	12.16	7.16	1.53	0.34	6.77	1.77	2.71	0.35
0.93	11.35	6.35	0.72	0.32	5.42	0.42	1.36	0.31
1.00	10.63	5.63	0.00	0.30	4.06	-0.94	0.00	0.28
Determining Optimum Rotary Blade Design for Wind-Powered Water-Pumping Systems for Local Selected Sites	219
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 214-222
0.50
0.40 --
0.30 --
020 --
0.10 --
0.00
						
						
						
				Metehara		
				Ziway		
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 Tip speed ratick }Ji\} [deg]
Fig. 7. The vs. of rotor performed blades for three selected sites
and Ci/Cd graphs characteristics of SG6043 with Re of 0.23x106, 0.25x106 and 0.29x106 are evaluated, as illustrated in Fig. 4, using QBlade, for Abomsa, Metehara, and Ziway, respectively. Fig. 5 shows the general flow diagram for determining blade geometries and performing the designed blades by iterating axial and tangential induction factors.
The blade geometries are determined with the number of segments or elements (V) equal to 15.
The local chord distributions, geometries set angles, and solidities of optimal and linearized blades for the three selected sites are determined and summarized in Tables 4, 5 and 6. The comparison of blade chord distributions and sectional angles of twist for optimal and linearized blades for the three selected sites are illustrated in Fig. 6.
Based on the BEM flow diagram of iteration in Fig. 5, each local axial and tangential induction factors have converged after many iterations by MATALB coding and the coefficients of performance are illustrated as shown in Fig. 7 for these three local sites. The geometry profiles (shape) of designed and optimized rotor blades are modelled by QBlade software, and the resulting blades are shown in Fig. 8. As can be observed from this figure, the maximum coefficient of performance values of CP(ri) are 0.448, 0.459 and 0.463 at the maximum tip speed ratio of A(ri) 3.0, 3.5 and 4.0 for the three sites Abomsa, Metehara and Ziway respectively. The geometry shape of optimal, linearized and optimized designed blades modelled by QBlade software and the resulting blades are shown in Fig. 8.
Fig. 8. Geometry shape of optimal, linearized and optimized blades modelled by QBlade for A) Abomsa, B) Metehara and C) Ziway sites.
Table 6. Geometric properties of Optimum and linearized rotor blades for Ziway site
Optimal blade
Linearized blade
ri i r	¥(rd	m	m	c(rd	¥(ri)	m	m	C(rd
0.07	50.05	45.05	40.69	0.47	21.41	16.41	17.29	0.67
0.13	41.29	36.29	31.93	0.66	20.17	15.17	16.05	0.64
0.20	34.23	29.23	24.87	0.69	18.94	13.94	14.82	0.61
0.27	28.77	23.77	19.41	0.65	17.70	12.70	13.58	0.58
0.33	24.58	19.58	15.22	0.60	16.47	11.47	12.35	0.55
0.40	21.34	16.34	11.98	0.54	15.23	10.23	11.11	0.52
0.47	18.79	13.79	9.43	0.49	14.00	9.00	9.88	0.49
0.53	16.74	11.74	7.38	0.45	12.76	7.76	8.64	0.46
0.60	15.08	10.08	5.72	0.41	11.53	6.53	7.41	0.43
0.67	13.70	8.70	4.34	0.38	10.30	5.30	6.18	0.39
0.73	12.55	7.55	3.19	0.35	9.06	4.06	4.94	0.36
0.80	11.57	6.57	2.21	0.32	7.83	2.83	3.71	0.33
0.87	10.73	5.73	1.37	0.30	6.59	1.59	2.47	0.30
0.93	10.00	5.00	0.64	0.28	5.36	0.36	1.24	0.27
1.00	9.36	4.36	0.00	0.26	4.12	-0.88	0.00	0.24
220
Mohammed, A. - Lemu, H.G. - Sirahbizu, B.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 214-222
3 CONCLUSION
In this study, the optimum rotary blade for a wind-powered water pumping system for local three selected sites has been designed and analysed. From these designs, the optimal blade parameters are determined, linearized, and optimized. The geometry profiles of optimized blades are modelled using QBlade software for the three local selected sites. The key results of the design and analysis work show that the maximum coefficient of performance values Cp(r) are 0.4512, 0.4587, and 0.4627 at the maximum tip speed ratio values of 3.0, 3.5 and 4.0 for the three sites Abomsa, Metehara, and Ziway, respectively. Furthermore, the average wind speeds at a height of 20 m are 3.50 m/s, 3.60 m/s and 4.60 m/s, respectively.
This article is part of a PhD research whose further work will focus on developing the design model of the blade with the obtained design specification and building a prototype wind turbine to test at the three locations.
4 NOMENCLATURE
	Axial induction factor, [-]
à(r,)	Angular induction factor, [-]
	Constant coefficient, [-]
Ar '	Reference area, [m2]
c(r i)	Chord distribution, [-]
Cd,i	Local drag coefficient, [-]
Cd	Drag coefficient, [-]
Ci,i	Local lift coefficient, [-]
C	Lift coefficient, [-]
Cp(ri)	Local power coefficient, [-]
CAr,)	Local thrust coefficient, [-]
F(r i)	Blade tip loss factor, [-]
H	Total pumping head, [m]
Hfr	Friction losses, [m]
Hst	Static water head, [m]
Nb	Total number of beneficiaries [-]
Phyd	Hydraulic power requirement [W]
p 1 wind	Wind power potential [W/m2]
Qd	Daily water requirement, [m3/day]
qpp	Water consumption per capita, [l/day]
R	Rotor blade radius, [m]
ri	Local radius, [m]
V * m	Wind velocity, [m/s]
ad	Design attack angle, [deg]
m	Pitch angle, [deg]
xW	Local tip speed ratio, [-]
X	Tip speed ratio, [-]
P car	Air density, [kg/m3]
ofc)	Local solidity, [m-1]
$(rt) Twist angle, [deg]
^(r.) Relative wind angle, [deg]
5 REFERENCES
[1]	Ahamed, E.Z., Gowtham, A.R. (2015). Aero-design analysis for modified Darrieus based-straight bladed VAWT systems. International Journal of Recent Technology and Engineering, vol. 4, no. 5, p. 510-514, DOI:10.15623/ijret.2015.0405095.
[2]	Supreeth, R., Arokkiaswamy, A, Nagarjun, J.R., Prajwal, H.P, Sudhanva, M. (2019). Geometrical design of a rotor blade for a small-scale horizontal axis wind turbine. International Journal of Recent Technology and Engineering, vol. 8, no. 3, p. 33903400.
[3]	Prasad, S.S., Virupaxi, A. (2012). Optimized design of rotor blade for a wind pump. International Journal of Renewable Energy Research, vol. 2, no. 4, P. 746-749.
[4]	Lopez-Lopez, A., Robles-Ocampo, J.B., Sevilla-Camacho, P.Y., Lastres-Danguillecourt, O., Muniz, J., Perez-Sarinana, B.Y., de la Cruz, S. (2020). Dynamic instability of a wind turbine blade due to large deflections: An experimental validation. Strojniški vestnik - Journal of Mechanical Engineering, vol. 66, no. 9, p. 523- 533, DOI:10.5545/sv-jme.2020.6678.
[5]	Lee, S.-L., Shin, S.J. (2020). Wind turbine blade optimal design considering multi-parameters and response surface method. Energies, vol. 13, no. 7, art. ID 1639, DOI:10.3390/ en13071639.
[6]	Shubham, R., Shubham, S., Rohan, S., Navjyot, S., Chaudhary, M.K. (2017). Simulation of micro wind turbine blade in QBlade. International Journal for Research in Applied Science and Engineering Technology, vol. 5, no. IV, p. 256-262.
[7]	Chan, Z.M. (2019). Design calculation of 1200W horizontal axis wind turbine blade for rural applications. International Journal of Trend in Scientific Research and Development, vol. 3, no. 5, p. 500-507.
[8]	Deepak, J.N., Chandan, R., Doddanna, K. (2017). Design & structural analysis of a small wind turbine blade for operation at low wind speed. International Research Journal of Engineering and Technology, vol. 4, no. 11, p. 2188-2191.
[9]	Muhsen, H. Al-Kouz, W., Khan, W. (2020). Small wind turbine blade design and optimization. Symmetry, vol. 12, no. 1, art. ID 18, DOI:10.3390/sym12010018.
[10]	Noronha, N.P., Krishna, M. (2020). Design and analysis of micro horizontal axis wind turbine using MATLAB and QBlade. International Journal of Advanced Science and Technology, vol. 9, no. 10S, p. 8877-8885.
[11]	Giguer'e, P., Selig, M.S. (1998). New airfoils for small horizontal, axis wind turbines. Journal of Solar Energy Engineering, vol. 120, no. 2, p. 108-114, D0I:10.1115/1.2888052.
[12]	Osei, E.Y., Opoku, R., Sunnu, A.K., Adaramola, M.S. (2020). Development of high performance airfoils for application in small wind turbine power generation. Journal of Energy, vol 2020, art ID. 9710189, DOI:10.1155/2020/9710189.
[13]	Mohammed, A., Lemu, H.G., Sirahbizu, B. (2020). Statistical analysis of Ethiopian wind power potential at selected sites. Proceedings of 8th EAI International Conference on Advancements of Science and Technology, Bahir Dar.
Determining Optimum Rotary Blade Design for Wind-Powered Water-Pumping Systems for Local Selected Sites
221
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 214-222
[14]	Ronak, D.G., Pramod, K., Debarshi, S., Bhushan, K., Shubham, C. (2015). Design and development of windmill operated water pump. International Journal of Innovative Research in Science, Engineering and Technology, vol. 5, no. 1, p. 61-74.
[15]	Girma, M., Molina, M., Assefa, A. (2015). Feasibility study of a wind powered water pumping system for rural Ethiopia. AIMS Energy, vol. 3, no. 4, p. 851-868, DOI:10.3934/ energy.2015.4.851.
[16]	African Water Facility (AWF) (2008). Ethiopia-Utilization of solar and wind energy for rural water supply In Ethiopia. Appraisal report, Ministry of Water Resources, Addis Abeba.
[17]	Sathyajith, M. (2006). Wind Energy, Fundamentals, Resource Analysis and Economics. Springer Verlag, Berlin Heidelberg, DOI:10.1007/3-540-30906-3.
[18]	Ahmmed, S. (2002). Investigation and Analysis of Wind Pumping System for Irrigation in Bangladesh. Bangladesh University of Engineering and Technology, PhD Thesis, Dhaka.
[19]	Yang, K. (2020). Geometry design optimization of a wind turbine blade considering effects on aerodynamic performance by linearization. Energies, vol. 13, no. 9, art. ID 2320, DOI:10.3390/en13092320.
[20]	Burton, T., Sharpe, D. (2006). Wind Energy Handbook. John Wiley & Sons Ltd., Chichester.
[21]	Bangga, G. (2018). Comparison of blade element method and CFD simulations of a 10 MW Wind Turbine. Fluids, vol. 3, no. 4, art ID 73, DOI:10.3390/fluids3040073.
222
Mohammed, A. - Lemu, H.G. - Sirahbizu, B.
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 223-234 © 2021 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2021.7089	Original Scientific Paper
Received for review: 2021-01-05
Received revised form: 2021-04-20
Accepted for publication: 2021-05-03
Investigation on the Modeling and Dynamic Characteristics of a Novel Hydraulic Proportional Valve Driven by a Voice Coil Motor
Mingxing Han1* - Yinshui Liu2 - Yitao Liao1 - Shucai Wang1 iHuazhong Agricultural University, College of Engineering, China 2Huazhong University of Science and Technology, School of Mechanical Science and Engineering, China
As the key control component of the water hydraulic systems, the water hydraulic proportional valve has a significant influence on the control performance of the systems. Due to the poor viscosity and lubricity of water, the valve spool resistance is large and non-linear. In this study a novel fast-response water hydraulic proportional valve is presented. The actuator of the valve adopts a voice coil motor (VCM), which has the advantages of fast response, high control precision and small volume. In order to realize the fast control of the valve, a lever amplifier is designed to obtain enough actuation force. A detailed and precise non-linear mathematical model of the valve considering both the valve's structural parameters and VCM electromagnetic characteristics is developed. A comprehensive performance simulation analysis has been carried out, mainly divided into an electromagnetic simulation, an analysis of the characteristics of the lever magnifier, and a dynamic performance simulation of the valve. The simulation results show that the adjusting time is about 28 ms, and the maximum overshoot is about 5 %. The step response rise time is about 15 ms. The test rig of the valve and VCM have been built. The test results of the prototype show that the optimal stroke range of VCM is 4 mm to 15 mm. The maximum overshoot of the valve is around 10 %; the adjusting time is about 30 ms in the opening process and 35 ms in the closing process. The test results prove that the valve has good static and dynamic control performance. Keywords: proportional control valve, voice coil motor, optimization design, dynamic performance
Highlights
•	A novel water hydraulic proportional valve with fast response is proposed.
•	The actuator of the valve adopts a voice coil motor (VCM), and a lever amplifier is designed to obtain enough actuation force to realize the fast and accurate control.
•	The detailed and precise mathematical models of the valve have been developed, and the comprehensive optimization design has been carried out.
•	The test rig has been built, and the results indicate that the step response time of 0 % to 100 % stroke is about 30 ms.
0 INTRODUCTION
Compared with oil hydraulic systems, water hydraulic systems could be a good solution for the environmental and safety problems [1]. Water hydraulic systems are widely used in steel and glass production, ocean exploration, food and medicine processing, and coal mining [2] to [4]. The water hydraulic control valve is one of the most important components of hydraulic systems, and its dynamic performance has an important impact on the performance of hydraulic systems [5]. However, it is quite difficult to design a high-performance water hydraulic proportional valve due to the low viscosity and oxidative corrosion of water. Many problems, such as poor lubrication and large friction, sealing, and leakage, need to be solved in the design process. Some studies [6] have shown that while the oil directional control sliding valve worked perfectly, there is a strengthening oscillatory movement of the same water valve spool, caused by non-optimal valve material combination (steel/steel) in water-lubricated conditions.
The current research indicates that the dynamic performance of electro-mechanical devices is one of the key factors affecting the response speed of the control valve [7]. Considering that the high-speed on/ off solenoids have fast response, some water hydraulic control valves are designed with them as the actuators [8] and [9]. However, they are usually used in highspeed switching valves with small flow capacity due to the small actuation force and stroke. Park et al. [10] and Park [11] proposed a water pressure proportional valve with the high-speed switching valve as a pilot stage, and the experiments show that the maximum flow rate can reach 9 l/min and the pressure drop is 7 MPa. This kind of valve can realize proportional control by pulse width modulation (PWM) signal but cannot achieve a good control effect. Piezoelectric actuators have become a solution for improving the valve dynamic performance because of their high response speed and high output force. However, their stroke usually is very small even at a large applied voltage, and the properties are greatly affected by temperature [12] and [13]. These drawbacks restrict
*Corr. Author's Address: College of Engineering, Huazhong Agricultural University, Wuhan, China, 430070. hanmingxing@mail.hzau.edu.cn
223
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 223-234
the valve opening (which in turn limits maximum flow rate), and it is difficult to adapt to the working environment. Compared with the high-speed on/off solenoids and piezoelectric actuators, a voice coil motor (VCM) has the obvious advantages of fast response and large stroke-actuation force. Therefore, as a high-performance actuator, a voice coil motor has attracted the attention of researchers. Dahoon et al. [14] designed a high force voice coil motor. VCM usually can provide 0.7 N to 1000 N actuation force and 0 mm to 100 mm stroke. The motion frequency of VCM can reach 1000 Hz. It is widely used in medical, aviation, vibration platform control, etc. A comparison of these existing electro-mechanical devices shows that VCM is an appropriate choice for electromechanical actuators in water hydraulic control valve.
As an ideal electric-mechanical conversion device, VCM can be used to quickly and accurately control the movement of the valve spool. Li et al. [15] and [16] proposed a pneumatic servo-valve with control pressure up to 20 MPa, which was driven by VCM. Experiments show that the valve has good control performance, and the step response frequency can reach 100 Hz. Zhang et al. [17] proposed a direct-drive water hydraulic control valve driven by a voice coil motor, and the rated flow is 14 l/min. The test results show that the control valve has the performance of approximate linear proportional control. In contrast, other researchers are committed to improving the dynamic performance of the valve by optimizing structural parameters. Liu et al. [18] studied a novel throttle poppet valve based on the hydraulic feedback principle. Filo et al. [19] proposed a novel control valve that consists of a throttle valve and a differential valve controlled by the pressure difference between the supply line and the return line to increase the speed of cylinder piston rod movement. Xu et al. [20] studied the mathematical modelling and simulation of a novel hydraulic variable valve time system. Simic and Herakovic [21] pointed out that non-optimized
hydraulic valve geometry is the main cause for many problems related to response time, actuation force and energy consumption. The results of the simulation analyses show that the axial component of the flow forces could be reduced significantly by optimizing the geometry of the valve spool and housing.
Due to the poor viscosity and lubricity of the water medium, the conventional hydraulic valve structure and the optimization design method will be not suitable for the water hydraulic proportional valve. In this study, a novel fast-response water hydraulic proportional valve that driven by VCM is presented. A comprehensive performance simulation analysis has been carried out, which mainly includes electromagnetic simulation of VCM, characteristics analysis of lever magnifier and dynamic performance simulation of the valve.
1 STRUCTURE AND PRINCIPLE OF THE VALVE
The water hydraulic valve proposed in this study uses a lever amplifier to amplify the VCM actuation force (as shown in Fig. 1). A smaller VCM can be used to drive the spool to move quickly. The traditional oil hydraulic control valve mostly adopts the slide valve structure. Due to its small sealing clearance, the working fluid cleanliness level has a significant impact on the valve performance [22]. The existing research shows that water hydraulic control valve with non-gap sealing structures, such as the seat valve (ball valve, disc valve, poppet valve etc.), can effectively solve the leakage problem caused by the low viscosity of water and improve the anti-pollution ability of the valve. The valve orifice adopts a seat valve structure, which is mainly composed of seat and ball. The displacement sensor detects the spool position, and the displacement signal will be fed back to the PID controller for closed-loop control.
The water hydraulic control valve with non-gap sealing structures such as seat valve (ball valve, disc
224
Fig. 1. Structure of the water hydraulic proportional valve
Han, M. - Liu, Y. - Liao, Y. - Wang, S.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 223-234
Fig. 2. The valve control schematic
valve, poppet valve etc.) can effectively solve the leakage problem caused by the low viscosity of water and high working pressure [7] and [11]. The valve orifice adopts a seat valve structure, which is mainly composed of valve seat and ball. The displacement sensor detects the spool position, and the displacement signal will be fed back to the controller for closed-loop control. The main parameters of the prototype are shown in Table 1.
Table 1. The main parameters of the valve
Parameter	Symbol	Value	unit
Fluid density	P	1000	kg/m3
Flow coefficient	Cd	0.82	/
Valve seat diameter	ds	10	mm
Valve ball diameter	r s	12.7	mm
Return spring stiffness	k	18	N/mm
Valve opening	Xy	0 to 1.2	mm
The closed-loop control has been used to improve the control precision and anti-interference ability of the valve. The control schematic is shown in Fig. 2. The control system is composed of a PC, controller, linear power amplifier (linear motor driver, AMC A12A100), displacement sensor and the water hydraulic proportional valve, etc. The PC sends out the control command, and the controller compares the control command with the feedback displacement signal. Through the comparison operations, ±10V control signal will be output to the VCM linear driver (AMC A12A100). The linear driver linearly amplifies the control signal into 0 A to 10 A current to drive VCM so as to drive the spool to move until reaching the specified position. At that point, the valve spool is finally kept at the specified position.
2 THE VALVE MODEL
As the intermediate structure of transfer force and displacement, the lever amplifier significantly influences the performance of the valve. Based on the VCM model and the valve model, the influence of lever amplifier parameters on the dynamic performance is studied in detail.
2.1 VCM model
Fig. 3 shows the principle of VCM, and Table 2 shows the main VCM parameters. It is composed of the support, coil, magnet yoke, soft iron and permanent magnet. The real permanent magnetic material and soft iron material are N40H and DT3, respectively. The permanent magnet provides a uniform and constant magnetic field inside the motor. According to Lorentz force law, the Lorentz force will be produced when the coil is electrified in the magnetic field. The Lorentz force is the actuation force for the VCM and drives the coil to move.
Fig. 3. VCM composition and arrangement of internal magnetic poles
The arrangement of the inner magnetic pole is shown in Fig. 3. A total of three permanent magnets are designed and installed in VCM. Among them, the two permanent magnet poles at the bottom are assembled in reverse to form a strong magnetic gathering effect at the working air gap. This will also help to reduce magnetic leakage and improve magnetic induction strength. Then the electromagnetic efficiency will be greatly improved.
Fig. 4. The VCM equivalent electrical circuit and equivalent dynamic model
Investigation on the Modeling and Dynamic Characteristics of a Novel Hydraulic Proportional Valve Driven by a Voice Coil Motor
225
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 223-234
Table 2. The main parameters of VCM
Parameter	Value	unit
Maximum voltage	24.5	V
Maximum current	5.8	A
Continuous actuation force	51.5	N
Peak actuation force	160	N
Stroke	0 to 18	mm
Force constant	27.5	/
When the coil is electrified, Lorentz force will be produced and drive the coil to move. The dynamic model of VCM can be simplified as a mass damping mechanical motion system. Therefore, VCM can be equivalent to the coil equivalent circuit model and mass damping motion model, as shown in Fig. 4. According to the force balance equation and voltage balance equation of VCM, the mathematical model can be obtained as Eq. (1).
F = BX Nia
di„
(1)
e = B l x
a	g c c
Fe = mc xc + kc xc
i R + L —
a a a 7.
dt
where Fe is the electromagnetic force of VCM [N]; kc is the damping coefficient [N/(m/s)]; ia is the current in the coil [A]; ua is the voltage of the coil armature [V]; ea is the back electromotive force generated by the coil motion [V]; La is the coil inductance [H]; mc is the mass of the coil [kg].
2.2 Valve Dynamic Model
Compared with the poppet valve, the ball valve can realize a self-centring function. Therefore, it has better sealing performance than the poppet valve. In order to obtain good sealing performance, the valve adopts the
structure of the ball valve. The flow rate of the valve can be calculated with the following equations:
q = CA( xv)
2Aps
A( xv) = Kds \
(1+2r} xv
(2)
I + (hs0 + xv)2
where Cd is flow coefficient of valve orifice; Aps is the pressure difference of the valve [MPa]; pis the water density [kg/m3]; ds is the diameter of the valve seat [m]; Xy is the valve opening [m]; hs0 can be calculated with the following formula:
hs0 =■
(3)
When the spool is moving, the valve will be subjected to the comprehensive effect of flow force, spring force, friction force, inertia force, and viscous resistance (as shown in Fig. 5). VCM needs to provide enough driving force to overcome these resistances so that the valve spool has enough acceleration to meet the requirements of rapid response. Assuming that the opening direction of the valve is positive, the kinetic equation can be established as follows:
m x = F H—— (d2 - d5) +
S V	S	4
XPk 4
(d5 - d2) - Fa - f - Fksi
Fksi = K(L + xv) Fa	d^4ro2 - 0.25Ah2 Apo
f =pqY v c°s«o
(4)
226
Fig. 5. The force analysis of the valve
Han, M. - Liu, Y. - Liao, Y. - Wang, S.
2
2
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 223-234
where ms is equivalent mass of the valve moving parts [kg]; ps is the inlet pressure [MPa]; di is the diameter of the valve stem (i =1, 2, 3, 4, 5) [m]; pk is the pressure acting on the stem (MPa); Ff1 is the friction force [N]; Fks1 is then spring force [N]; Ffv1 is flow force (N]; ksv is spring stiffness [N/m]; lv0 is initial compression of spring [m]; p is the friction coefficient of o-ring; ro is the o-ring radius [m]; Ah is O-ring compression [m]; Apo is difference pressure acting on o-ring [MPa]; qv is the flow rate [l/min]; v is the valve orifice flow velocity [m/s]; a0 is the valve orifice jet angle [°]. Where Fs is the driving force acting on the ball, and it can be got by:
F =
x h1
Fe cos(arctan(cot a---)) -s-	xc < hsl cot a
e hs1 hs2	c s
x	h 1
Fe cos(arctan(—- - cot a )) —s-	xc > hs1 cot a
e h K	° s
, (5)
In the same way, the relationship between the stroke of VCM and the stroke of the valve spool can be obtained:
xh.
(cos(arctan(cot a--))hs1)
x < h, cot a
x h,
, (6)
x > h, cot a
(cos(arctan(--cot a))hs1 )
h, s
where Fe is Lorentz force of VCM [N]; a is the swing angle of leverage [°]; hs1 is the centre distance of VCM [m]; hs2 is the centre distance of ball [m]; xc is the coil displacement [m].
3 SIMULATION ANALYSIS
Based on the mathematical model of the valve, the simulation model can be established in the MATLAB/ Simulink and Maxwell 3D, respectively. The influence of the lever parameters on the performance of the output actuation force and displacement is compared and analysed, following which the best amplification coefficient is obtained.
3.1 Electromagnetic Simulation
The simulation model of VCM is established in Maxwell 3D, as shown in Fig. 6. Since the VCM is of symmetrical structure, half of the model can be used for simulation calculation to improve the calculation
efficiency. The main electromagnetic simulation parameters of VCM is shown in Table 3.
Fig. 6. VCM FEM model
Table 3. Electromagnetic simulation parameters
Name Soft iron
Permanent magnet
Coil
Parameter
Material: 10#
Material: N40H; Coercivity: 970 kA/m; Remanence: 1.3 mT
Copper wire diameter:0.5 mm; Turns - number: 550; Inside diameter of coil:41 mm; Outside diameter of coil:50 mm; Coil length:50.6 mm_
e) xc=\2 mill	I) xc= I iTmriT
Fig. 7. Electromagnetic simulation of VCM in the whole stroke, a) xc = 0 mm, b) xc = 2 mm, c) xc = 4 mm, d) xc = 6 mm, e) xc = 12 mm, f) xc = 18 mm
The magnetic field distribution can be obtained by the electromagnetic simulation, as shown in Fig. 7. The simulation results indicate that there will be a small amount of magnetic leakage at the top of the coil and the outside of the magnetic steel, respectively. The effect of magnetic concentration is obvious when the two permanent magnets are assembled in reverse,
si
x

Investigation on the Modeling and Dynamic Characteristics of a Novel Hydraulic Proportional Valve Driven by a Voice Coil Motor
227
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 223-234
and the majority of the Lorentz force is generated here. In the whole stroke of VCM, the maximum vector magnetic potential increases slightly from 1.6 x 10-4 WB/m to 1.6171 x 10-4 WB/m; this ensures the stability of the actuation force output of the voice coil motor in the stroke.
Fig. 8. Simulation results of VCM actuation force
Fig. 8 shows the simulation results of the VCM actuation force. The Lorentz force Fe with the variation of parameters (Uc, xc) has been shown in Fig. 8, of which the coil current ia is determined by coil control voltage Uc and xc is the coil displacement. With the movement of the coil, the maximum actuation force is generated in the middle section, and the minimum actuation force is generated at both sides of the stroke. The distribution of electromagnetic actuation force is small on both sides and large in the middle section. When the continuous actuation force is 51.5 N, the current is about 1.673 A to 1.804 A.
3.2 Characteristics Analysis of Lever Magnifier
When VCM pushes the lever to drive the spool, the vertical distance between the VCM push rod and the lever fulcrum will change with the rotation of the lever. When the lever is perpendicular to
5	10
aj	VCM sroke xc [mm]	¡jj
Fig. 9. Relationship between a) VCM stroke and spool stroke; and b) residual sum of squares
Fig. 10. Relationship between VCM stroke and magnification; a) hs1 = 40 mm, b) hs1 = 50 mm, c) hs1 = 60 mm, a) hs1 = 70 mm, a) hs1 = 80 mm, a) hs1 = 90 mm
228
Han, M. - Liu, Y. - Liao, Y. - Wang, S.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 223-234
the horizontal direction, the magnification is the maximum. The variation of magnification will lead to the non-linearity of actuation force and stroke that acts on the spool. Therefore, it is necessary to select the reasonable distance hs1 to obtain the linear output performance.
Fig. 9a shows the relationship between VCM stroke and spool stroke when hs1 is 40 mm to 90 mm. Fig. 9b is the residual sum of squares (SSE) of the spool stroke in Fig. 9a. With the increase of hs1, SSE decreases gradually. When hs1 = 70 mm, a good linear relationship can be obtained. The magnification is between 5.44 and 5.46, and the fluctuation amplitude of magnification is 0.0208, as shown in Fig. 10.
The results show that with the increase of amplification, the amplification fluctuation (A^t is the amplification fluctuation, as shown in Fig. 10) will gradually increase from 0.0117 to 0.0263. Based on the comprehensive comparison and analysis, the optimized value of hs1 is 70 mm. and the optimal magnification can be obtained.
4 EXPERIMENTS AND DISCUSSION
A prototype of the water hydraulic proportional valve is developed. The dynamic performance test rig of
the valve and VCM have been built, respectively. The dynamic and static performances of VCM and the valve are tested.
4.1 Experimental Study on the Performance of VCM
As shown in Fig. 11a, the actuation force test rig of VCM is mainly composed of support frame, force sensor, displacement sensor, step motor, ball screw, adjustable linear power supply, computer, etc. Fig. 11b shows the schematic diagram of the test system. The industrial computer sends control commands to the stepper motor driver and VCM linear driver through the acquisition card, respectively. The stepper motor drives the linear motion platform to move with VCM through the ball screw. The force sensor is fixed on the base and remains stationary. Therefore, by moving the linear motion platform and adjusting the control signal, the VCM actuation force under different stroke can be tested. The displacement sensor acquisitions the displacement signal, and the force sensor records the VCM actuation force signal every 1 mm interval.
The results presented in Figs. 12 and 13 show the actual test results of VCM actuation force and the difference between simulation and test results. With the movement of the coil, the maximum actuation force
a)
b)
Fig. 11. VCM actuation force testing rig; a) Testing rig of VCM, b) Schematic diagram of test system
' y 02»
Fig. 12. VCM actuation force test results
Fig. 13. Difference between simulation and test of VCM actuation force
Investigation on the Modeling and Dynamic Characteristics of a Novel Hydraulic Proportional Valve Driven by a Voice Coil Motor
229
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 223-234
is generated in the middle section, and the minimum actuation force is generated at both sides of the stroke. The distribution of electromagnetic actuation force is small on both sides and large in the middle section. The distribution characteristics of actuation force are essentially consistent with the simulation results. The simulation value is slightly larger than the real test value, and the maximum difference between them is nearly 8 N (the maximum error is 8.4 %). The distribution of error is irregular and significantly fluctuates with the increase of actuation force and displacement. However, in general, the error increases with the increase of actuation force. Both simulation results and test results show that the actuation force decreases in the initial stroke (0 mm to 4 mm) and the end stroke (15 mm to 18 mm). Therefore, the optimal stroke range is 4 mm to 15 mm.
4.2 Experimental Study on Dynamic Performance of the Valve
As shown in Fig. 14, the test system is mainly composed of pump, relief valve, pressure gauges, pressure transducer, displacement transducer, throttle valve, hydraulic multimeter, etc. During the test, the water pump is used to generate rated pressure and rated flow for the prototype valve. The test pressure can be set by adjusting the relief valve.
O&V) tWULLpQ
r
Fig. 14. Principle of the hydraulic test system; 1. Relief valve, 2. Water pump, 3. Filter, 4. Pressure gage, 5. Water hydraulic valve, 6. Flowmeter, 7. Throttle valve, 8. Pressure transducer, 9. Water tank
Pressure sensors are installed at the inlet and outlet of the valve. The two different ranges of the pressure sensors are 0 MPa to 20 MPa and 0 MPa to 35 MPa, respectively. The output signals of the pressure
sensors are both 4 mA to 20 mA. The range of the eddy current displacement sensor is 0 mm to 2 mm. The output signal of the eddy current displacement sensor is 4 mA to 20 mA. The pressure sensors are used to measure the pressure at the inlet and outlet of the valve. The pressure difference (Ap) can be obtained. The eddy current displacement sensor is installed on the valve and used to measure the valve spool stroke (xv). The PC sends the control commands to the controller. The controller sent control signals into the linear driver (AMC A12A100), and then the control signals will be linearly amplified to drive VCM, in order to control the movement of the valve spool. The displacement signal and pressure signals can be collected in real time with a hydraulic multimeter.
According to the above principle of the hydraulic test system, the test rig is built (as shown in Fig. 15, which is composed of hydraulic power supply, the valve test bench, DC power supply and electrical control system. The linear DC power supplies power to the controller, VCM and sensors. The electric control system realizes the closed-loop control position.
Elcctrical control system ™'	DC power supply
Fig. 15. The valve dynamic performance test system
As shown in Fig. 16, the step response (the spool displacement of 1.0 mm, 0.5 mm and 0.1 mm) is compared with the corresponding simulation results. The simulation results show that the step response adjusting time of the valve is about 28 ms and the maximum overshoot is about 5 %. The step response rise time is about 15 ms. Furthermore, the test results indicate that the adjusting time of the opening process is about 30 ms, and the adjusting time of the closing process is about 35 ms. The test results show that there is no obvious delay in the process of opening and closing. It takes more time for the valve spool to close than to open. The actuation force is provided
230
Han, M. - Liu, Y. - Liao, Y. - Wang, S.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 223-234
by VCM when the valve is opening. As the return spring closes the valve when the VCM returns. This is the main reason for the difference between the opening and closing time. The simulation results are in good agreement with the experimental results. An approximately 10 % overshoot occurs when the spool displacement is 1.0 mm, while there is no overshoot when the spool displacement is 0.5 mm and 0.1 mm. The rising time is about 12 ms to 15 ms. The valve has good dynamic performance, and the valve opening and closing process is controlled smoothly.
) response lime [
Fig. 16. Step response test and simulation
Control signal U [V]
Fig. 17. Control signal and flow characteristics
The control signal-spool displacement curve and control signal-flow characteristic curve are shown in Fig. 17. During the test, the inlet pressure of the valve is kept at 1 MPa by adjusting the relief valve. The outlet is connected to the water tank, and the pressure difference between the inlet and outlet of the valve is kept at 1 MPa. The maximum displacement of the valve spool is 1.2 mm in the range of 0 V to 10 V control signal. The control signal has the good linear relationship with the displacement of the valve spool, and the position hysteresis of the valve is less than 1 %. By comparing the experimental flow rate
with the simulation results, it can be found that the test flow rate is in good agreement with the simulation flow rate when the control signal is less than 6.5 V. However, the greater the flow rate, the higher the deviation. When the control signal is more than 6.5 V, there is a certain deviation between them, which is mainly caused by the machining error of the valve orifice and the pressure fluctuation of the test system. Generally, the control signal of the valve has a good linear relationship with the flow rate.
Ap [MPa]
Fig. 18. q and Ap characteristic
Fig. 19. The valve input power characteristics
Fig. 18 shows the pressure difference (Ap) and flow rate (qv) characteristic under different valve opening (xv = 0.2 mm, 0.4 mm, 0.6 mm, 0.8 mm, 1.0 mm) in detail. There is irregular change when Ap is over 4 MPa at the valve opening xv = 0.2mm. This is probably caused by the machining errors at the valve orifice. In general, the results indicate that the Ap and qv characteristics curve is smooth and linear. The valve has good control performance. The relationship between the flow rate and VCM power under different valve opening is shown in Fig. 19. With the increase
Investigation on the Modeling and Dynamic Characteristics of a Novel Hydraulic Proportional Valve Driven by a Voice Coil Motor
231
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 223-234

o.
00 ,
1.025 h'.ZZZZ 0.9R5	Displacement — Pressure - - k	
		\
' 0.854 h—^	I4" 28ms 33ms -0.882 i52, 4 \	4— " 4.558 _
Time t [s]
a)
E
£
tl.O-
fo^
o e.
t/>
0.0
UI53
I.021
Displacement — Pressure ___
0.679
jr/v''1—'
<- 31 ms 0.71
3 6m s
1,957
1.993

2.0 I J*
s «
o.s| 0.0

IO.sh
o o
CL
to
0.0
P053
i7(âï"
Displacement Pressure
0,667
- 30ms
0.697
/
35ms -< 2.445
12.480
14 „
I2r
— 1.0-
10

1 Time 1 [s]2
c)
£0.5
"O
"c
o
"Vo
b ' A Time r Jsf b)		
_ L065_	Displacement—pressuré	
■ 1.011	r------—	
—i 1.03! —J*	1 /' if- 33ms 38ms-* \\m 3.08 * X	¡3.118 ^-1„ mu i
6£
2
P
VI t/5
4J
It 0
22 20
16 a
3
2 Time 1 [s] ^
d;
Fig. 20. Step response under different pressure, a) Ap = 1 MPa, b) Ap = 4 MPa, c) Ap = 10 MPa, d) Ap = 16 MPa
Table 4. Step response characteristics comparison
Manufacturer	MOOG	Tested valve	ATOS	Danfoss
Model	D633/D634	Prototype	QVKZOR-A*-10	VOH30PE
Type	Direct drive servo valve	Proportional valve actuated by VCM	Proportional valve actuated by solenoid	Proportional valve actuated by solenoid
Hydraulic medium	Mineral oil	Water	Mineral oil	Water
Max Pressure	35 MPa	25 MPa	21 MPa	14 MPa
Response time (0 % to 100 % Stroke)	<12 ms <20 ms	30 ms	45 ms	<150 ms
of the flow rate, the power required to drive the spool will also increase. When the flow rate is 70 l/min, the corresponding VCM power is 60 W, which can meet the requirements of long-term stable and reliable operation of the valve.
The step response test of the valve under different pressures is shown in Fig. 20. During the test, the inlet pressure of the valve is maintained at 1 MPa, 4 MPa, 10 MPa and 16 MPa by adjusting the relief valve, while the outlet is connected with the water tank. Thus, the outlet pressure of the valve is always kept at 0. According to the test results, it can be seen that the maximum overshoot will increase slightly with the increase of the inlet pressure. Nevertheless, the valve can maintain the rapid response performance. The maximum overshoot is approximately 10 %. In general, the adjusting time of the valve is about 30
232
ms in the opening process and 35 ms in the closing process.
Table 4 shows the comparison between the new water hydraulic valve and the other proportional valves in the market. The step response characteristics of the new valve presented in this paper are better than those of the similar type of oil hydraulic proportional valve. While compared with the servo valve, there is a gap in the dynamic performance. However, compared with the traditional water hydraulic proportional valve (such as Danfoss V0H30PE), the new water hydraulic valve actuated by VCM has better dynamic performance. The results prove that VCM is a good solution instead of solenoids for the hydraulic control valve.
Han, M. - Liu, Y. - Liao, Y. - Wang, S.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 223-234
5 CONCLUSIONS
In this study, a water hydraulic proportional valve with fast response is proposed. Given that VCM has the advantages of high speed and high control accuracy, it has been used as the electrical-mechanical conversion device for the valve. Due to the poor viscosity and lubricity of the water medium, the valve will need a large force to push the spool. Thus, a lever amplifier is used to amplify the VCM actuation force. The position feedback closed-loop control has been used to improve the dynamic performance and anti-interference ability of the valve.
A detailed mathematical model of the valve has been developed. The simulation models are established in the MATLAB/Simulink and Maxwell 3D platform, respectively. A comprehensive optimization design method has been proposed. The test rig of the valve and VCM have been built. The dynamic and static performances of VCM and the valve have been tested. Both the VCM actuation force test and simulation results show that the optimal stroke range is 4 mm to 15 mm. According to the dynamic response test, the maximum overshoot of the valve is approximately 10 %, the adjusting time is about 30 ms in the opening process, and 35 ms in the closing process. The test results prove that the valve can maintain a fast response speed under different pressures. Compared with the traditional water hydraulic proportional valve, the new water hydraulic valve actuated by VCM significantly improved the dynamic performance. VCM is a good solution instead of solenoids for the water hydraulic proportional valve. The water hydraulic proportional valve designed in this paper has good static and dynamic control performance.
6 ACKNOWLEDGEMENTS
This work was supported by the National Key R&D Program of China (2018YFB2004001). Project 2662019QD023 supported Fundamental Research Funds for the Central Universities.
7 REFERENCES
[1]	Liao, Y.Y., Lian Z.S., Feng J.L., Yuan, H.B., Zhao, R.H. (2018). Effects of multiple factors on water hammer induced by a large flow directional valve. Strojniški vestnik - Journal of Mechanical Engineering, vol. 64, no. 5, p. 329-338, DOI:10.5545/sv-jme.2017.5109.
[2]	Tandon, S., Divi, S., Muglia, M., Vermillion, C., Mazzoleni, A. (2019). Modeling and dynamic analysis of a mobile underwater turbine system for harvesting marine hydrokinetic
energy. Ocean Engineering, vol. 187, art. ID. 106069, DOI:10.1016/j.oceaneng.2019.05.051.
[3]	Woodacre, J.K., Bauer, R.J., Irani, R. (2018). Hydraulic valve-based active-heave compensation using a model-predictive controller with non-linear valve compensations. Ocean Engineering, vol. 152, p. 47-56, DOI:10.1016/j. oceaneng.2018.01.030.
[4]	Liao, Y.Y., Yuan, H.B., Lian, Z.S., Feng, J.L., Guo, Y.C. (2015). Research and analysis of the hysteresis characteristics of a large flow directional valve. Strojniški vestnik - Journal of Mechanical Engineering, vol. 61, no. 6, p. 355-364, DOI:10.5545/sv-jme.2015.2487.
[5]	Han, M.X., Liu, Y.S., Wu, D.F., Zhao, X.F., Tan, H.J. (2017). A numerical investigation in characteristics of flow force under cavitation state inside the water hydraulic poppet valves. International Journal of Heat and Mass Transfer, vol. 111, p. 1-16, DOI:10.1016/j.ijheatmasstransfer.2017.03.100.
[6]	Majdic, F., Pezdirnik, J. (2010). Oil- and water-based continuous control valve. Industrial Lubrication and Tribology, vol. 62, no. 3, p. 136-143, DOI:10.1108/00368791011034511.
[7]	Han, M.X., Liu, Y.S., Zheng K., Ding, Y.C., Wu, D.F. (2020). Investigation on the modeling and dynamic characteristics of a fast-response and large-flow water hydraulic proportional cartridge valve. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 234, no. 22, p. 4415-4432, DOI:10.1177/0954406220922860.
[8]	Zhu, K.W., Gu, L.Y., Chen, Y.J., Li, W. (2012). High speed on/off valve control hydraulic propeller. Chinese Journal of Mechanical Engineering, vol. 25, no. 3, p. 463-473, DOI:10.3901/CJME.2012.03.463.
[9]	Mahrenholz, J., Lumkes, J.J. (2010). Analytical coupled modeling and model validation of hydraulic on/off valves. Journal of Dynamics Systems, Measurement, and Control, vol. 132, no. 1, p. 1-10, DOI:10.1115/1.4000072.
[10]	Park, S.-H., Kitagawa, A., Kawashima, M., (2004). Water hydraulic high-speed solenoid valve. Part I: development and static behavior. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, vol. 218, no. 5, p. 399-499, DOI:10.1243/0959651041568560.
[11]	Park, S.-H. (2009). Development of a proportional poppet-type water hydraulic valve. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 223, no. 9, p. 2099-2107, DOI:10.1243/09544062JMES1380.
[12]	Roberts, D.C., Li, H., Steyn, J.L., Yaglioglu, O., Spearing, S.M., Schmidt, M.A., Hagood, N.W. (2003). A piezoelectric microvalve for compact high-frequency, high-differential pressure hydraulic micropumping systems. Journal of Microelectromechanical Systems, vol. 12, no. 1, p. 81-92, DOI:10.1109/JMEMS.2002.807471.
[13]	Rogge, T., Rummler, Z., Schomburg, W.K. (2004). Polymer micro-valve with a hydraulic piezo-drive fabricated by the AMANDA process. Sensors and Actuators A: Physical, vol. 110, no. 1-3, p. 206-212, DOI:10.1016/j.sna.2003.10.056.
[14]	Ahn, D.H., Hong, D.P., Gweon, D.G. (2014). Design of a high force voice coil motor. Applied Mechanics and Materials,
Investigation on the Modeling and Dynamic Characteristics of a Novel Hydraulic Proportional Valve Driven by a Voice Coil Motor
233
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 223-234
vol. 483, p. 559-562, D0l:10.4028/www.scientific.net/ AMM.483.559.
[15]	Li, B.R., Gao, L.L., Yang, G. (2013). Modeling and control of a novel high-pressure pneumatic servo valve direct-driven by voice coil motor. Journal of Dynamic Systems, Measurement, and Control, vol. 135, no. 1, p. 1-5, D0I:10.1115/1.4007702.
[16]	Li, B.R., Gao, L.L., Yang, G. (2013). Evaluation and compensation of steady gas flow force on the high-pressure electro-pneumatic servo valve direct-driven by voice coil motor. Energy Conversion and Management, vol. 67, p. 92-102, D0I:10.1016/j.enconman.2012.11.004.
[17]	Zhang, Z.M., Gong, Y.J. Hou, J.Y., Wu, H.P. (2014). Simulation on linear-motor-driven water hydraulic reciprocating plunger pump. Advanced Materials Research, vol. 842, p. 530-535, D0I:10.4028/www.scientific.net/AMR.842.530.
[18]	Liu, W., Wei, J.H., Fang, J.H., Li, S.Z. (2015). Hydraulic-feedback proportional valve design for construction machinery. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 229, no. 17, p. 3162-3178, D0I:10.1177/0954406214568822.
[19]	Filo, G., Lisowski, E., Kwiatkowski, D., Rajda J., (2019). Numerical and experimental study of a novel valve using the return stream energy to adjust the speed of a hydraulic actuator. Strojniški vestnik - Journal of Mechanical Engineering, vol. 65, no. 2, p. 103-112, D0I:10.5545/sv-jme.2018.5823.
[20]	Xu, Y. L., Nie, H.W., Zhao, H.L., Liu, J. (2020). Mathematical modelling and simulation of a novel hydraulic variable valve timing system. International Journal of Simulation Modelling, vol. 19, no. 2, p. 303-312, D0I:10.2507/IJSIMM19-2-C06.
[21]	Simic, M., Herakovic, N. (2015). Reduction of the flow forces in a small hydraulic seat valve as alternative approach to improve the valve characteristics. Energy Conversion and Management, vol. 89, p. 708-718, D0I:10.1016/j. enconman.2014.10.037.
[22]	Karanovic, V., Jocanovic, M., Baloš, S., Kneževic, D., Mačužic, I. (2019). Impact of contaminated fluid on the working performances of hydraulic directional control valves. Strojniški vestnik - Journal of Mechanical Engineering, vol. 65, no. 3, p. 139-147, D0I:10.5545/sv-jme.2018.5856.
234
Han, M. - Liu, Y. - Liao, Y. - Wang, S.
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 235-244 © 2021 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2021.7113	Original Scientific Paper
Received for review: 2021-01-28 Received revised form: 2021-04-13 Accepted for publication: 2021-05-04
Virtual Minimization of Residual Stress and Deflection Error in the Five-Axis Milling of Turbine Blades
Mohsen Soori* - Mohammed Asmael
Eastern Mediterranean University, Department of Mechanical Engineering, Turkey
To simulate and analyse the real machined parts in virtual environments, virtual machining systems are applied to the production processes. Due to friction, chip forming, and the heat produced in the cutting zone, parts produced using machining operation have residual stress effects. The machining force and machining temperature can cause the deflection error in the machined turbine blades, which should be minimized to increase the accuracy of machined blades. To minimize the residual stress and deflection error of machined parts, optimized machining parameters can be obtained. In the present research work, the application of a virtual machining system is presented to predict and minimize the residual stress and deflection error in a five-axis milling operations of turbine blades. In order to predict the residual stress and deflection error in machined turbine blades, finite element analysis is implemented. Moreover, to minimize the residual stress and deflection error in machined turbine blades, optimized parameters of machining operations are obtained by using a genetic algorithm. To validate the research work, experimentally determining residual stress by using a X-ray diffraction method from the machined turbine blades is compared with the finite element results obtained from the virtual machining system. Also, in order to obtain the deflection error, the machined blades are measured by using the CMM machines. Thus, the accuracy and reliability of machined turbine blades can be increased by analysing and minimizing the residual stress and deflection error in virtual environments.
Keywords: virtual machining, residual stress, deflection error, parameter optimization, turbine blade Highlights
•	The application of a virtual machining system is presented to predict and minimize the residual stress and deflection error in the five-axis milling operations of turbine blades.
•	In order to minimize the residual stress and deflection error in machined turbine blades, an optimization technique based on genetic algorithms is used.
•	The accuracy and reliability of machined turbine blades can be increased by analysing and minimizing the residual stress and deflection error in virtual environments.
0 INTRODUCTION
The residual stress in machined parts can be generated as a result of mechanical, thermal, and chemical effects in chip forming of metal cutting operations. The generated residual stress in the machined parts can impair the performance of components, such as fatigue life, corrosion resistance, and part distortion in actual working conditions. Due to the cutting forces and cutting temperature, the machined blades have deflection errors, which can cause inaccuracy in the machined turbine blades. As a result, the residual stress and deflection error should be analysed and decreased to increase the accuracy and reliability of parts produced using machining operations.
To improve the precision of machined parts, the evaluation of surface error characteristics in thin-walled constructs during peripheral milling has been studied by Wimmer and Zaeh [1]. The process parameter optimization of thin-wall machining for wire arc additive manufactured parts is investigated by Grossi et al. [2] to decrease the deformation error in the thin-walled manufactured components. The method of error compensation during the milling operations
of flexible thin-wall parts is presented by Ratchev et al. [3] to reduce the deflection error of machined parts. To improve the accuracy of thin-walled machined parts, an adaptive toolpath methodology for three-axis milling is presented by Grossi et al. [4]. The finite element method (FEM) based cutting velocity selection for thin-walled part machining is presented by Scippa et el. [5] to enhance the precision of machined components. Finite-Element modelling of workpiece vibrations is analysed by Bolsunovskiy et al. [6] to optimize the machining parameters and increase accuracy in the milling operations of thin-walled components.
Jiang et al. [7] analysed the effects of cutting forces and cutting zone temperature to the residual stress of machined components. The application of response surface methodology in obtaining the optimized machining parameters, such as depth of cut and spindle speed, are investigated by Masmiati et al.
[8]	in order to minimize the residual stress, cutting force, and surface roughness in the end milling of S50C medium carbon steel. Mohammadpour et al.
[9]	investigate the effects of machining parameters on the residual stress of machined parts in milling
*Corr. Author's Address: Department of Mechanical Engineering, Eastern Mediterranean University, Famagusta", North Cyprus, Mersin 10, Turkey, mohsen.soori@emu.edu.tr
235
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 235-244
operations. To analyse the influence of the machining parameters, such as cutting speed as well as feed rate, on the surface roughness and residual stresses in produced parts using milling operation, the finite element method is used by Zhang and Wu [10].
To decrease the residual stress in machined parts using turning operations, the effects of cutting tool parameters, such as tool nose radius and tool wear on residual stresses, are investigated by Lin et al. [11]. Cutting tool parameters, such as tool edge radius, are investigated by Yang et al. [12] in order to decrease the residual stress in machined parts. To study the effects of cutting-edge radius and cutting forces on residual stresses of machined parts, a finite element model is developed by Nasr et al. [13].
To decrease the residual stress in machined parts and effects of cutting tools materials, machining parameters, such as cutting speeds and depth of cut to residual stress of machined parts, are investigated by Arunachalam et al. [14]. The influence of cutting tool parameters, such as tool nose radius on surface integrity and residual stresses, are investigated by Sharman et al. [15] to decrease the residual stress in machined parts. A variable depth-of-cut milling strategy for thin-walled workpiece is proposed by Yan et al. [16] in order to decrease the deflection error in the machined turbine blades. A five-axis adaptive flank milling of flexible thin-walled parts based on the on-machine measurement is developed by Huang et al. [17] to modify the cutting tool paths in terms of deflection error minimization. The time domain flutter analysis based on the blade element momentum (BEM) theory for bend-twist coupled large composite wind turbine blades is investigated by Shakya et al. [18] to increase the performances of wind turbines in actual working conditions. The blisk vibration in aircraft engine using similitude models is analysed by Luo et al. [19] to predict and decrease the vibration of the engine systems in the virtual environments. To analyse and modify the machining operations in virtual environments, virtual machining systems and applications are presented by Soori et al. [20] to [24]. To analyse and develop the process of part manufacturing using welding operations, a review in recent development of friction stir welding process is presented by Soori et al. [25]. To increase efficiency in process of part production using the friction stir spot welding, effects of machining parameters to the mechanical properties of aluminium sheet alloys is presented by Nasir et al. [26].
According to the analysis of previous published papers, the area of residual stress as well as deflection error prediction and minimization in milling
operations by using virtual machining systems has been insufficiently explored. Moreover, to minimize the residual stress and deflection errors in machined parts, the applications of the virtual machining system is not studied.
In this research work, the application of a virtual machining system is presented to predict and minimize the residual stress and deflection error in a five-axis milling operations of turbine blades. The cutting forces as well as cutting temperatures for each position of cutting tool along machining paths are calculated in order to obtain residual stress and deflection error of blades due to machining operations by using the finite element analysis (FEA). The optimization technique of a genetic algorithm is used to calculate the optimized machining parameters in terms of residual stress and deflection error minimization. To measure the residual stresses on the surface and in-depth of the machined blades, the X-ray diffraction technique is used in the study. The coordinate measuring machines (CMM) machines are sued to obtain the deflection error in the machined blades. Finally, the obtained data are compared by using diagrams.
1 CUTTING FORCE MODEL
To calculate the cutting forces in the five-axis machining operations, the cutting force model developed by Zhang et al. [27] is applied.
The XYZ planes are defined in order to describe the tool motion process, workpiece geometry and tool path in machining operations. Also, XrYrZr is created in order to describe the tool rotation coordinate system in machining operations. Thus, the cutting forces in the differential format, which is applied to the jth axial disk element of the ith tooth at rotation angle 0,;(t), can be presented as:
' dFT ., (t) = (Ktsh(i, j, t)db + Ktpdb)W(0(0,j(t))) < dFRi:j ( t) = (Krsh(i, j, t)db + Kpdb)W(0(0,j ( t))), (1)
dFA . j (t) = (Kash(i, j, t)db + Kpdb)W(0(0,j (t)))
where h(i,j,t) is the instantaneous uncut chip thickness of the jth axial disk element of the ith tooth at t moment and the db is the axial height of the cutting disk element. Moreover, the coefficients of shearing specific cutting force in different directions as tangential, radial, and axial can be shown as:
kqs = aqs + bqse«M'' (q = t, r, a).
236
Soori, M. - Asmael, M.
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 235-244
Also, the coefficients of ploughing specific cutting force in different directions as tangential, radial, and axial direction can be presented as:
kqp = aqp + KeC,'h(t) (q = t, r, a).
To determine whether the condition of disk element is in or out of cutting operation, the W(0) can be presented as:
fi e,. - <d(é. . (t)) <e . -MfKh - (t))) = ]n stiJ (V,-J ()) J— (2)
[ 0	others
where 8exÂJ stand for entry and exit angle of jth disk element on ith tooth respectively, as shown in Fig. 1 [27].
Fig. 1. Entry and exit angle model of cutting tool [27]
Thus, the overall cutting forces that are applied to the t moment can be calculated by summing three cutting forces acting on disk element within axial depth of cut as is shown in Eq. (3).
Fc, F
(t ) C,, (t)
F ,z (t )
= Z, I Rot
dFT dF„
,j (t ) ,, (t )
dFA,, (t)
(3)
where Rot is the cutter rotation matrix. With the aid of transformation relationship matrix in the five-axis machine tool, the total cutting force components acting on the workpiece can be obtained by:
.x (t )
, (t )
Fw Fw Fw,z (t )
= - RotC ■ RotB
F x (t )
Fc,y (t )
F z (t )
(4)
where RotC RotB is rotation matrix for B-axis and C-axis.
2 CUTTING TEMPERATURE MODEL
In the machining operations, the chip formation process generates heat in the cutting zone. To model the cutting temperature regarding the cutting parameters, the response surface methodology (which is statistical) and mathematical techniques can be applied. The experimental values of cutting temperature corresponding to the control variables of machining parameters are incorporated into the response surface methodology in terms of cutting temperature model modification [28]. As a result, The equation of the cutting temperature regarding the cutting parameters based on response surface methodology is presented by Shaw [29] as,
1	jrbl rbl b3 b4
X = CXV J z ae ap ,
(5)
where X is temperature [°C], CX is temperature coefficient determined by workpiece material, machine tool, and cutting tool geometry parameter. The bb b2, b3, and b4 are exponents influencing the machining parameters V, cutting speed (m/min), fz, feed rate (mm/rotation), ae, radial feed [mm] and ap, axial feed [mm].
To obtain the experimental values of cutting temperature and calculate the cutting temperature during machining operation, the tool-workpiece thermocouple device with three brushes is used. The electromotive force produced in the tool-workpiece thermocouple circuit during machining is proportional to the temperature difference between the tool and room temperature. Extruded round bars (101 mm x 2000 mm) of the aluminium alloys 1350-0 and 7075-T6 were used as the workpiece. The cutting tool used in the experiment is carbide bar (310 x 10 x 4) mm K15 grade, and the tool holder is (25 x 20 x 150) mm made of steel. The cutting speed 600 m/minute, feed rate 0.3 mm/rotation and axial depth of cut 3 mm are considered as machining parameters of the experiments. The average of the chip-tool temperatures measured by the tool workpiece thermocouple system (average of two replicates) are calculated using analysis of variance of the 2k factorial design as input variables. The coefficients of the temperature model are then calculated using a statistic analysis based on multiple nonlinear regression of the chip-tool interface temperature over 32 experiments. As a result, the equations of the cutting tool and workpiece temperatures in the machining operations of Al alloys based on multiple nonlinear regression analysis are presented by Santos et al. [30] as:
Virtual Minimization of Residual Stress and Deflection Error in the Five-Axis Milling of Turbine Blades
237
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 235-244
1 1 1/C cm T/0.211 /-0.181 0.0464 0.00391
= 116.503 K fz ae a p , (6)
1 Al 1 1 n T/0.365 /-0.244 0.0423 0.019
Xw = 43.319K f a° ap , (7)
where X is temperature [°C], V is cutting speed [m/ min], f is feed rate [mm/rotation], ae and ap are radial feed and axial feed [mm], respectively.
3 DEFLECTION ERROR ANALYSIS OF THIN-WALL WORKPIECE IN MACHINING OPERATIONS
The machining operations of thin-walled parts have deflection errors as a result of cutting forces as well as cutting temperature. The workpiece deflects to a new position due to cutting forces and cutting temperatures which creates dimensional errors and inaccuracies in the process of machining operations.
The corresponding surface dimensional error can be presented as Eq. (8) [31],
e = 8
p t ,p
Jf, p'
(8)
where, Stp and f are the normal projections of the cutting force and cutting temperature-induced deflection error for the Point P, respectively.
For the convenience of study, the distance between the initial surface to be machined and the desired machined surface is named as the nominal radial depth of cut notated by RN. In the actual milling process, to ensure that the surface dimensional error does not violate the tolerance, RA is often specified to be different from RN4RA. In this case, Eq. (9) is no longer applicable to the calculation of surface dimensional error and must be corrected as [31]:
ep = S,,p +Sf,p + rn -ra.
(9)
Note that RN and RA are the nominal and specified radial depth of cut, respectively.
4 VIRTUAL MACHINING SYSTEM
The Visual Basic programming language is used in this study to develop the presented virtual machining system. The developed virtual machining system can obtain the cutting forces regarding the cutting tool details and parameters of machining operations. Thus, cutting forces at each position of the cutting tool along machining paths can be calculated.
The software calculates the cutting temperature regarding the cutting process parameters at each position of the cutting tool along machining paths. Then, the software is linked to the Abaqus R2016X FEM analysis software to analyse the residual stress and deflection error due to machining operations.
The calculated cutting forces and cutting temperature at each position of cutting tool along the machining paths in turbine blade machining are input to Abaqus software to be used for residual stress and deflection error calculations. As a result, the residual stress and deflection error due to cutting forces and cutting temperatures at each position of the cutting tool can be calculated and presented. A flowchart and the strategy of the virtual machining system in cutting force calculation, residual stress and deflection error prediction are shown in Fig. 2.
Fig. 2. Flowchart of the virtual machining system
To obtain the optimized machining parameters in terms of residual stress and deflection error minimizations, the genetic algorithm is applied. The genetic algorithm is an advanced optimization method based on natural selection by which the fittest individuals are selected for reproduction to produce offspring of the next generation. The natural process of evolution and set of chromosomes is introduced in terms of optimization process. The initial population for the optimization process is created using the binary encoding process. To provide evaluation criteria in the optimization process and rank the chromosomes in the population, the fitness function is calculated as is presented in Eq. (10) [32].
F (x) = -
1
1 + f (*)'
where fitness function and objective function are F(x) and fx), respectively. The main operators of the algorithms are reproduction, crossover, and mutation. The operators are applied to the initial population of
238
Soori, M. - Asmael, M.
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 235-244
the optimization process in order to provide a faster convergence to the optimized parameters.
The mathematical equation of the surface roughness in end milling operations is presented in Eq. (11) [33].
f
R =318-—, ° 4d
(11)
where f is the feed rate, and d is the diameter of the cutting tool. Moreover, the time of machining operation should be minimized in order to decrease the cost of machined parts. The mathematical equation of machining time is presented in Eq. (12) [33].
= A = 1'
(12)
where k is the distance of cutting tool to reach to the operational zone and f is feed rate.
Also, the cutting tool life should be maximized in the optimization process in order to decrease the cost of machining operation. The cutting tool life can be shown as Eq. (13) [33].
* -l.f
C
G
V (A)W
(13)
where Q is the contact proportion of cutting edge with workpiece per revolution, C is 33.98 for the highspeed steel (HSS) tools and 100.05 for the carbide tools, g = 0.14, V is cutting speed [mm/min], w = 0.28, m is 0.15 for HSS tools while it reaches a maximum of 0.30 for carbide tools, G and A are slenderness ratio and chip cross-section, which can be shown as Eqs. (14) and (15), respectively [33].
g=a, f
A = a • f,
(14)
(15)
where a is a depth of cut and f is feed rate in machining operations.
According to Eqs. (6) and (7), the cutting temperature can be decreased by reducing the feed rate and spindle speed. Also, according to Eq. (11), the surface roughness can be decreased by reducing the feed rate. However, the machining time will be increased by reducing the feed rate according to Eq. (12). Moreover, increasing the spindle speed can decrease the cutting forces according to Eq. (4) and cutting tool life according to Eq. (13). It is clear that the presented mathematical models of the cutting
forces, cutting temperature, time of machining operations, surface roughness, feed rate and cutting tool life are related, which should be considered as an optimization problem. The objective function of the optimization process is minimizing the residual stress and deflection error by calculating the optimized machining parameters. To obtain the optimized machining parameters in order to minimize the residual stress and deflection error using the minimized cutting forces and cutting temperature in machined turbine blade, the Matlab programming language is used. As a result, the optimized feed rate and spindle speed are obtained in order to minimize the residual stress as well as deflection error in the machined parts. The complete methodology of the study is shown in Fig. 3.
Fig. 3. Methodology of the study
5 VALIDATION AND COMPARISON
In order to validate the presented virtual machining system in the study, real turbine blades are machined by using the A 5-axis computer navigated control (CNC) milling machine tool Kondia HM 1060. The turbine blade material is AL 2618, which is used in the compressor section of gas turbines. Then, the test turbine blades are machined by using the five-axis Kondia HM 1060 CNC machine tool. The cutting tool used in the experiment is a carbide flat end mill with 10 mm diameter, helix angle 30°, flute number 4, overall length 60 mm and flute length 30 mm. The feed rate and spindle speed are 100 mm/min and 1000 rpm, respectively. The machining process of the turbine blade is shown in Fig. 4.
To calculate the cutting forces in virtual environments, the cutting force model of five-axis
5
Virtual Minimization of Residual Stress and Deflection Error in the Five-Axis Milling of Turbine Blades
239
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 235-244
CNC machine tools presented by Zhang et al. [27] is used in this study.
Fig. 4. Machining process of the turbine blades
In order to calculate the coefficients of specific cutting force, the average of cutting forces for twenty slot milling tests with 1.5 mm axial depth of cut using 5-axis CNC milling machine tool Kondia HM 1060 are measured with a Kistler dynamometer. The spindle rotating speed is 5000 rpm and the feed per tooth and feed rate are 0.5 mm and 100 mm/min, respectively. The cutting tool overhang is 60 mm. As result, the specific cutting force coefficients are obtained as Fig. 6.
The real machined Al turbine blades are shown in Fig. 5.
Fig. 5. The real machined turbine blades
The cutting forces and cutting temperatures for each position of the cutting tool along machining paths are obtained by using the virtual machining system. Then, the cutting forces and cutting temperatures are
entered into Abaqus software to obtain the residual stress and deflection error of machined blades. Thus, the calculated residual stresses by Von Mises yield criterion due to machining operation of turbine blade in MPa unit are shown in Fig. 7.
a)
o
^ 200
		
— Krp	Ktp	KUp
b)
0JÄ	OJ	0.15	0 i	'J 23
Uncut chip thickness [mm]
Fig. 6. Specific cutting force coefficients, a) Kts, Krs, Kas,
and b) Krp, Ktp Kap
Fig. 7. The calculated residual stresses by the Von Mises yield criterion due to machining operation of turbine blade in MPa unit
Fig. 8. The calculated deflection error of machined turbine blade
240
Soori, M. - Asmael, M.
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 235-244
The calculated deflection error of the machined turbine blade by using the FEM method is presented in Fig. 8.
Residual stresses of the machined blade are measured by using an X-ray diffraction method. In order to obtain the residual stress in the machined turbine blades, two-circle Bragg Brentano-Diffractometer of the type Philips X'PERT MPD is used by Cr ka radiation with 35 kV and 35 mA. The sin2^ method is the most widely used methodology for stress assessment. Several XRD measurements are performed by considering the different tilts angle of In this experiment, the tilting angles are set at 0°, 13.8°, 18.7°, 24.3°, 29.4°, 32.2°, 34.8°, 36.5°, and 40.4°. The inter-planar spacing or two-theta peak is calculated in order to obtain a plot or a curve for the measured data. As a result, the stress can be determined using the obtained plot by measuring the line gradient and with simple knowledge of the material's elastic properties. In order to measure the residual stress in the machined turbine blades, the X-ray diffraction method is used, as shown in Fig. 9.
Fig. 9. The X-ray diffraction method to measure the residual stress of machined turbine blade
The residual stresses in the depth of machined blade are measured and presented in the Fig. 10.
Fig. 10. Measured residual stress in depth of produced blade
A 90.4 % compatibility is obtained in the comparison of the results of the experimental test and FEM simulation. Optimized machining parameters are calculated using the optimization techniques developed in the study based on the genetic algorithm. To measure the cutting forces in the machining operation, the Kistler dynamometer is used. To provide input data for the optimization process, the measured cutting forces and cutting temperature, cutting condition and cutting tool geometry are entered in the optimization algorithm. In the optimization process, a population size of 28 is selected with the iterated for 228 generations. Also, the probability of crossover of 0.85 and mutation of 0.001 are selected. As a result, the feed rate 87 mm/min, radial feed ae 0.2 mm and axial feed ap 12 mm and the spindle speed 1269 rpm are obtained. Measured residual stress in the depth of machined blade using optimized machining parameters is shown in Fig. 11.
Fig. 11. Measured residual stress in depth of machined blade using optimized machining parameters
In order to measure the deflection error in the machined turbine blades, the CMM machine is used. The used probe in the error measuring is the Renishaw RSH 250 probe while its repeated accuracy is 1 ^m in touching directions. The sweep scan method is used to obtain the surface condition of the machined blade.
A NURBS surface "patch" is generated for each of the sweep scans in order to obtain the deflection error in the machined turbine blades. Then, to obtain the deflection error, sections are taken along the generated surfaces at any required height. The amount of error between the CAD model of turbine blade and fitted NURBS surfaces from the measured data of the CMM machine are calculated in each section of generated surface as the deflection error of machined turbine blade.
The process of surface generation and deflection error calculation due to cutting forces and cutting
Virtual Minimization of Residual Stress and Deflection Error in the Five-Axis Milling of Turbine Blades
241
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 235-244
temperature for the machined blade are shown in Fig. 12.
Fig. 12. The deflection error calculation from the CAD model of turbine blade and generated NURBS surfaces of machined turbine blades
As a result, the measured and predicted deflection errors in the machined turbine blade without and with the optimized machining parameters for the L7 line in Fig. 12 is presented in Fig. 13.
Fig. 13. The measured and predicted deflection errors in the machined turbine blade without and with optimized machining parameters for the L7 line in the Fig. 12
The predicted deflection error due to the cutting forces and cutting temperatures without optimized machining parameters for the L7 line in Fig. 12 is presented in Fig. 14.
The optimized machining parameters can decrease the unnecessary cutting forces and the cutting temperature in the machining operation. Therefore, a reduction in the residual stress of produced blades using optimized machining parameters can be achieved, which can increase accuracy as well as reliability of machined turbine blades. Moreover, the deflection error due to machining forces and cutting temperature is minimized to increase the accuracy of machined turbine blades.
Turbine blade length [mm]
Fig. 14. The predicted deflection error due to the cutting forces and cutting temperatures without optimized machining parameters for the L7 line in the Fig. 12
6 CONCLUSIONS
To predict and minimize residual stress and deflection error in a five-axis turbine blade milling operations, a virtual machining system is developed in the study. The system can calculate the cutting forces and cutting temperature at each position of the cutting tool along machining paths with regard to the parameters of machining operations, geometries, and material properties of cutting tools as well as workpieces. Then, FEA is used to obtain residual stress and deflection error of the turbine blades in machining operations. The real turbine blades are produced and experimentally tested using a five-axis CNC milling machine tool in order to validate the developed system in the study. Then, the residual stress in the machined turbine blades are measured using X-ray diffraction machines. To measure the deflection error in the machined blades, the CMM machine is used. Finally, the obtained results from the virtual machining system and experimental tests are compared to present the rate of compatibility in the developed system. Therefore, a 90.4 % and 94.3 % compatibility are obtained in comparison to the experimental and virtual machining system results for the residual stress and deflection error, respectively.
To minimize the residual stress and deflection error in machining operation, optimized machining parameters using the genetic algorithm is obtained. As a result, 24.2 % and 23.9 % reductions in residual stress can be obtained in the real machined turbine blade and virtual machining system results. respectively. Moreover, 26.3 % and 24.1 % reductions
242
Soori, M. - Asmael, M.
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 235-244
in the deflection error of real and virtual machining systems results are respectively measured by using the optimized machining parameters in machining operations.
The accuracy and the reliability of machined turbine blades can be increased by using the application of the virtual machining system in the minimization of residual stress and deflection error due to machining operations. Furthermore, the amount of the deflection errors due to cutting forces and cutting temperatures can be accurately predicted in virtual environments using the FEM analysis. The system developed in this study can be used in the machining operations of free-form surfaces in aircraft airfoils by using the five-axis CNC machine tools in order to decrease the residual stress of machined parts. These are the concepts of future research works for the authors.
7 REFERENCES
[1]	Wimmer, S., Zaeh, M. (2018). The prediction of surface error characteristics in the peripheral milling of thin-walled structures. Journal of Manufacturing and Materials Processing, vol. 2, no. 1, art. ID 13, DOI:10.3390/jmmp2010013.
[2]	Grossi, N., Scippa, A., Venturini, G., Campatelli, G. (2020). Process parameters optimization of thin-wall machining for wire arc additive manufactured parts. Applied Sciences, vol. 10, no. 21, art ID. 7575, DOI:10.3390/app10217575.
[3]	Ratchev, S., Liu, S., Becker, A.A. (2005). Error compensation strategy in milling flexible thin-wall parts. Journal of Materials Processing Technology, vol. 162-163, p. 673-681, DOI:10.1016/j.jmatprotec.2005.02.192.
[4]	Grossi, N., Scippa, A., Croppi, L., Morelli, L., Campatelli, G. (2019). Adaptive toolpath for 3-axis milling of thin walled parts. MM Science Journal, no. 04, p. 3378-3385, DOI:10.17973/ MMSJ.2019_11_2019096.
[5]	Scippa, A., Grossi, N., Campatelli, G. (2014). FEM based cutting velocity selection for thin walled part machining. Procedia CIRP, vol. 14, p. 287-292, DOI:10.1016/j.procir.2014.03.023.
[6]	Bolsunovskiy, S., Vermel, V., Gubanov, G., Kacharava, I., Kudryashov, A. (2016). Thin-walled part machining process parameters optimization based on finite-element modeling of workpiece vibrations. Procedia CIRP, vol. 8, p. 276-280, DOI:10.1016/j.procir.2013.06.102.
[7]	Jiang, X., Li, B., Wang, L., Wang, Z., Li, H. (2016). An approach to evaluate the effect of cutting force and temperature on the residual stress generation during milling. The International Journal of Advanced Manufacturing Technology, vol. 87, p. 2305-2317, DOI:10.1007/s00170-016-8605-5.
[8]	Masmiati, N., Sarhan, A.A., Hassan, M.A.N., Hamdi, M. (2016). Optimization of cutting conditions for minimum residual stress, cutting force and surface roughness in end milling of S50C medium carbon steel. Measurement, vol. 86, p. 253265, DOI:10.1016/j.measurement.2016.02.049.
[9]	Mohammadpour, M., Razfar, M.R., Saffar, R.J. (2010). Numerical investigating the effect of machining parameters on residual stresses in orthogonal cutting. Simulation Modelling Practice and Theory, vol. 18, no. 3, p. 378-389, DOI:10.1016/j. simpat.2009.12.004.
[10]	Zhang, X., Wu, H. (2018). Effect of tool angle on cutting force and residual stress in the oblique cutting of TC21 alloy. The International Journal of Advanced Manufacturing Technology, vol. 98, p. 791-797, DOI:10.1007/s00170-018-2324-z.
[11]	Lin, K., Wang, W., Jiang, R., Xiong, Y. (2019). Effect of tool nose radius and tool wear on residual stresses distribution while turning in situ TiB2/7050 Al metal matrix composites. The International Journal of Advanced Manufacturing Technology, vol. 100, p. 143-151, DOI:10.1007/s00170-018-2742-y.
[12]	Yang, D., Xiao, X., Liang, X. (2019). Analytical modeling of residual stress in orthogonal cutting considering tool edge radius effect. The International Journal of Advanced Manufacturing Technology, vol. 103, p. 2965-2976, DOI:10.1007/s00170-019-03744-9.
[13]	Nasr, M.N.A., Ng, E.-G., Elbestawi, M.A. (2007). Modelling the effects of tool-edge radius on residual stresses when orthogonal cutting AISI 316L. International Journal of Machine Tools and Manufacture, vol. 47, no. 2, p. 401-411, DOI:10.1016/j.ijmachtools.2006.03.004.
[14]	Arunachalam, R.M., Mannan, M.A., Spowage, A.C. (2004). Residual stress and surface roughness when facing age hardened Inconel 718 with CBN and ceramic cutting tools. International Journal of Machine Tools and Manufacture, vol. 44, no. 9, p. 879-887, DOI:10.1016/j.ijmachtools.2004.02.016.
[15]	Sharman, A.R.C., Hughes, J.I., Ridgway, K. (2015). The effect of tool nose radius on surface integrity and residual stresses when turning Inconel 718™. Journal of Materials Processing Technology, vol. 216, p. 123-132, DOI:10.1016/j. jmatprotec.2014.09.002.
[16]	Yan, Q., Luo, M., Tang, K. (2018). Multi-axis variable depth-of-cut machining of thin-walled workpieces based on the workpiece deflection constraint. Computer-Aided Design, vol. 100, p. 14-29, DOI:10.1016/j.cad.2018.02.007.
[17]	Huang, N., Bi, Q., Wang, Y., Sun, C. (2014). 5-Axis adaptive flank milling of flexible thin-walled parts based on the on-machine measurement. International Journal of Machine Tools and Manufacture, vol. 84, p. 1-8, DOI:10.1016/j. ijmachtools.2014.04.004.
[18]	Shakya, P., Sunny, M.R., Maiti, D.K. (2020). Time domain flutter analysis of bend-twist coupled large composite wind turbine blades: a parametric study. Mechanics Based Design of Structures and Machines, p. 1-23, DOI:10.1080/15397734 .2020.1824796.
[19]	Luo, Z., Wang, Y., Zhai, J., Zhu, Y., Wang, D. (2019). Prediction of vibration characteristics of blisks using similitude models. Mechanics Based Design of Structures and Machines, vol. 47, no. 2, p. 121-135, DOI:10.1080/15397734.2018.1481427.
[20]	Soori, M., Arezoo, B., Habibi, M. (2017). Accuracy analysis of tool deflection error modelling in prediction of milled surfaces by a virtual machining system. International Journal of Computer Applications in Technology, vol. 55, no. 4, p. 308321, DOI:10.1504/ijcat.2017.10006843.
Virtual Minimization of Residual Stress and Deflection Error in the Five-Axis Milling of Turbine Blades
243
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 235-244
[21]	Soori, M., Arezoo, B., Habibi, M. (2014). Virtual machining considering dimensional, geometrical and tool deflection errors in three-axis CNC milling machines. Journal of Manufacturing Systems, vol. 33, no. 4, p. 498-507, D0I:10.1016/j.jmsy.2014.04.007.
[22]	Soori, M., Arezoo, B., Habibi, M. (2013). Dimensional and geometrical errors of three-axis CNC milling machines in a virtual machining system. Computer-Aided Design, vol. 45, no. 11, p. 1306-1313, D0I:10.1016/j.cad.2013.06.002.
[23]	Soori, M., Arezoo, B. (2021). Virtual Machining Systems for CNC Milling and Turning Machine Tools: A Review. International Journal of Engineering and Future Technology, vol. 18, no. 1, p. 56-104.
[24]	Soori, M., Arezoo, B., Habibi, M. (2016). Tool deflection error of three-axis computer numerical control milling machines, monitoring and minimizing by a virtual machining system. Journal of Manufacturing Science and Engineering, vol. 138, no. 8, art. ID. 081005, D0I:10.1115/1.4032393.
[25]	Soori, M., Asamel, M., Solyali, D. (2020). Recent development in friction stir welding process: A review. SAE International Journal of Materials and Manufacturing, vol. 14, no. 1, p. 6380, D0I:10.4271/05-14-01-0006.
[26]	Nasir, T., Kalaf, O., Asmael, M. (2021). Effect of rotational speed, and dwell time on the mechanical properties and microstructure of dissimilar AA5754 and AA7075-T651 aluminum sheet alloys by friction stir spot welding. Materials Science. Early Access, D0I: 10.5755/j02.ms.26860.
[27]	Zhang, X., Zhang, J., Pang, B., Zhao, W. (2016). An accurate prediction method of cutting forces in 5-axis flank milling of sculptured surface. International Journal of Machine Tools and Manufacture, vol. 104, p. 26-36, D0I:10.1016/j. ijmachtools.2015.12.003.
[28]	Mia, M., Dhar, N.R., (2016). Response surface and neural network based predictive models of cutting temperature in hard turning. Journal of Advanced Research, vol. 7, no. 6, p. 1035-1044, D0I:10.1016/j.jare.2016.05.004.
[29]	Shaw, M.C. (2005). Metal Cutting Principles. 2"d ed., Oxford University Press, Oxford.
[30]	Santos, M.C., Machado, A.R., Barrozo, M.A. (2018). Temperature in machining of aluminum alloys. in: temperature sensing. IntechOpen, D0I:10.5772/intechopen.75943.
[31]	Wan, M., Zhang, W.H., Qin, G., Wang, Z.P. (2008). Strategies for error prediction and error control in peripheral milling of thin-walled workpiece. International Journal of Machine Tools and Manufacture, vol. 48, no. 12-13, p. 1366-1374, D0I:10.1016/j.ijmachtools.2008.05.005.
[32]	Palanisamy, P., Rajendran, I., Shanmugasundaram, S. (2007). Optimization of machining parameters using genetic algorithm and experimental validation for end-milling operations. The International Journal of Advanced Manufacturing Technology, vol. 32, p. 644-655, D0I:10.1007/s00170-005-0384-3.
[33]	Tolouei-Rad, M., Bidhendi, I. (1997). On the optimization of machining parameters for milling operations. International Journal of Machine Tools and Manufacture, vol. 37, no. 1, p. 1-16, D0I:10.1016/S0890-6955(96)00044-2.
244
Soori, M. - Asmael, M.
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 245-255 © 2021 Journal of Mechanical Engineering. All rights reserved.
D0I:10.5545/sv-jme.2020.7078	Original Scientific Paper
Received for review: 2020-12-22
Received revised form: 2021-04-02
Accepted for publication: 2021-05-03
Ultrasonic Scattering Attenuation in Nodular Cast Iron: Experimental and Simulation Studies
Mario Santos* - Jaime Santos
University of Coimbra, CEMMPRE, Department of Electrical and Computer Engineering, Portugal
This work evaluates the ultrasonic scattering attenuation of structures with complex scatterer distributions via experimental and simulation studies. The proposed approach uses experimental attenuation knowledge to infer the scatterer size and its concentration in the studied structures, which are important for the effective construction of simulated models. The MATLAB k-Wave toolbox has been used to implement the simulator. Several cast-iron samples have been used to demonstrate the importance of simulation in the characterization of such structures. First, the scattering attenuation was evaluated using the Truell and Papadakis models, and then the results were compared with experimental ones. Emphasis was given to the Papadakis approach because it takes into account the scatterer size distribution. It is demonstrated that both analytical models provide results that are far from the experimental ones. The developed simulator for the studied samples led to a predictive model, in which the attenuation was proportional to the fifth power of the scatterer size, and the corresponding formulation is close to the one proposed by the analytical models. Keywords: modelling, anisotropy, pulse-echo, simulation, ultrasonic attenuation
Highlights
•	Experimental attenuation in cast iron samples was carried out.
•	Scattering attenuation theoretical models do not apply to complex nodular cast-iron structures.
•	Simulation models as a strategy to predict the experimental performance of cast iron.
•	A k-Wave simplified simulation model is used to characterize complex structures.
0 INTRODUCTION
There are many applications for nodular cast iron due to its castability, high thermal conductivity, and good mechanical properties, specifically tensile strength and ductility. The mechanical properties of a metal greatly depend on the microstructure; in the case of nodular cast iron, which is produced by adding, shortly before solidification, a small amount (lower than 0.04 %) of substances such as magnesium or cerium are present. These substances give rise to the growth of nodular graphite, whose shape and distribution are of fundamental importance in the behaviour of the metal [1] to [3]. Thus, the non-destructive evaluation of such structures is very important for the identification of the nodularity and matrix phases.
Ultrasonic characterization offers great advantages when compared with destructive metallographic methods. The interaction of ultrasound waves with the material microstructure can be evaluated, measuring the acoustic parameters, including velocity and attenuation. Two ultrasound attenuation mechanisms are generally identified: absorption and scattering. Absorption is related to thermal conduction loss, hysteresis, and a viscous loss mechanism [4]. Scattering is due to heterogeneities such as grain boundaries, voids, inclusions, second-phase particles or porosity [5] to [9]. This attenuation
mechanism is commonly accepted as the most important in heterogeneous materials, such as the case of cast iron [10] to [13]. It is important to take into account the fact that the scattering effects of the matrix grains are too small when compared to the nodular scattering effects and can be ignored [6], [7], [14] to [16].
Several authors have extensively studied scattering attenuation. At the beginning of the last century, Rayleigh presented a scattering formula [17], later adapted by Mason and McSkimin [18] and [19]; the case of polycrystalline aluminium. Huntington [20] used a stochastic theory to explain the scattering effects in polycrystalline structures. Lifshitz and Parkhomovskii [21] proposed a theory that considers the mode conversion at the grain boundaries. Moreover, a great contribution was made by Papadakis [5], [6], [22] to [25] with several published works related to that topic. The author classified the scattering in three classical regimes depending on the relation between the grain size (D) and the wavelength (X): (1) Rayleigh regime (for X>>D), where the attenuation is proportional to the fourth power of frequency; (2) stochastic regime (for XsD), where the attenuation is proportional to the square of frequency; (3) and geometrical regime (for X<<D), where the attenuation is frequency independent. For each regime, the microstructure is assumed to be
*Corr. Author's Address: DEEC, Polo II, Pinhal de Marrocos, 3030-290 Coimbra, Portugal, marioj@deec.uc.pt
245
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 245-255
composed of spherical grains of the same size filling the medium. As the Mason and McSkimin models [18] are very complex, due to multiple scattering, their theoretical treatment is a laborious task. That led other authors to restrict the scattering analysis to a limited frequency range [5], [24] to [27]. Later, Hirsekorn [28] and [29], Stanke and Kino [30], and Weaver [31] developed general solutions valid for all grain size to wavelength ratios and for cubic symmetry polycrystalline materials. However, other researchers demonstrated that those general solutions failed in more complex microstructures, such as the commercial aluminium alloys [32] and two-phase sintered powders [33]. More recent works deal with advances related to the interaction of ultrasound waves with polycrystalline materials. Yang et al. [34] studied the shape effect of elongated grains on attenuation. Arguelles and Turner [35] evaluated the errors resulting from neglecting the grain size distribution and using an average grain size. Ryzy et al. [36] used a semi-analytical attenuation model that considers the grain morphology and incorporates an exact spatial two-point correlation function. The authors concluded that the grain shape has a strong effect on attenuation in the Rayleigh-stochastic transition region. Rokhlin et al. [37] presented a model that includes second-order multiple scattering, applicable to all frequency ranges, providing small relative errors on both longitudinal and transverse attenuations for low anisotropy. Sha
[38]	extended the second-order attenuation (SOA) model for elastic waves in texture-free materials to textured polycrystals with ellipsoidal grains of arbitrary crystal symmetry. The predicted attenuation results of this work agree well with the literature on a textured stainless steel polycrystal.
Numerical or grid-based methods are also powerful tools in the analysis of the propagation in scattering media. The improvement of computational resources makes possible the implementation of those methods to study increasingly complex interactions
[39].	Recent works used finite element methods (FEM) to accurately mimic wave propagation through a rather complex material [40] to [44]. Norouzian, et al. [45] use DREAM.3D tool to construct material volumes with lognormal grain-size distributions. The results show that the correlation between attenuation and distribution width can be modelled with a power law, and the frequency dependence of attenuation has been shown to be strongly depend on the distribution width.
The main advantage of these methods, when compared to analytical approaches, is related to the fact that they allow microscopic scale analysis instead
of an averaged effective medium. No assumptions have to be made concerning grain statistics, anisotropy degree, or multiple scattering events [36]. Limitations related to frequency range, model dimensions and grid discretization tend to be strongly reduced, in the future, with the increase of computational resources.
One technique widely established for the generation of polycrystalline material morphologies is the Voronoi tessellation [40] to [50] for 2-D or 3-D dimensions. Details about that approach can be found in [50]. The major drawback of elastic wave models based on low-order finite difference or finite element schemes is the large number of grid points per wavelength required to avoid numerical dispersion. The k-Wave toolbox for Matlab can be used as an alternative to those approaches. It uses a Fourier domain pseudospectral method for the simulation and reconstruction of photoacoustic wave fields in a faster way. Also, it uses less memory and is user friendly [51] and [52]. Fewer spatial and temporal grid points are needed for accurate simulations [51]. Recently many authors have used the k-Wave toolbox in applications related to attenuation in ultrasonic computed tomography [53], guided waves in layered structures [54] and [55] time domain, power law attenuation in breast and liver tissues [56], one-sided ultrasonic nondestructive evaluation [57], high intensity focused ultrasound [58], nonlinear ultrasound propagation in absorbing media [59], ultrasonic transducers field modelling [60], and 3-D ultrasound imaging [61].
In the present work, the authors propose a simulation model to evaluate the scattering attenuation in structures with complex scatterer distributions. It uses the experimental attenuation knowledge in nodular cast iron to infer the scatterer size and its power dependence, as well as the concentration that best mimics the structure behaviour. The k-Wave toolbox is used to construct the simulator. Additionally, the results provided by the theoretical and simulation models are presented and a discussion is also made about the feasibility of the Truell et al. [62] and Papadakis [22] models for attenuation evaluation in nodular cast iron.
1 ULTRASONIC SCATTERING ATTENUATION 1.1 Theory
For the Rayleigh domain, characterized by a wavelength much higher than the scatterer size, Ying and Truell [27] presented a model for the total scattering cross-section, assuming the scatterers are spherical, solid, elastic, isotropic, have uniform size
246
Santos, M. - Santos, J.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 245-255
and are embedded in a solid elastic matrix. Based on those assumptions, Truell et al. [62] presented an expression for the ultrasonic scattering attenuation evaluation of longitudinal waves:
2n
a,
9 vi
(1)
where N is the number of scatterers per unit volume, m is the angular frequency, vl represents the longitudinal propagation velocity, r is the scatterer radius, and g denotes a factor related to the elastic properties of the medium. The original g expression in [27] can be changed to encompass the longitudinal (vl) and transversal (vt) velocities:
g =

v v'i y
V
-1
v
V V 2
- 4
+ 4
1
+
3
+40
f	f Vl '	3 ^	f \ V,iM2 1	2
1 + 2				+
v	v vti y	y	V V 2 Mi y	
2 + 3
n5 >
v vti y
Ml Mi
-1
f \2 vii
+ 2
V v
Ml
Mi
2
vii
- 4
(2)
y y
where f is the shear modulus given by vt2 p ; p is the density, and the indexes 1 and 2 are related to the surrounding medium and scatterer, respectively.
Based on a similar work developed by Bathia and Moore [63], Papadakis [6] and [22] introduced a correction factor for the attenuation considering the effect of the grain size distribution. The author also showed a way of obtaining the correction factor in solids using two-dimension micrographs. The scatterer volume correction factor (T) is:
T=zL
3 3
3 r
(3)
where rn6 and rn3 are the average scatterer radius values of the sixth and third power, respectively, obtained from micrograph samples. From Eq. (3), it is clear that the scattering effects of a single large
scatterer are much greater than the effect summation of many small scatterers having the same volume. Taking the correction factor into account, the scattering attenuation is given as [6]:
= NT
4n*f4
(f
1 Ap^_ AL2 8 ALAp 32 Ap
HL2
2 A
3 Pl2
Ll
Mi
'2 ( 2 Ap2
3 Pi2
3 L 8
H--
15
15 Ll2 Ap
\
Li2
(4)
where f is the frequency, Lj and pj are the longitudinal and the shear modulus of surrounding medium, respectively, and Ap, AL, Ap are the differences in density, longitudinal modulus, and shear modulus between scatterers and the surrounding medium, respectively. Although Eq. (4) uses a different approach for the attenuation calculation, the results are similar to those provided by the Truell model for uniform scatterer sizes.
1.2 Experimental Attenuation
The attenuation is the result of acoustic wave interactions with the propagation medium. Experimentally, the attenuation can be calculated by collecting the front surface reflection and two back wall ultrasound pulses from samples with parallel faces, as illustrated in Fig. 1. Following that approach, the authors studied different methodologies to calculate the ultrasonic attenuation in nodular cast iron samples [64], and concluded that Eq. (5) provided more reliable results, because it does not take into account the reflection coefficient between water and sample, which is frequency-dependent, as demonstrated experimentally [64]:
a = —ln 2L
f A AD ^
A 2 V Ai
- A2 As
(5)
In Eq. (5), L is the sample thickness, Aj andA2 are the first and second back wall signals from samples, As is the reflected signal on the sample front face, and Dc is the diffraction correction coefficient for the sample path. The attenuation (a) unit is [Np/m], where Np is the symbol of neper, which is a logarithmic unit for ratios of measurements of physical quantities. The neper and dB are related by the following relationship: 1 Np = 8.686 dB
2
3
2
v
2
v vti y
Ultrasonic Scattering Attenuation in Nodular Cast Iron: Experimental and Simulation Studies
247
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 245-255
is A, .h
Ao Wilier						
Sample			S N		O1	
Water
Fig. 1. Setup for attenuation measurement 1.3 k-Wave Simulation
The simulation of elastic wave propagation has many applications in ultrasonic non-destructive testing. The open-source k-Wave toolbox is an easy-to-use time-domain forward model based on a k-space pseudospectral time-domain solution to couple first-order acoustic equations for homogeneous and heterogeneous media in one, two, and three dimensions [51]. That tool makes possible the modelling of arbitrary sources, detecting surfaces with directional elements, and wave propagation that can account for nonlinearity, acoustic heterogeneities and power-law absorption. It is also possible to use optional input parameters to adjust the visualization and performance, including making wave propagation movies for presentations and running simulations on graphics processing units (GPU). The available geometry creation functions allow both Cartesian and grid-based geometries, such as circles, arcs, disks, spheres, shells or balls. As an example, Fig. 2 shows a three-dimensional model for a circular transducer (top) radiating into a sample with randomly distributed scatterers. The main drawback of k-Wave in ultrasonic non-destructive testing is related to the existence of large contrasts between media. Fourier-based methods such as k-Wave are effective when everything is smooth, and the contrasts are not too big (the code was originally designed for modelling ultrasound in biological tissue, where the contrasts are low). For large contrasts, oscillatory errors can be accumulated, because the steep edges cannot be represented very well by a small number of Fourier components. Some caution must be taken in such situations. Smoothing the density and propagation velocity or keeping the time step small will help to reduce that effect.
y [voxels]	x [voxels]
Fig. 2. Three-dimensional simulation model
2 RESULTS AND DISCUSSION
2.1 Samples Granulometry Analysis
Six machined cast iron samples with parallel surfaces within 15 ^m and 13.1 ± 0.05 mm in thickness were used in this study (Fig. 3). The samples were prepared using standard metallographic methods. First, they were mechanically polished using metallographic carbon silicate sandpaper with decreasing granulometry (P180, P1000, and P2500) and a final polishing by a 3-^m diamond suspension. The microstructure of the processed zones was examined using optical microscopy, with a 200* magnification.
The typical microstructure of such samples is shown in Fig. 4, where the precipitated graphite nodules (in black) embedded in a ferrite/perlite matrix are easily identified, as are the grain boundaries. It was observed in all samples that the nodules present non-uniform shapes and have a large wide range of sizes that can vary from less than one micrometre to several tens of micrometres (the larger measured radius was 35 |im).
Fig. 3. Cast iron samples
248
Santos, M. - Santos, J.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 245-255
computer for further processing. The transducer is moved using a computer-controlled micro-positioning system.
Fig. 4. Optical micrograph of a nodular cast iron sample
Ten micrograph images were taken from each sample, corresponding to an analysis section of 10.32 mm2. The ImageJ free package software was used to determine the nodules' size [65], which were considered to be spherical. The nodule size distribution is presented in Fig. 5. Most of the nodules have a radius less than 10 ^m. The spatial variability of the nodules' size in different regions of each sample is low. The standard deviation of the nodule size divided by its average value obtained for all images of each sample was about 8.9 %. Also, for each sample, the total number of nodules was used to obtain the lognormal distribution parameters. The mean and standard deviation are presented in Table 1.
2.2 Experimental Attenuation Evaluation Setup and Results
The experimental immersion setup for the attenuation evaluation is presented in Fig. 6. A pulserireceiver is used to excite a broadband 8 MHz central frequency transducer and to collect, amplify, and filter the reflected signals. Then, the acquired signals are displayed in an oscilloscope and transferred to a
Nodule r;idius [/ym]
Fig. 5. Samples' nodule size distribution
Fig. 6. Experimental setup for attenuation evaluation
All samples characterized in section 2.1 were used in the study. Four acquisitions were accomplished in different regions of each sample to take into account their structural variation. Then, using Eq. (5) an average attenuation of 11.96 Np/m was obtained,
Table 1. Lognormal fitting parameters for each sample
Sample	1	2	3	4	5	6
Mean ^m]	0.2335	0.1178	0.4310	0.1007	0.0671	0.2420
Standard deviation [^m]	0.9492	0.8026	0.7862	0.8751	0.8545	0.9161
Table 2. Acoustic properties of cast iron components
	Longitudinal velocity	Transversal velocity	Longitudinal modulus	Shear modulus	Density
	[m/s]	[m/s]	[GPa]	[GPa]	[kg/m3]
Pearlite/Ferrite	5830	3090	269	75.6	7920
Graphite	4210	2030	38	8.9	2170
Ultrasonic Scattering Attenuation in Nodular Cast Iron: Experimental and Simulation Studies	249
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 245-255
with a standard deviation (SD) of 1.22 Np/m and root mean square error (RMSE) of 1.10 Np/m for the 24 measurements.
2.3 Attenuation Evaluation Using Analytical Models
In order to evaluate the scattering attenuation in the cast iron samples, important parameters as longitudinal and shear velocities, longitudinal and shear moduli, and density related to the matrix (pearlite/ferrite) and scatterers (graphite), must be known. They are represented in Table 2 [24].
The scatterer concentration knowledge is also required to estimate the attenuation. Thus, the two-dimensional scatterer concentration is obtained relating the total amount of nodules (nT) with the area of micrograph (AT), as follows:

(6)
Then, the scatterer concentration (N can be calculated using the following equation [66]:
N =
i
^ i
6/
{^nnA )
(7)
g J
where fg is the graphite nodule fraction that is about 3.33 % for the analysed images, and aN = 1.25, which is a parameter related to the width of the nodule size distribution.
To infer the real scatterer effect in the attenuation, the scatterers that contribute to the attenuation should be determined. To do so both the Truell and Papadakis models will be used and their performance analysed.
Due to the wide nodular distribution observed in the samples, the Truell model given by Eq. (1) cannot be used to calculate the scattering attenuation, because it uses the same size scatterers. Thus, in order to encompass all scatterer sizes found in the samples, the authors present a new equation for the attenuation, based on Eq. (1), which takes into account the summation of the different classes of scattering nodules:
-A 2n m 6 asrm =	g — Niri .
i=i 9 vi
(8)
In Eq. (8), n is the total number of scatterer classes, r, and N are the average radius and scatterer concentration of each class, respectively. Thirty-five scatterer classes with 1 ^m step were considered for the attenuation calculation. The scatterer concentration of each class is a fraction of the one obtained by Eq.
(7). From Eq. (8), resulted an attenuation of 2.67 Np/m.
The different attenuation values provided by the Truell model and experimental approach could be due to the fact that the Truell model makes use of spherical, solid, elastic, and isotropic scatterers. In addition, the referred model assumes that there is a sharp variation of the properties between the scatterers and the surrounding media. However, the nodular cast iron structure is rather more complex, especially its matrix, which is usually formed by ferrite around the nodules, with a homogeneous structure, and by pearlite in other regions with a lath-type structure with some degree of heterogeneity [24]. Therefore, the boundary effect between scatterers and matrix is certainly more complex than the one presented in the model.
The Papadakis model takes into account the scatterer size distribution, which is an improved version of the Truell theory. Generally, in practical materials, such as cast iron, the scatterers are not all the same size and follow a distribution that can be evaluated. The correction factor (see Eq. (3)), requires
the knowledge of r6 and r3:
=zZL,
ti n
n 3
(9)
(10)
Using the experimental scatterer distribution
(see Fig. 5), resulted rn6 = 3.5 x 105 ^m6 and
rn3 = 54.73 ^m3 for the whole sample area analysed, and T = 2.68 x 10-14 m3.
Then, using Eq. (4), where the right side is essentially composed by the physical parameters of matrix and scatterers presented in Table 1, resulted asP = 370 Np/m, which is very large compared to the experimental attenuation. In addition to the explanation made for the Truell model deviation in relation to the experimental attenuation, which also applies to the Papadakis model, other important factors can contribute to the observed discrepancy between theory and experiment. Perhaps, the most important one deals with the scatterer distribution obtained from the micrograph images. In the Papadakis model, all scatterer sizes and shapes were taken into account for the attenuation calculation that originates a very high scatterer concentration value. However, the authors concluded the attenuation decreases sharply if the lower size scatterers are discarded, as illustrated in Fig. 7. That tendency is expected because the
6
6
r
r
250
Santos, M. - Santos, J.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 245-255
scatterer concentration (N) falls more quickly than the correction factor (T). For instance, if all scatterers lower than 7 |im are discarded, the attenuation is 15 Np/m, which is very close to the experimental one.
Finally, the Papadakis model also presents limitations in predicting the scattering attenuation in the cast iron samples.
Fig. 7. Attenuation behaviour versus scatterer size 2.4 Simulation Model
The simulation of the cast iron nodular structure shown in Fig. 4 is difficult due to the non-uniformity of scatterers in size and shape. Here, the authors propose a simplified model assuming all scatterers with the same size and known concentration that leads to the experimental attenuation value. For that purpose, it is assumed that the attenuation is directly proportional to the scatterer concentration (see Eq. (11)), which is in accordance with the theoretical models [21] and [60] and is a common-sense assumption. The equivalent scatterer size (req) and its exponent dependence (n) will be determined based on the inverse problem, using the experimental attenuation. The simplified attenuation expression based on the previous assumptions is as follows,
= CN(req J
(11)
where C and n are constants to be determined. Eq. (11) follows the Rayleigh scattering model, and only the geometric parameters are of interest, such as the scatterer concentration, equivalent scatterer size, and its exponent dependence. The frequency dependence and the matrix and scatterer properties are included in the constant C, because the study related to the attenuation variation with these parameters is outside the scope of this work.
Using k-Wave, a simulation model with dimensions x = 6 mm, y = 6 mm and z = 40 mm, similar to the one illustrated in Fig. 2, was defined by a computational grid for simulation purposes. For the scattering attenuation evaluation, an ultrasonic probe is located at the bottom of the model (z = 0), working as a transmitter and another at z = 30 mm working as a receiver. Two signals are collected at the receiver: a reference signal a0 (model without scatterers) and a signal as considering the model filled up with randomly distributed scatterers. The scattering attenuation is then given by:
a,,
= — ln
d
( „
(12)
where d = 30 mm is the distance between the probes. The two probes are 6 mm in diameter, like the one used in the experimental measurements (see section 2.2). The probe excitation was made by a tone burst, as shown in Fig. 8, to mimic the experimental transducer response illustrated in Fig. 9.
The first set of simulations were carried out for scatterers of req = 60 ^m, corresponding to a grid point spacing of 120 ^m (cuboid edge). The scatterer concentration was varied from 0 mm-3 to 100 mm-3 with a step of 10 mm-3, and ten signals as were collected for each concentration value. The corresponding attenuation was then calculated using Eq. (12), and the results are shown in Fig. 10. An average attenuation value for each concentration was considered. As expected, the attenuation increases as N also increases. That behaviour is more pronounced for lower concentrations. In the boxplot of the simulated data (Fig. 10) the central mark is the median, the edges of the box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points (not considering the outliers, which are plotted individually as crosses).
	0JS
(LI	D.fi
	
	li.il
c.	
£	0J
	
	tl
rj	
	
■ —	-11,2
13	
§	II. 4
Q	
	-flji
	■(IK
1 A	I
1	V
15	I	L5
Time |//s|
Fig. 8. Simulated tone burst
0
a
V V /
Ultrasonic Scattering Attenuation in Nodular Cast Iron: Experimental and Simulation Studies
251
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 245-255
Fig. 9. Experimental transducer response
For the simulated scatterer size, the concentration value that leads to the experimental attenuation calculated in section 2.2 (aexp = 11.96 Np/m), is about 3.2 mm-3. This result is more clearly observed in Fig. 11, which corresponds to an expansion of Fig. 10.
Fig. 10. Simulated attenuation for req = 60 ym versus scatterer concentration
N [mm" ]
Fig. 11. Zoomed version of Fig. 10 for lower concentrations
In order to establish a correlation between the scatterer concentration (AO and size (req), additional attenuation simulations were carried out using different scatterer sizes, resulting in the concentrations (presented in Table 2) that lead to the experimental attenuation. A power fitting applied to the values in Table 2 gave rise to Eq. (13), whose behaviour is illustrated in Fig. 12, where a remarkable goodness of fit with R2 = 0.99 is observed:
N = 4.527 x 10"
'h )"
(13)
Table 2. Scatterers dimensions and concentration for aexp = 11.96 Np/m
req M	40	50	60	70	80
N [mm-3]	25.5	8.0	3.2	1.4	0.68
Fig. 12. Power fitting of scatterer concentration N
as function of re
eq
Solving the proposed simplified model given by Eq. (11), in order to extract the scatterer concentration:
N =
11.96 C
k F
(14)
and comparing the Eqs. (13) and (14), results C=2.65 x 1013 and n = 5.173, which allows writing Eq. (11) as:
a = 2.65 x1013 N (req
(15)
where the attenuation [Np/m] is a function of the fifth power of the scatterer size, which is not far from the trend of Eq. (1), as expected. It should be pointed out that Eq. (15) resulted from the studied cast iron samples illustrated in Fig. 4, whose attenuation was experimentally calculated by Eq. (5). Thus, based on the previous knowledge of experimental
252
Santos, M. - Santos, J.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 245-255
attenuation, simulation models can be created varying the parameters N and req, according to Eq. (14). This allows mimicking structures with complex scatterer distributions in terms of uniform scatterer concentration and size.
The presented simulation model considers the geometrical factors (size and concentration of scatterers) and can be used to mimic the microstructure of cast iron for a specific excitation frequency. The scatterer equivalent size can be defined a priori and the concentration calculated by Eq. (14) to establish the simulation model.
3 CONCLUSIONS
This work aimed to evaluate the ultrasonic scattering attenuation of structures with complex scatterer distributions. The idea was to predict the attenuation behaviour for such complex structures by developing a simulation model based on the experimental attenuation results. That was motivated by the limitations presented by the Truell and Papadakis theoretical models that provide ultrasonic scattering attenuation results quite inconsistently when compared with the experimental one. That conclusion was observed for nodular cast-iron structures, and the authors believe that the same results are expected for materials having inhomogeneous scatterer sizes and concentrations. The k-Wave simulation model developed in this work allows an easy and customized implementation, where the scatterer size and concentration can be varied to mimic the experimental ultrasonic scattering attenuation measured in inhomogeneous structures.
4 ACKNOWLEDGEMENTS
This research is sponsored by FEDER funds through the program COMPETE - Programa Operacional Factores de Competitividade - and by national funds through FCT - Fundaçao para a Ciência e a Tecnologia - under the project UIDB/00285/2020.
5 REFERENCES
[1]	Gur, C., Aydinmakina, B. (2010). Microstructural Characterisation of Ductile Irons by Measuring Velocity and Apparent Attenuation of Ultrasonic Waves. Insight, vol. 43, no. 11, p. 731-734.
[2]	Prakash, P., Mytri, V.D., Hiremath, P.S. (2011). Digital microstructure analysis system for testing and quantifying the ductile cast iron. International Journal of Computer Applications, vol. 19, no. 3, p. 22-27, DOI:10.5120/2342-3062.
[3]	Cheng, T.-W., Lui, T.-S., Chen, L.-H. (2012). Microstructural features and erosion wear resistance of friction stir surface hardened spheroidal graphite cast iron. Materials Transactions, vol. 53, no. 1, p. 167-172, D0l:10.2320/ matertrans.M2011224.
[4]	Ristic, V.M. (1983). Principles of Acoustic Devices. John Wiley & Sons, New York.
[5]	Papadakis, E.P. (1965). Ultrasonic attenuation caused by scattering in polycrystalline metals. Journal of Acoustical Society of America, vol. 37, no. 4, p. 711-717, D0I:10.1121/1.1909401.
[6]	Papadakis, E.P. (1981). Ultrasonic attenuation caused by Rayleigh scattering by graphite nodules in nodular cast iron. Journal of Acoustical Society of America, vol. 70, no. 3, p. 782787, D0I:10.1121/1.386916.
[7]	Kruger, S.E., Rebello, J.M.A., Charlier, J. (2004). Broadband ultrasonic backscattering applied to nondestructive characterization of materials. IEEE Transaction on Ultrasonic and Frequency Control, vol. 51, no. 7, p. 832-838, D0I:10.1109/TUFFC.2004.1320742.
[8]	Palanichamy, P., Joseph, A., Jayakumar, T., Baldev, R. (1995). Ultrasonic velocity measurements for estimation of grain size in austenitic stainless steel. NDT&E International, vol. 28, no. 3, p. 179-185, D0I:10.1016/0963-8695(95)00011-L.
[9]	Zeng, F., Agnew, S.R., Raeisinia, B., Myneni, G.R. (2010). Ultrasonic attenuation due to grain boundary scattering in pure niobium. Journal of Nondestructive Evaluation, vol. 29, no. 2, p. 93-103, D0I:10.1007/s10921-010-0068-2.
[10]	Thompson, R. (2002). Elastic-wave propagation in random polycrystals: Fundamentals and application to nondestructive evaluation. Fink, M., Kuperman, W.A., Montagner, J.P., Tourin, A. (eds.). Imaging of Complex Media with Acoustic and Seismic Waves. Topics in Applied Physics, Springer-Verlag, Berlin, Germany, vol. 84, p. 233-257, D0I:10.1007/3-540-44680-X_9.
[11]	Thompson, R.B., Margetan, F.J., Haldipur, P., Yu, L., Li, A., Panetta, P., Wasan, H. (2008). Scattering of elastic waves in simple and complex polycrystals. Wave Motion, vol. 45, no. 5, p. 655-674, D0I:10.1016/j.wavemoti.2007.09.008.
[12]	Pawlicki, A.D., O'Brien, W.Jr. (2013). Method for estimating total attenuation from a spatial map of attenuation slope for quantitative ultrasound imaging. Ultrasonic Imaging, vol. 35, no. 2, p. 162-172, D0I:10.1177/0161734613478695.
[13]	Lobkis, O.I, Yang, L., Li, J., Rokhlin, S.I. (2012). Ultrasonic backscattering in polycrystals with elongated single phase and duplex microstructure. Ultrasonics, vol. 52, no. 6, p. 694-705, D0I:10.1016/j.ultras.2011.12.002.
[14]	Lee, S.C., Suen, J.M. (1989). Ultrasonic nondestructive evaluation of matrix structures and nodularity in cast irons. Metallurgical Transactions A, no. 20, p. 2399-2407, D0I:10.1007/BF02666675.
[15]	Liu, X., Takamori, S., Osawa, Y. (2007). Effect of matrix structure on ultrasonic attenuation of ductile cast iron. Journal of Materials Science, vol. 42, p. 179-184, D0I:10.1007/ s10853-006-0845-6.
[16]	Kruger, S.E., Rebello, J.M.A,, Charlier, J. (2002). Measuring cast iron graphite size by ultrasonic attenuation. AIP Conference Proceedings, vol. 21, p. 1531-1537, D0I:10.1063/1.1472975.
[17]	Rayleigh, L. (1929). Theory of Sound. McMillan, New York.
Ultrasonic Scattering Attenuation in Nodular Cast Iron: Experimental and Simulation Studies
253
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 245-255
[18]	Mason, W.P, McSkimin, H.J. (1947). Attenuation and scattering of high frequency sound waves in metals and glasses. Journal of Acoustical Society of America, vol. 19, no. 3, p. 464-473, D0l:10.1121/1.1916504.
[19]	Mason, W.P., McSkimin, H.J. (1948). Energy losses of sound waves in metals due to scattering and diffusion. Journal of Applied Physics, no. 19, p. 940-946, D0I:10.1063/1.1697900.
[20]	Roth, W. (1948). Scattering of ultrasonic radiation in polycrystalline metals. Journal of Applied Physics, no. 19, no. 10, p. 901-910, D0I:10.1063/1.1697896.
[21]	Lifshitz. I., Parkhomovskii, G. (1950) Theory of propagation of ultrasonic waves in polycrystals. Journal of Experimental and Theoretical Physics, vol. 20, p.175-182.
[22]	Papadakis, E.P. (1961). Grain size distribution in metals and its influence on ultrasonic attenuation measurement. Journal of Acoustical Society of America, vol. 33, no. 11, p. 16161621, D0I:10.1121/1.1908517.
[23]	Papadakis, E.P. (1964). From micrograph to grain-size distribution with ultrasonic applications. Journal of Applied Physics, vol. 35, no. 5, p. 1586-1594, D0I:10.1063/1.1713671.
[24]	Papadakis, E.P. (1965). Revised grain-scattering formulas and tables. Journal of Acoustical Society of America, vol. 37, p. 703-710, D0I:10.1121/1.1909399.
[25]	Papadakis, E.P. (1963). Rayleigh and Stochastic Scattering of Ultrasonic Waves in Steel. Journal of Applied Physics, vol. 34, no. 2, p. 265-269, D0I:10.1063/1.1702596.
[26]	Waterman, P., Truell, R. (1961). Multiple scattering of waves. Journal of Mathematical Physics, vol. 2, no. 4, p. 512-537, D0I:10.1063/1.1703737.
[27]	Ying, C.F., Truell, R. (1956). Scattering of a plane longitudinal wave by a spherical obstacle in an isotropically elastic solid. Journal of Applied Physics, vol. 27, no. 9, p. 1086-1097, D0I:10.1063/1.1722545.
[28]	Hirsekorn, S. (1982). The scattering of ultrasonic waves by polycrystals. Journal of Acoustical Society of America, vol. 72, no. 4, p. 1021-1031, D0I:10.1121/1.389206.
[29]	Hirsekorn, S. (1983). The scattering of ultrasonic waves by polycrystals. II. Shear waves. Journal of Acoustical Society of America, vol. 73, no.4, p. 1160-1163, D0I:10.1121/1.389206.
[30]	Stanke, F.E, Kino, G.S. (1984). A unified theory for elastic wave propagation in polycrystalline materials. Journal of Acoustical Society of America, vol. 75, no. 3, p. 665-681, D0I:10.1121/1.390577.
[31]	Weaver, R. (1990). Diffusivity of ultrasound in polycrystals. Journal of the Mechanics and Physics of Solids, vol. 38, no. 1, p. 55-86, D0I:10.1016/0022-5096(90)90021-U.
[32]	Man, C., Paroni, R., Xiang, Y., Kenik, E.A. (2006). On the geometric autocorrelation function of polycrystalline materials. Journal of Computational and Applied Mathematics, vol. 160, no. 1-2, p. 200-210, D0I:10.1016/j.cam.2005.01.044.
[33]	Liu, D., Turner, JA. (2008). Influence of spatial correlation function on attenuation of ultrasonic waves in two-phase materials. Journal of Acoustical Society of America, vol. 123, no. 5, p. 2570-2376, D0I:10.1121/1.2896757.
[34]	Yang, L., Lobkis, O.I., Rokhlin, S.I. (2011). Shape effect of elongated grains on ultrasonic attenuation in polycrystalline materials. Ultrasonics, vol. 51, no. 6, p. 697-708, D0I:10.1016/j.ultras.2011.02.002.
[35]	Arguelles, A.P., Turner, J.A. (2017). Ultrasonic attenuation of polycrystalline materials with a distribution of grain sizes. Journal of the Acoustical Society of America, vol. 141, no. 6, p. 4347-4353, D0l:10.1121/1.4984290.
[36]	Ryzy, M., Grabec, T., Sedlak, P., Veres, I.A. (2018). Influence of grain morphology on ultrasonic wave attenuation in polycrystalline media with statistically equiaxed grains. Journal of the Acoustical Society of America, vol. 143, no. 1, p. 219-229, D0I:10.1121/1.5020785.
[37]	Rokhlin, S.I., Li, J., Sha, G. (2015). Far-field scattering model for wave propagation in random media. Journal of the Acoustical Society of America, vol. 137, no. 5, p. 2655-2669, D0I:10.1121/1.4919333.
[38]	Sha, G. (2020). Attenuation and phase velocity of elastic wave in textured polycrystals with ellipsoidal grains of arbitrary crystal symmetry. Acoustics, vol. 2, no. 1, p. 51-72, D0I:10.3390/acoustics2010005.
[39]	Van Pamel, A., Sha, A., Rokhlin, S.I., Lowe, M.J.S. (2016). Finite-element modelling of elastic wave propagation and scattering within heterogeneous media. Proceedings of the Royal Society A, vol. 473, no. 2197, p. 1-21, D0I:10.1098/rspa.2016.0738.
[40]	Ghoshal, G., Turner, J.A. (2009). Numerical model of longitudinal wave scattering in polycrystals. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 56, no. 7, p. 1419-1428, D0I:10.1109/TUFFC.2009.1197.
[41]	Lan, B., Lowe, M., Dunne, F.P.E. (2014). Experimental and computational studies of ultrasound wave propagation in hexagonal close-packed polycrystals for texture detection. Acta Materialia, vol. 63, p. 107-122, D0I:10.1016/j. actamat.2013.10.012.
[42]	Van Pamel, A., Brett, C., Huthwaite, P., Lowe, M. (2015). Finite element modelling of elastic wave scattering within a polycrystalline material in two and three dimensions. Journal of the Acoustical Society of America, vol. 138, no. 4, p. 23262336, D0I:10.1121/1.4931445.
[43]	Shahjahan, S., Rupin, F., Aubry, A., Chassignole, B., Fouquet, T., Derode, A. (2014). Comparison between experimental and 2-D numerical studies of multiple scattering in inconel600® by means of array probes. Ultrasonics, vol. 54, no. 1, p 358367, D0I:10.1016/j.ultras.2013.06.012.
[44]	Chassignole, B., Duwig, V., Ploix, M.-A., Guy, P., El Guerjouma, R. (2009). Modelling the attenuation in the ATHENA finite elements code for the ultrasonic testing of austenitic stainless steel welds. Ultrasonics, vol. 49, no. 8, p. 653-658, D0I:10.1016/j.ultras.2009.04.001.
[45]	Norouzian, M., Islam S., Turner. J.A. (2020). Influence of microstructural grain-size distribution on ultrasonic scattering. Ultrasonics, vol. 102, p. 1-7, D0I:10.1016/j. ultras.2019.106032.
[46]	Voronoi, G. (1908). Nouvelles applications des paramètres continus à la théorie des formes quadratiques. Premier mémoire. Sur quelques propriétés des formes quadratiques positives parfaites. Journal fur die Reine und Angewandte Mathematik, no. 133, p. 198-287, D0I:10.1515/ crll.1908.133.97.
[47]	O'Rourke, J. (1988). Computational Geometry in C, Cambridge University Press, Cambridge.
254
Santos, M. - Santos, J.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 245-255
[48]	Zhang, P., Balint, D., Lin, J. (2011). An integrated scheme for crystal plasticity analysis: virtual grain structure generation. Computational Materials Science, vol. 50, no. 10, p. 28542864, D0I:10.1016/j.commatsci.2011.04.041.
[49]	Zhang, K., Holmedal, B., Hopperstad, O., Dumoulin, S., Gawad, J., Van Bael, A., Van Houtte, P. (2015). Multi-level modelling of mechanical anisotropy of commercial pure aluminium plate: crystal plasticity models, advanced yield functions and parameter identification. International Journal of Plasticity, vol. 66, p. 3-30, D0I:10.1016/j.ijplas.2014.02.003.
[50]	Quey, R., Dawson, P.R., Barbe, F. (2011). Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing. Computer Methods in Applied Mechanics and Engineering, vol. 200, no. 17-20, p. 1729-1745, D0I:10.1016/j.cma.2011.01.002.
[51]	Treeby, B.E., Cox, B.T. (2010). k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave-fields. Journal of Biomedical Optics, vol. 15, no. 2, p. 21314-121314-12, D0I:10.1117/1.3360308.
[52]	Treeby, B., Jaros, J., Rendell, A., Cox, B. (2012). Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method. Journal of the Acoustical Society of America, vol. 131, no. 6, p. 4324-4336, D0I:10.1121/1.4712021.
[53]	Pérez-Liva, M., Herraiz, J.L., Udías, J.M., Miller, E., Cox, B.T., Treeby, B.E. (2017). Time domain reconstruction of sound speed and attenuation in ultrasound computed tomography using full wave inversion. Journal of the Acoustical Society of America, vol. 141, no. 3, p. 1595-1604, D0I:10.1121/1.4976688.
[54]	Brill, T.M., Klieber, C. (2016). Reflection and mode-conversion of ultrasonic leaky Lamb waves at inaccessible discontinuities in layered structures. 2016 IEEE International Ultrasonics Symposium, D0I:10.1109/ULTSYM.2016.7728587.
[55]	Klieber, C., Brill, T.M. (2017). Mapping of ultrasonic Lamb-wave field in elastic layered structures using laser probes. Proceedings of Meetings on Acoustics, vol. 30, art. ID 065013, D0I:10.1121/2.0000603.
[56]	Zhao, X., McGougha, R.J. (2018). Time-domain analysis of power law attenuation in space-fractional wave equations. Journal of the Acoustical Society of America, vol. 144, no.1, p. 467-477, D0I:10.1121/1.5047670.
[57]	Almansourin, H., Venkatakrishnan, S., Bouman, C., Santos-Villalobos, H. (2019). Model-based iterative reconstruction
for one-sided ultrasonic nondestructive evaluation. IEEE Transactions on Computational Imaging, vol. 5, no. 1, p. 150164, DOI:10.1109/TCI.2018.2883218.
[58]	Matsui, K., Azuma, T., Fujiwara, K., Takeuchi, H., Itani, K., Wang, J., Iwahashi, T., Kobayashi, E. Sakuma, I. (2017). Improving high-intensity focused ultrasound beam imaging via a backscattering suppression algorithm. Japanese Journal of Applied Physics, vol. 56, no. 5, p. 057301-1:057301-12, D0l:10.7567/JJAP.56.057301.
[59]	Wang, K., Teoh, E., Jaros, J., Treeby, B. (2012). Modelling nonlinear ultrasound propagation in absorbing media using the k-Wave toolbox: experimental validation. 2012 IEEE International Ultrasonics Symposium, D0I:10.1109/ ULTSYM.2012.0130.
[60]	Acquaticci, F., Guarracino, J.F., Gwirc, S.N., Lew, S.E. (2019). A polydimethylsiloxane-based axicon lens for focused ultrasonic brain stimulation techniques. Acoustical Science and Technology, vol. 40, no. 2, p. 116-126, DOI:10.1250/ ast.40.116.
[61]	Cui, W., Qin, K. (2018). Fast 3-D ultrasonic imaging using time-domain synthetic aperture focusing techniques based on circular scan conversions. IEEE Transactions on Computational Imaging, vol. 4, no. 4, p. 632-639, DOI:10.1109/TCI.2018.2870303.
[62]	Truell, R., Elbaum, C., Chick, B.B. (1969). Ultrasonic Methods in Solid State Physics, Academic Press, New Your, D0I:10.1016/C2013-0-12565-2.
[63]	Bhatia, B.A. Moore, A.R. (1959). Scattering of High-Frequency Sound Waves in Polycrystalline Materials. Journal of Acoustical Society of America, vol. 31, p. 1140-1141, D0I:10.1121/1.1907843.
[64]	Santos, M., Santos, J. (2019). Evaluation of three different approaches for the ultrasound attenuation coefficient measurement in nodular cast iron. Materials Evaluation, vol. 77, no. 5, p. 605-612.
[65]	Schneider, C.A., Rasband, W.S., Eliceiri, K.W. (2012). NIH Image to ImageJ: 25 years of image analysis. Nature Methods, vol. 9, p. 671-675, D0I:10.1038/nmeth.2089.
[66]	Pedersen, K.M., Tiedje, N.S. (2008). Graphite nodule count and size distribution in thin-walled ductile cast iron. Materials Characterization, vol. 59, no. 8, p. 1111-1121, D0I:10.1016/j. matchar.2007.09.001.
Ultrasonic Scattering Attenuation in Nodular Cast Iron: Experimental and Simulation Studies
255
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271 © 2021 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2020.7094	Review Scientific Paper
Received for review: 2021-01-15 Received revised form: 2021-04-22 Accepted for publication: 2021-05-06
Operating Performance of External Non-Involute Spur and Helical Gears: A Review
Ivan Okorn* - Marko Nagode - Jernej Klemenc
University of Ljubljana, Faculty of Mechanical Engineering, Ljubljana, Slovenia
In practical use, most gears have an involute shape of tooth flanks. However, external involute gears have some drawbacks, such as unfavourable kinematic conditions at the beginning and end of meshing, a limited minimum number of teeth, and the highly loaded convex-convex (i.e., non-conformal) contact. Researchers have developed and analysed various non-involute forms of tooth flanks, but they have not been widely accepted. The main reasons are higher manufacturing costs and sensitivity to manufacturing and assembly errors. Analyses of non-involute forms of teeth are mostly theoretical (analytical and numerical), while there is a lack of experimental confirmations of theoretical assumptions. This paper reviews external non-involute shapes, their operating characteristics and possibilities of use compared to involute gears. Established criteria, such as Hertzian pressure, oil film thickness, bending stress at the root of the tooth, contact temperature, and gear noise, were used for assessment. The results of analytical studies and experimental research on S-gears are presented in more detail. S-gears have a higher surface durability and a lower heat load when compared to involute gears. The usability of non-involute gears is increasing with the development of new technologies and materials. However, the advantages of non-involute shapes are not so significant that they could easily displace involute gears, which are cheaper to manufacture.
Keywords: non-involute gears, tooth profile, path of contact, Novikov gears, S-gears, gearing load capacity Highlights
•	The advantages and disadvantages of established involute gears are analysed.
•	A systematic overview of non-involute spur and helical gears is presented.
•	A comparison between involute and non-involute gears is performed based on established criteria.
•	Results of research on S-gears are given in more detail.
0 INTRODUCTION
To ensure smooth transmission of rotation between two gears, the flanks of the teeth must be designed to comply with the law of gearing: the common normal of the contacting flanks of teeth must pass through the pitch point C (Fig. 1b) at each point of contact. Involute gears with involute-shaped flanks of teeth are predominantly used in industry. Since non-involute gears will be compared with involute ones, the main properties of involute gears and definitions of the basic concepts will be summarised in Chapter 1. Our overview is limited to external gears. Involute gears are described in more detail in gear-books [1] to [3].
The principal disadvantages of external involute gears are:
•	limited minimum number of pinion teeth,
•	convex-convex contact of the tooth flanks and therefore a higher contact load than in the concave-convex contact, which occurs in most non-involute gears,
•	unfavourable kinematic conditions (high sliding speed) at the beginning and end of meshing. Many researchers have attempted to eliminate
these disadvantages by changing the shape of the tooth profile. Non-involute gears can be designed so
that the shape of the tooth flanks of the basic rack is defined first. The path of contact and the tooth flanks of pinion and wheel are defined using the known methods described in the literature [1]. In this way, the S-gears [4] and cosine gears [5] are defined. Another possibility is to define the shape of the tooth flanks based on a pre-selected path of contact. This procedure is used when defining a gearing with a parabolic path of contact [6].
In the follow-up, the theoretical criteria for assessing the load capacity of gears, their efficiency, and noise will be presented first. Next, an overview of non-involute tooth flanks and their comparison with involute teeth will be made. The results of theoretical analyses and tests from the cited literature will be given. The results of research on S-gears that were developed and tested at the Faculty of Mechanical Engineering, University of Ljubljana, will be presented in more detail [7] and [8]. The purpose of this paper is to give a systematic overview of non-involute gears and their advantages and disadvantages. This reference-backed review is useful for both researchers of non-involute gears and for users in practical applications. The geometry and properties of individual non-involute forms are given in Chapter 3. A summary of properties of different non-involute
256
*Corr. Author's Address: University of Ljubljana, Faculty of Mechanical Engineering, Aškerčeva 6, 1000 Ljubljana, Slovenia, ivan.okorn@fs.uni-lj.si
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
gears compared to involute gears, areas of application and references are given in Tables 2 and 3 in Chapter 4.
1 INVOLUTE GEARS
In involute gears, the tooth flank profile has a shape of an involute, which is a curve drawn by a point on a line, which is rolled on the so-called base circle. The involute is mathematically defined by the involute function inv a (Fig. 1a). The curvature of the involute decreases when the radius of the base circle increases. If the radius of the basic circle is infinite, the involute is a straight line. The rack cutter for involute gears has a trapezoidal shape. Tools for manufacturing gears, therefore, have a simple straight-tooth shape, which is one of the most significant advantages of involute gears. The creation of the involute and the meshing of the teeth are shown in Fig. 1b, while the meshing of rack and pinion is shown in Fig. 2. The curve that connects the points of contact of the tooth flanks is called the path of contact. If the shape of the path of contact is defined, the shapes of the rack cutter profile and the flanks of teeth are also defined. In the case of involute gears, the path of contact is a line segment between the initial contact point A and the final contact point E. The extension of the path of contact runs tangentially to the base circle of the gears. The angle between the direction of the normal force to the tooth flank and the direction of the tangent through the pitch point defines the pressure angle a. In involute gearing, the direction of the force on the tooth does not change during meshing, which is favourable with regard to dynamic forces. Involute gears are not sensitive to centre distance errors. If the centre distance in the gearbox is greater than the theoretical one, the backlash increases, the pressure angle increases slightly, while kinematic conditions in the contact of the tooth flanks do not change significantly.
The operating properties of involute gears can be influenced by profile shift (Fig. 2). In a gear without a profile shift, the datum line of the rack runs tangentially to the pitch circle. The profile shift V can be negative or positive displacement of the tooth profile with respect to the axis of rotation of the gear. An undercut appears in an involute gear with a small number of teeth. The undercut limit depends on the helix angle, pressure angle, and height of the flat part of the rack. In the case of spur gears with a pressure angle of 20° and a rack as in Fig. 2, the undercut limit is at 17 teeth. To make a functional gear with a smaller number of teeth, a positive profile shift is needed. The choice of profile shifts is limited by the sharpness of
teeth, so the number of teeth of an involute pinion has its lowest limit.
Fig. 1. Involute gears; a) the formation of an involute, and b) meshing of teeth
Fig. 2. Profile shift at involute gears
The sum of profile shifts is essential for the operating properties and the centre distance of the gear pair. When the sum of profile shifts is positive, the centre distance, pressure angle and the root load capacity of gears increase, while the contact ratio of gears decreases and the noise of gears increases [1] and [2]. The load capacity of involute gears can be increased, and noise and vibrations can be reduced by a tip relief profile modification [9] to [11] and crowning
Operating Performance of External Non-Involute Spur and Helical Gears: A Review
257
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
[12] to [14]. Special software (e.g., KISSsoft) enables the optimization of operating properties of involute gears.
A special design of involute gearing is asymmetric involute gearing. Its geometry and gear design are discussed by Kapelevich [15]. The gearing is not sensitive to centre-distance errors. The rack cutter is not symmetrical and not standardized. The pressure angle at the loading side is greater than on the unloading side, so the asymmetrical gearing has a higher tooth surface durability and higher root load capacity than the symmetrical gearing [16] to [18]. Sekar [19] discusses the friction between the tooth flanks as a function of geometric parameters and load. The author proves that power losses are smaller than with classical involute gears.
2 CRITERIA FOR ASSESSING THE OPERATING PROPERTIES OF GEARS
2.1 Surface Durability of Gears
During the expected lifespan of gears, there should be no tooth fractures or damages to the flanks of teeth. Flank damages include pitting, scuffing, micro pitting and abrasive wear [20]. Pitting is damage made due to fatigue of the tooth flank material. It starts with a crack that can be on the surface or underneath. Because of the cyclic load, the crack widens until a chip is formed, which is flushed out by the oil. Due to the pits formed on the surface of the flanks, the bearing surface of the flanks is reduced, and the contact load on the remaining surface increases. The criterion for the occurrence of pitting is the contact load, defined by the Hertz equation [2]. Evaluation of involute gears according to the ISO 6336-2 standard [21] is based on the following form of the Hertz equation:
o„ =
F,
ld-P-X-
'i-V E
1 -v.
2 v
(1)
J2 J
At each point of contact, the line contact of teeth can be simulated with two cylinders with radii p1 and p2 or with a cylinder and a ring if one of the radii is negative and the contact is concave-convex. The relative radius of curvature is:
P =
Pi 'Pi
(2)
Pi +Pi
The Hertzian pressure between two cylindrical surfaces depends on the normal force Fn, material properties (modulus of elasticity E, Poisson number v), contact length ld and relative radius of curvature.
The relative radius of curvature depends on the shape of the tooth flanks. When designing non-involute gears, the trend is to achieve the largest possible relative radius of curvature, which consequently means a reduction of the Hertzian pressure. In the case of S-gears [4], the relative radius of curvature in the area at the beginning and end of meshing is significantly larger than in involute gears due to the concave-convex contact. However, in the vicinity of the pitch point C, the convex-convex contact also exists in S-gears and, therefore, the relative radius of curvature is approximately the same as in involute gears. It depends on the pressure angle at point C. The same is true for cosine gears [5] and gears with a parabolic path of contact [6]. An analytical method for calculating the curvature radii for gears with a curved path of contact was developed in [7].
For the test conditions, materials, and geometry of the test gears, Hertzian pressures are defined at which pitting occurs after a certain number of tooth flank meshes. Values for involute gear materials are listed in the ISO 6336-5 standard [22]. When evaluating involute gears, the differences between the test conditions (lubricant, speed and the geometry of the test gears) and the actual conditions and geometry are considered by applying additional coefficients [21].
Due to the friction between the tooth flanks, the oil temperature increases, but the oil viscosity and the thickness of the oil film decrease. When the oil temperature exceeds the allowable value, the oil film breaks. When the protective layer is crushed, the two surfaces weld. Due to the sliding of the flanks, the joint collapses, and a gap is formed in the sliding direction. This type of damage is called "scuffing of gears". The evaluation of involute gears on scuffing is based on temperature criteria. Two methods are available in ISO/TS technical specification: the flash temperature method [23] and the integral temperature method [24]. The calculation of the flash temperature is based on Blok's formula [25] and [26].
9fla = 0.62 )0
P
(4BM1 ■ vn
"VBM 2 • Vt 2 )
(3)
Heat generation is influenced by the coefficient of friction p, load (wn=Fn/ ld), relative radius of curvature and sliding speed vg=vt1 -vt2. The denominator of the last term in Eq. (3) takes into account the heat dissipation from the contact. The substitute elasticity
258
Okorn, I. - Nagode, M. - Klemenc, J.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
modulus E' depends on the elasticity modules of both gears and their Poisson numbers.
1
E'
1 -y' 1 -v.
2 \
E
(4)
E
The frequently cited equation [1] and [27] is used for calculating the coefficient of friction.
( \0-2
H = 0.045 ■
V(l + V2 )'P.
• n"0'05 ■ Y
■no ■ yr
(5)
The coefficient of friction depends on load, relative radius of curvature, sum of tangential velocities vt1 + vt2, viscosity of oil n0 and surface roughness (taken into account by the XR coefficient). Non-involute gears with a curved path of contact, such as S-gears [4] and [7], cosine gears [5], and gears with a parabolic path of contact [6], have larger relative radii of curvature and lower sliding speeds at points that are critical for scuffing (at the beginning of meshing). The thermal load is lower at these points, so they are less susceptible to scuffing than involute gears are. Xue et al. [28] combined the dynamic load and transient thermal elastohydrodynamic lubrication (TEHL) to study the scuffing load capacity of spur gears. The hydrodynamic pressure, oil film thickness and flash temperature were calculated with a numerical method. The comparison between the TEHL and Blok's theory was implemented.
The minimum thickness of the oil film at the contact point is crucial for the occurrence of micro-pitting and abrasive wear of the tooth flanks. In general, the thickness of the oil film is calculated by solving the Reynolds differential equation. This procedure was used by Xu et al. [29] in the case of harmonic gear transmission. Typically, the Dawson-Hamrock equation [30] is used for the calculation of the minimum oil-film thickness. The general notation of the equation applies to the elliptical contact.
h = 3.63 + V^
lo
E-p
•(• E')°
/ „ \-0.073 ,	.
•((ft) - 0.61^" ).
(6)
Coeficient k = a/ b is the ellipticity parameter where a and b stand for ellipse axes. If the contact of the tooth flanks is a straight line, the expression in the last bracket of Eq. (6) is equal to 1. The equation can then be rearranged to the following form.
h0 =(3.63-C8-a0'49
E'-0117 )•
,p0A66. F-oal\	(7)
Velocities, the relative radius of curvature, and load vary along the path of contact. Non-involute gears with a curved path of contact, such as S-gears [4] and [7], have a larger relative radius of curvature and a greater sum of tangential velocities at the beginning and end of meshing than the involute ones. The greater thickness of the oil film has a positive effect on both the damage of the flanks and the dynamic forces.
The evaluation of involute gears regarding micro-pitting is covered by the ISO/TS 6336-22 standard [31]. To assess the hazard of the micropitting, the relative thickness of the oil film is used, defined as:

h
0.5 -{Ral + Ra
(8)
where hy is the local film thickness on the assumption of smooth surfaces, while Ra1 and Ra2 are the arithmetic mean roughness values. Bergstedt et al. [32] have demonstrated experimentally that the probability of micropitting is lower on a smoother surface. Clarke et al. [33] and Liu et al. [34] provided an overview of the results of the latest micropitting research.
2.2	Root Load Capacity of Gears
The bending stress at the root of a tooth is an indication of load; the stress depends on the tooth thickness and the tooth fillet at the root. When evaluating involute gears for root load capacity, the bending stress at the root of the tooth is compared with the root fatigue strength, which is determined by gear tests. Analytical evaluation of involute gears is defined in the ISO 6336-3 standard [35]. Numerical methods are predominantly used to calculate the stress at the root of the teeth of non-involute gears, and the finite element method (FEM) is the prevailing numerical method. Analytical methods (such as beamlike models for cycloid gearings) are rare. Numerical analyses are used in many papers dealing with non-involute gears [36] and [37].
2.3	Power Losses and Gear Efficiency
Friction between the flanks of teeth causes power losses. The friction power at the contact point i depends on the normal force Fn, the coefficient of friction f and the sliding speed vg=vt1 - vt2.
P = F
rfri 1 ni
'Pi '
(9)
In non-involute gears with a curved path of contact, such as S-gears [4] and [7], cosine gears [5], and gears with a parabolic path of contact [6], the
2
2
Operating Performance of External Non-Involute Spur and Helical Gears: A Review
259
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
sliding speeds are lower at the beginning and end of meshing; therefore, the power losses are also lower. The friction work during the meshing of one pair of teeth between the initial contact point A and the final contact point E is:
= ¿P Jri-Ati =J Pfr ■ dt.
(10)
With the known speed of the pinion n1 in min-1 and the number of teeth of the pinion z1, the average friction power can be determined.
- Z ■ n
Pfr =-
60
(11)
For involute gears, there is a theoretical method for determining power losses during the meshing of teeth. The procedure is described in the literature [1] and [2]. In the previous decade, an elastohydrodynamic lubrication (EHL) -based friction coefficient was developed and applied to the prediction of gear efficiency [38] and [39]. Ziegltrum et al. [40] compared the simulated load-dependent gear power losses of a transient thermal EHL model with the experimental results in the Forschungsstelle fur Zahnrader und Getriebebau (FZG) gear test rig. Results show a very strong correlation when taking into account the mixed lubrication.
2.4 Vibration and Noise of Gears
The main geometric parameter that affects gear noise is the contact ratio. Mesh stiffness depends on the contact ratio, and it has a decisive influence on vibration and noise level [28]. The gearing is more silent and quieter if the contact ratio is higher. For helical gears, the total contact ratio is the sum of the transverse contact ratio and the face contact ratio.

-p-
(12)
The transverse contact ratio is the ratio between the arc length of action and the angular pitch.
£„ =
C1AC1E C1AC1E
(13)
For involute gears (with points of contact on a line), it is defined as the ratio between the length ACE and the base pitch (Fig. 1). The distance between the initial contact point A and the pitch point C is greater in involute gears, so their transverse contact ratio is higher. Non-involute gears with a curved path of contact, such as S-gears [4] and [7], cosine gears [5], and gears with a parabolic path of contact [6], have a
shorter distance between the initial meshing point A and the pitch point C, so the contact ratio is lower than that of involute gears. The contact ratio of involute gears can be changed by making a profile shift. If the sum of profile shifts is negative, the contact ratio is higher than if the sum is equal to 0 or positive. The overlap ratio is higher if a gear is wider and if the helix angle is larger.
The force on the tooth acts perpendicularly to the tooth flank profile. Because the path of contact is a straight line at involute gears, the direction of the force does not change during meshing, which has a good effect on vibrations. Vibration and noise of gears are decisively affected by tooth deviations [41]. The noise of involute and non-involute gears can only be adequately compared if the gears have been made in the same accuracy class. Greater manufacturing accuracy significantly reduces vibration and noise levels.
3 NON-INVOLUTE GEARS
Various forms of non-involute gears are presented in this chapter. Operating properties are described for each gear type, supported by theoretical and experimental research. A comparison with comparable involute gears is given.
3.1 Cycloid Gearing
The flank of the cycloid gear tooth consists of a hypo- and an epicycloid [3]. A cycloid is a curve produced by tracing a path of a chosen point on the circumference of a circle that rolls around a fixed base circle. The shape of the addendum flank of the tooth is obtained by rolling the circle on the outer side of the base circle, while the shape of the dedendum flank of the tooth is obtained by rolling a circle inside the base circle (Fig. 3). During meshing, the hypocycloid and the epicycloid touch, the touch is concave-convex, and the path of contact has a shape of circular arcs. The rack has flanks curved in the shape of an orthocycloid, which is obtained by rolling a circle along a straight line. The shape of the hypocycloid depends on the ratio of diameters of the rolling circle and the base circle (kc). The hypocycloid is a straight line if the value of kc=0.5. The contact of the hypo- and the epicycloid is concave-convex when kc< 0.5, so the contact pressure is lower than that of involute gears. In these cases, the relative sliding is also lower than with comparable involute gears. The critical point of the gearing is the pitch point C, where a transition from the epicycloid to the hypocycloid occurs. In the
A
m -n
n
260
Okorn, I. - Nagode, M. - Klemenc, J.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
area of the pitch point, the relative curvature radii are very small, while the contact pressure is higher than in involute gears.
Fig. 3. Cycloid gears
A pinion can have a very small number of teeth (theoretically only two or three). Cycloid gears with two or three teeth are used for rotors in Roots blowers. Gearing is very sensitive to centre distance errors and does not allow a profile shift. If there is an error in the centre distance, vibrations increase. Manufacturing gears is a complicated and expensive process. This type of gears is used in precision mechanics (e.g., mechanical clocks).
i of contact
Fig. 4. Lantern gears
One special example of cycloid gearing is lantern gearing (Fig. 4) with pins on one of the gears. The flank of the gear tooth has the shape of an equidistant cycloid. The contact point on the pin slides along the flank of the gear tooth during meshing. The lantern gearing is only suitable for small circumferential speeds, e.g., drives with large driven gears in transport
devices. A special shape of the cycloid gear is installed in the rotary vector cycloidal-pin drive [42], used in robots. The main feature of such a gearbox is a large gear ratio.
Peng et al. [43] developed a new form of cycloid gearing named arc-tooth-trace cycloid gear (Fig. 5).
Fig. 5. Arc-tooth-trace cycloid gear; a) tooth profile, and b) arc-tooth-trace; adapted from [43]
The tooth trace of the gear is a segment of a circular arc. The profile of the tooth flank consists of two parts: above the reference circle is a circular arc segment, while a cycloid segment is below the circle. Equations for describing the geometry of the gearing and calculation of the contact ratio are derived in the paper. The transverse contact ratio is a constant value (0.5), while the overlap ratio depends on the module and the radius of a circular arc. 3D models of the new gears and comparable involute gears were made. The results of the numerical analysis show that the new gearing has 17 % higher tooth surface durability and 35 % higher root load capacity than involute gearing.
3.2 Helical gears with circular arc teeth
In 1926, Wildhaber [44] invented helical gearing with circular arc teeth (in a normal cross-section). Independently, a similar shape (a tooth profile with a circular arc in the front face) was invented 30 years later by Novikov [45]. Due to a similar shape, researchers in the past had named the gearing with the
Operating Performance of External Non-Involute Spur and Helical Gears: A Review
261
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
sides of teeth in the form of circular arcs as Wildhaber-Novikov gearing, or shorter W-N gearing. Radzevich [46] and [47] clarifies differences between the two systems and explains why Wildhaber and Novikov gearing should be addressed separately. In his opinion, the name "W-N gearing" is incorrect. The geometry of the original Novikov gearing is shown in Fig. 6. This kind of gearing was applied in the final reduction stage of a Westland Lynx helicopter gearbox [48].
Fig. 6. Original Novikov gearing
The touch of the teeth is concave-convex. The contact surface has (in theory) an elliptical shape under load. When two gears rotate, the contact surface travels in the axial direction, parallel to the axis of the gear rotation. Transmission of rotation is continuous if the overlap ratio is higher than 1. To ensure the lateral backlash, the radius of the pinion's circular arc must be slightly larger than the circular arc radius of the gear wheel.
Dyson et al. [49] and [50] discuss the gearing geometry and kinematics in detail and present design instructions. Contact stress during meshing was analysed in more detail by Coulbourne [51]. Due to the concave-convex contact, contact stress is lower than in involute gearing. According to Niemman's research [52], W-N gearing has 1.5 to 3 times higher surface durability than comparable involute gearing. These values apply to a gear shape when px = p2. Differences decrease with an increasing helix angle. The influence of pressure angle and root fillet on the stress in the tooth root is discussed by Tsay [53]. His analysis was performed using FEM. Results of analyses indicate
that the W-N gearing has a slightly lower root load capacity than the involute gearing. The conditions for the formation of the oil film are better than in involute gearing [54].
Over the previous decade, Markovski and Batsch [55] have conducted numerical and experimental research on Novikov gearing. They incorporated errors that occur during manufacturing and assembly into their mathematical model. On the basis of numerical contact analysis, they analysed the impact of errors on the wear of tooth flanks. It turned out that the contact surface is not a real ellipse as it was theoretically assumed. They developed a method for measuring and displaying the contact surface of Novikov gearing [56]. They performed fatigue tests on Novikov gears (mn = 3 mm, zj = 30, z2 = 47, b = 30 mm) and comparable involute gears [57]. With the same torque on the pinion, the Novikov gears lasted three times longer till the first signs of pitting. This confirms that the surface durability of the Novikov gearing is about three times higher than the durability of involute gearing. They also measured vibration amplitudes on the housing. They were approximately five times higher than in the involute ones. The main reasons for this are slightly worse manufacturing accuracy and a small contact ratio. Batsch extended the research of Novikov gears to bevel gears [58]. A mathematical model of meshing is derived. Contact stresses and maximum bending stresses in Novikov gears and comparable spiral bevel gears are calculated numerically. In Novikov gears, the contact area is larger, so the calculated contact stresses are lower by approximately 22 %. Root stresses are higher by 27 % on the Novikov pinion.
Litvin et al. [59] developed a new shape of the W-N helical gearing. The rack cutter has a parabolic shape instead of a circular shape (Fig. 7). The contact ellipse is longer and narrower than in the basic Novikov gearing. Advantages of the new W-N gearing compared to the basic Novikov gearing are lower noise and vibration, lower sensitivity to manufacturing and assembly errors, the possibility of grinding the gears, lower contact pressure and lower bending stress. The advantages were proven by a simulation of meshing and by numerical calculations of contact stresses and bending stresses.
W-N gearing can also be made with two meshing zones (Fig. 8a) [60]. Compared to gearing with a single-zone meshing, this one has lower meshing stiffness. The jerks at the beginning of the meshing are smaller, and the distribution of force on the teeth is more favourable. A rack cutter with two circular arcs for manufacturing such gearing is shown in Fig. 8a.
262
Okorn, I. - Nagode, M. - Klemenc, J.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
Fig. 7. New version of the W-N gear; a) rack-cutter in normal section, and b) transverse profiles; adapted from [59]
Fig. 8. Profiles of rack-cutter; a) W-N gears with two zones of meshing, and b) stepped triple circular-arc gears; adapted from [61]
The generated tooth flank has a singular point due to the sharp transition between the radii. To eliminate it, Yang [61] proposed a rack cutter with three circular
arcs, shown in Fig. 8b. He named it a stepped triple circular-arc gearing. He developed a mathematical model of the gearing and performed a stress analysis using FEM.
Ariga and Nagata [62] also suggested an improvement of Wildhaber-Novikov gearing. The top of the new tooth profile has the shape of a circular arc, and the dedendum consists of two involutes. Tests have shown that these gears are less sensitive to centre-distance errors than conventional W-N gears.
The gearing with a gear profile composed of circular arcs was also developed by Hlebanja and Hlebanja [63] (Fig. 9), called uniform power transmission gears (UPTG). A theoretical comparison with involute gearing is made. To confirm its functionality, a model of aluminium gears was made.
Fig. 9. UPT gears [63]; a) shape of the rack-cutter, and b) pair of UPT
Liang et al. [64] proposed gear transmission with a double circular arc-involute tooth profile. Three components are important in the generation of a tooth profile: the upper convex circular arc, the middle involute and the lower convex circular arc. Based on the equations of tooth surfaces, the software for calculating the data points of tooth surfaces was developed. The results were exported to the 3D drawing software. A numerical comparison of contact stresses and sliding coefficients between the involute
Operating Performance of External Non-Involute Spur and Helical Gears: A Review
263
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
gearing and the new gearing was made; both are lower in the new gearing. Proposed gears and involute gears were also compared experimentally. The transmission efficiency under different operating conditions is presented. The measured efficiency of the new gear was lower due to the lower manufacturing quality.
3.3	Pure Rolling Helical Gears
Chen et al. [65] developed novel circular arc helical gears (CAHGM) with pure rolling. Based on the meshing principle of space curve meshing, the parameter equations of the concave-convex circular arc profiles were established. A pair of CAHGM was manufactured via rapid prototyping technology. The contact ratio equation and the theoretical kinematic performance were validated. The gears have a large backlash which is harmful to the forward and reverse drives. Chen et al. [66] and [67] also developed pure rolling helical gears with convex-to-convex meshing type and pure rolling rack and pinion mechanism. Geometric design, meshing performance, and mechanical behaviour were presented. Parametric equations for contact curves and for the tooth surfaces were derived. A numerical comparison of contact stresses and bending stresses between the modified involute gearing and the new gearing was made. With new gearing, the maximum bending stresses are significantly lower, and the contact stresses are 1.5 to 2 times higher. The maximum contact stress can be decreased by increasing the tooth number and face width and decreasing the helix angle. Properties of pure rolling gearing have not yet been experimentally researched.
3.4	Convoloid Gearing
Berlinger and Coulbourne [68] developed a new type of gearing, named convoloid gearing. Its tooth's dedendum has a concave shape, while its addendum has a convex shape (Fig. 10). The transverse contact ratio is between 1.1 and 1.3, while the overlap ratio is an integer. A theoretical and experimental comparison with equivalent involute gears is made. An important geometric parameter is the ratio between the width of the gear b and centre distance a. If b / a< 0.2, the convoloid gearing has no advantages. The load capacity ratio between convoloid and involute gears increases with increasing b / a ratio. Test results show that 20 % to 35 % higher load capacity can be achieved with convoloid gears.
J>— outside diameter - start of tip break addendum mid point start of addendum relief
addendum transition point pitch radius
dedendum transition point start of dedendum relief
dedendum mid point start of root break profil control radius root radius
Fig. 10. Convoloid gearing; adapted from [68]
3.5 Gears with Curved Paths of Contact 3.5.1 S-Gears
Development of S-gears (invented by Hlebanja [69] and developed by Hlebanja et al. [70]) began in the 1970s. This type of gear was used in a rolling mill with a power of 1500 kW [71]. A path of contact in the shape of the letter S is a characteristic of this gearing. The first gearing had a small pressure angle at the pitch point C compared to today's shape. In the 1990s, the basic rack tooth profile equation was defined [72].
y = ap-(i-(1 - x)").	(14)
The following parameters were selected in [7]: ap = 1.3 and n = 1.9. The shape of the tooth profile, the path of contact, and both tooth flanks are shown in Fig. 11. A tool for shaping S-gears with a module of 4.575 mm was made. Research studies of load capacity of spur S-gears made of steel and comparable involute spur gears were performed. Involute gears had a module mn = 4.5 mm, width b = 20 mm, number of teeth z1 = 16 and z2 = 24, and the profile shift coefficients x1 = 0.233 and x2 = 0.12. Before the tests, the geometry of all gears was measured on a coordinate-measuring machine. S-gears were
264
Okorn, I. - Nagode, M. - Klemenc, J.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
manufactured in accuracy classes 7 and 8, and involute gears in class 6 (according to DIN 3962 [73]). The roughness of the tooth flanks was approximately the same in both types of gears.
Fig. 11. S-gears [7]; a) basic rack tooth profile, and b) tooth flanks of S-gears
The load capacity tests were made on a standard FZG test rig [74] with a centre distance of 91.5 mm. Theoretical calculations of Hertzian pressures, flash temperatures and oil film thicknesses were made for the geometry of the test gears and the test conditions. Figs. 12 and 13 show a comparison of Hertzian pressures and flash temperatures along the path of contact. Theories about the load capacity of S-gears were defined on the basis of the criteria described in Chapter 2.
The tests showed a significantly higher resistance of S-gears to scuffing, which is the expected result on the basis of the calculation of contact temperatures and oil film thickness at the beginning of meshing. In fatigue tests, pitting occurred in the pitch point area of S-gears, and in the dedendum of involute gears. S-gears withstood a slightly higher number of load cycles than involute ones (by approx. 20 %) [7] and [75]. The wear of S-gears was about one half lower during the same operating time due to a thicker oil film.
Fig. 12. Hertzian pressures during tests; a) involute gears, and b) S-Gears
Fig. 13. Flash temperatures during tests; a) involute gears, and b) S-Gears
The heat load of the gears was estimated on the basis of oil temperature measurements [76] and [77]. S-gears with a modulus of 4.575 mm and involute gears with a modulus of 4.5 mm, and the number of teeth of 16 and 24 were used in the tests. The oil temperature at the beginning of the test was 90 °C. Oil temperatures after 15 minutes of operation at a given load in the standard FZG scuffing test are shown in Table 1. The measured temperatures prove that the heat load in S-gears is lower than in involute ones.
Operating Performance of External Non-Involute Spur and Helical Gears: A Review
265
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
Table 1. Oil temperature in °C at the end of the test at each load step
Load step	7	8	9	10	11	12
Torque T [Nm]	183	239	302	372	450	543
S-gears	94	97	102	106	113	118
Involute gears	99	103	107	112	120	127
One of the disadvantages of S-gears is a lower contact ratio than with involute gears. This weakness can be eliminated by using helical S-gears depicted in Fig. 14a. A gear pair with helical S-gears is shown in Fig. 14b. These gears were manufactured with the same shaping tool as the test spur S-gears. No experimental research on these gears has been made yet, but it is expected that operating characteristics will be better with helical gears than with spur S-gears.
Fig. 14. Helical S-gears; a) contact of the teeth, and b) a gear pair
The S-gears can also be used in internal gears, which are widely used in planetary gears. The analysis of the load capacity of the internal S-gears is discussed by Hlebanja and Hlebanja [78]. Since the path of contact is a straight line in the pitch point region (which is similar to involute gears), the S-gears are significantly less sensitive to centre distance errors [79] than classic cycloid gearing. The operating properties of the S-gears depend on parameters ap and n. The influence of these parameters on the shape of teeth and operating properties is discussed by Kulovec and Duhovnik [80]. The S shape of the path of contact
was also used in worm gears. Theoretical analysis of the situation and demonstration of the operation is given by Hlebanja et al. [81].
In the previous decade, the concept of using polymer S-gears has emerged. Polymer gears have a low modulus of elasticity, so the deformations are significantly higher. Most gears are made by injection moulding, with which it is not possible to achieve high precision. The contact temperature (which depends on the friction coefficient) has a significant impact on damage to polymer gears. The theoretical contact temperature of polymer S-gears is lower than that of involute gears with comparable geometry (same modulus and number of teeth), so the expected load-bearing capacity of S-gears is higher. In the first phase of research, the gears were made by injection moulding. The advantages of S-gears have not been demonstrated due to poor manufacturing accuracy and large deformations. In the second phase of research, the gears were made by hobbing. In this case, a higher load capacity of the S-gears was measured. Results on research concerning polymer S-gears were published in several papers [82] to [85]. A comparison of the efficiency of S-gears and involute gears was made. The results show better efficiency of S-gears [86]. Problems of polymeric S-gears are large deformations, small size of teeth (test gears had a module of 1 mm) and manufacturing accuracy. The problem is general and specific to plastic S-gears. In steel gears, the deformations do not have such an effect on meshing. The accuracy class of the geometry is also higher for steel gears. Advantages of S-gears can be seen only if manufacturing accuracy is high enough, which can be achieved by hobbing. The question is, whether energy savings and higher load capacity justify the higher cost of manufacturing polymer S-gears.
3.5.2 Cosine Gearing
The shape of the flank of a pinion with a cosine gearing is described by a cosine function. Luo et al. [5] derived equations for the description of the tooth flank and the path of contact. A 3D model of the gearing was made, and a numerical analysis was performed with FEM. Contact ratio and sliding coefficient were calculated. Compared to involute gearing, the cosine gearing, which is discussed in [5], has the following characteristics: lower contact pressure (by approx. 22 %), lower root stress (by approx. 35 %), lower contact ratio (by 20 %) and considerably smaller sliding coefficient, both on gear and pinion. The stated properties of the gearing are valid generally, but the values in parentheses depend on the geometric
266
Okorn, I. - Nagode, M. - Klemenc, J.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
parameters of the gears. Wadagaonkar and Shinde [87] confirmed lower contact stresses (by approx. 30 %) and lower bending stresses (by approx. 50 %) for their geometry of cosine gears.
3.5.3	Gearing with a Parabolic Path of Contact
Wang et al. [6] worked on gears with a parabolic path of contact. A mathematical model of the gearing was derived on the basis of a known path of contact. The influence of the parabola parameters on the shape of teeth of the pinion and the gear was analysed. The number of teeth in which dedendum undercut occurs is significantly smaller than in the case of involute gearing. A numerical comparison of contact stresses and bending stresses between the involute gearing and the new gearing was made; both are lower in the new gearing. This general statement is valid for any form of a parabola, while the quantitative values depend on the parabola parameters that can be optimized. Properties of a gearing with a parabolic path of contact have not yet been experimentally researched.
3.5.4	Gears with a Constant Relative Radius of Curvature
The relative radius of curvature has a decisive effect on the Hertzian contact pressure (Eq. (1)). Liu et al. [88] developed a new form of non-involute gearing based on a controlled relative radius of curvature. It was named constant relative curvature (CRC) gears. Equations were developed to calculate the profile of the rack cutter and gears. A mathematical model for the gearing was defined with a constant relative radius of curvature along the path of contact. The tooth profile is shown in Fig. 15. CRC and involute gears were compared analytically. The contact ratio was calculated for both types of gears; it was slightly lower for CRC gears. Contact stress, sliding coefficients, and oil film thicknesses were determined. The gearing has similar advantages over the involute gears as the
involute gear CRC gear
Fig. 15. Tooth profile of CRC gear and involute gear; adapted from [88]
S-gears: lower contact stress, thicker oil film and less sliding of flanks. The advantages of CRC gears have not yet been confirmed experimentally.
4 DISCUSSION
Gears in modern gearboxes must have the best possible efficiency, high load capacity, and low noise level combined with the smallest possible size. In recent decades, research has focused mainly on improving the properties of involute gears. The effects of geometry, lubrication, manufacturing accuracy and new materials on gearing load capacity and noise were researched in detail. New findings in tribology have significantly contributed to the improvement of the operational characteristics of involute gears.
Research on non-involute forms is, therefore, less frequent and less up-to-date. The development of software increased the possibilities of accurate analyses of non-involute forms. Mathematical models for defining geometry parameters have been developed to describe various non-involute shapes mentioned in this paper. A mathematical model of the gearing enables the creation of a 3D model of the gear pair, which can be analysed numerically. Analyses are relatively inexpensive and quickly feasible on modern computers. Most analyses are performed for the ideal shape of the tooth flanks. The accuracy of manufacturing and deformations of teeth significantly affect the shape of teeth and the operating characteristics. Involute and non-involute gears can only be compared if their manufacturing accuracy is comparable.
Based on the research of gears with a curved path of contact and the results of research on non-involute gears conducted and published by other authors, Table 2 assesses the properties of non-involute gears. The comparison refers to a comparable involute gear (the same number of teeth, helix angle, modulus, and accuracy class).
The three black dots in Table 2 indicate that the non-involute gearing is slightly better than the involute one in terms of the given property, while four and five dots indicate that it is significantly better than the involute one. At two black dots or less, the non-involute gearing is worse than the involute one.
At the Faculty of Mechanical Engineering, University of Ljubljana, research was carried out on steel gears with a progressively curved path of contact [7]. Prototypes of S-gears and comparable involute gears were made. A theoretical comparison was performed according to the criteria described in Chapter 2. The theoretical assumptions were
Operating Performance of External Non-Involute Spur and Helical Gears: A Review
267
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
Table 2. Evaluation of the properties of non-involute gears compared to involute gears (•
> best; ooooo worst)
_ . n .. . _ „. . Contact ratio,	Sensitivity to
Surface Root load Gear efficiency, .... .. .	. . .■' .
Vibration and	manufacturing and
durability capacity Heat load
noise	assembly errors
References
Cycloid gears	• ••oo	•••oo	•••oo	••ooo	ooooo	[1], [3], [42]
Lantern gears	• •ooo	••ooo	ooooo	••ooo	•oooo	[1], [3]
Arc-tooth-trace cycloid gears	• ••oo	•••oo	•••oo	••ooo	••ooo	[43]
Original Novikov gears	• ••oo	•••oo	•••oo	•oooo	ooooo	[44] to [58], [60]
New W-N gears	••••o	•••oo	•••oo	••ooo	••ooo	[36], [59], [62]
Stepped triple circular-arc gears	•••oo	•••oo	•••oo	••ooo	•oooo	[61]
UPT Gears	•••oo	•••oo	•••oo	••ooo	••ooo	[63]
Double circular arc-involute tooth profile	•••oo	•••oo	••ooo	••ooo	••ooo	[64]
Pure rolling helical gears	•••oo	•••oo	•••••	••ooo	••ooo	[65], [66], [67]
Convoloid gears	••••o	•••oo	•••oo	••ooo	••ooo	[68]
S-Gears	••••o	•••oo	••••o	••ooo	••ooo	[4], [7], [8], [69] to [72], [75] to [86]
Cosine gears	•••oo	•••oo	•••oo	••ooo	••ooo	[5], [87]
Gears with a parabolic path of contact	•••oo	•••oo	•••oo	••ooo	•oooo	[6]
CRC-Gears	••••o	•••oo	•••oo	••ooo	••ooo	[88]
confirmed by tests at the FZG test rig. In the last decade, extensive research has also been conducted on plastic S-gears [8]. In papers on cosine gears [5], gears with a parabolic path of contact [6] and CRC gears [87], researchers listed similar properties as were found in S-gears.
Non-involute gears are relatively well researched analytically and numerically in terms of load capacity, while vibration and noise analyses remain lacking. We also lack experimental confirmations of theoretical assumptions regarding load capacity, efficiency, and noise.
Table 3. Areas of use of some non-involute gears
Areas of use
Cycloid gears Cycloidal-pin drive
watches, roots blowers, gear pumps robots (large gear ratio)
Lantern gears
transport devices (large gear ratio, small speeds)
W-N gearing
aircraft industry (helicopters), oil industry
Convoloid gearing
wind turbines
S-Gears
rolling mill, wind turbines, planetary gearboxes, precision mechanics
Cosine gearing
gear pumps
Gearing with a parabolic path of contact
gear pumps
Areas of use of some non-involute gears are listed in Table 3. The gears presented in recent years are still in the research phase and have not yet been used in practical applications.
5 CONCLUSIONS
The paper presents both non-involute gears with a long history (cycloid gears, classical Novikov gears) and gears first presented in the last decade. New technologies (3D printing) enable easy and fast prototyping and testing of the functionality of new designs. Unfortunately, prototypes are not enough to convince users that new designs are useful.
Different gears with a curved path of contact, which show similar operating properties compared with involute gears (higher root load capacity and surface durability, lower tooth sliding and heat load, lower contact ratio), are discussed in more detail. Theoretical findings for spur S-gears have also been experimentally confirmed. Vibration and gear noise research is lacking. A higher contact ratio can be achieved with helical S-gears, and it would make sense to research them experimentally in the future.
Experimental studies on gears manufactured in real-life accuracy and surface roughness are rare. We believe that it is not possible to succeed in the market with a gearing that is not supported by experiments.
Therefore, the geometry and manufacturing of involute gears are much simpler and, above all, cheaper than that of non-involute gears. The evaluation of involute gears is carried out according to standards and recommendations. The references only list ISO standards and ISO technical specifications, but national standards (e.g., AGMA, DIN, etc.) and recommendations also exist (e.g., VDI).
256
Okorn, I. - Nagode, M. - Klemenc, J.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
Because it is much more difficult to achieve the same accuracy with non-involute gears, they are still limited to specific applications. The current situation in development and research indicates that non-involute gears cannot displace involute gears on a large scale.
6 ACKNOWLEDGMENT
The authors acknowledge the financial support from the Slovenian Research Agency (research core funding No. P2-0182).
7 REFERENCES
[1]	Niemann, G., Winter, H. (2003). Maschinenelemente: Band 2. Springer Verlag, Berlin, DOI:10.1007/978-3-662-11873-3.
[2]	Linke, H., Börner, J., Heß, R. (2016). Cylindrical gears: Calculation - Materials - Manufacturing. Carl Hanser Verlag, München, D0l:10.3139/9781569904909.
[3]	Litvin, F.L., Fuentes, A. (2004). Gear Geometry and Applied Theory. Cambridge University Press, New York, D0I:10.1017/ CB09780511547126.
[4]	Hlebanja, J., Okorn, I. (1999). Characteristische Eigenschaften von Zahnräder mit statig gekrümter Eingriffslinie. Antriebstechnick, vol. 38, no. 12, p. 55-58.
[5]	Lou, S.M., Wu, Y., Wang, J. (2008). The generation principle and mathematical models of a novel cosine gear drive. Mechanism and Machine Theory, vol. 43, p. 1543-1556, D0I:10.1016/j.mechmachtheory.2007.12.007.
[6]	Wang, J., Hou, L., Luo, S., Wu, R.Y. (2013). Active design of tooth profiles using parabolic curve as the line of action. Mechanism and Machine Theory, vol. 67, p. 47-63, D0I:10.1016/j.mechmachtheory.2013.04.002.
[7]	Okorn, I. (2000). Research of Tooth Flank's Durability of Gears with Progressive Curved Path of Contact. PhD Thesis, University of Ljubljana, Ljubljana. (in Slovene)
[8]	Zorko, D. (2019). Engineering Design of Polymer Gears with S Form of Teeth. PhD Thesis, University of Ljubljana, Ljubljana. (in Slovene)
[9]	Wagaj, P., Kahraman, A. (2002). Influence of tooth profile modification on helical gear durability. Journal of Mechanical Design, vol. 124, no. 3, p. 501-510, D0I:10.1115/1.1485289.
[10]	Sánchez, M.B., Pleguezuelos, M., Pedrero, J.I. (2019). Influence of profile modifications on meshing stiffness, load sharing and transmision error of involute spur gears. Mechanism and Machine Theory, vol. 139, p. 506-526, D0I:10.1016/j.mechmachtheory.2019.05.014.
[11]	Wang, Z.G., Chen, Y.C. (2020). Design of helical gear set with adequate linear tip-relief leading to improved static and dynamic characteristics. Mechanism and Machine Theory, vol. 147, art. ID. 103742, D0I:10.1016/j. mechmachtheory.2019.103742.
[12]	Bergseth, E., Björklund, S. (2010). Logarithmical crowning for spur gears. Strojniški vestnik - Journal of Mechanical Engineering, vol. 56, no. 4, p. 239-244.
[13]	Wang, N., Li, X., Wang, K., Zeng, Q., Shen, X. (2017). A novel axial modification and simulation analysis of involute spur
gear. Strojniški vestnik - Journal of Mechanical Engineering, vol. 63, no. 12, p. 736-745, D0l:10.5545/sv-jme.2017.4307.
[14]	Weinberger, U., Otto, M., Stahl, K. (2020). Closed-form calculation of lead flank modification proposal for spur and helical gear stages. Journal of Mechanical Design, vol. 142, no. 3, art. ID. 031106, D0I:10.1115/1.4045396.
[15]	Kapelevich, A.L. (2000). Geometry and design of involute spur gears with asymmetric teeth. Mechanism and Machine Theory, vol. 35, no. 1, p. 117-130, D0I:10.1016/S0094-114X(99)00002-6.
[16]	Mo, S., Ma, S., Jin, G. (2019). Research on composite bending stress of asymmetric gear in consideration of friction. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 233, no. 8, p. 2939-2959, D0I:10.1177/0954406218797975.
[17]	Mo, S., Ma, S., Jin, G. (2019). Design principle and modeling method of asymmetric involute internal helical gears. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 233, no. 1, p. 244-255, D0I:10.1177/0954406218756443.
[18]	Marimuthu, P., Muthuveerappan, G. (2016). Design of asymmetric normal contact ratio spur gear drive through direct design to enhance the load carrying capacity. Mechanism and Machine Theory, vol. 95, p. 22-34, D0I:10.1016/j. mechmachtheory.2015.08.013.
[19]	Sekar, R.P. (2019). Determination of load dependent loss factor on asymmetric spur gear. Mechanism and Machine Theory, vol. 135, p. 322-335, D0I:10.1016/j. mechmachtheory.2019.02.011.
[20]	Snidle, R.W., Evans, H.P. (2009). Some aspects of gear tribology. Proceedings IMechE, Part C: Journal of Mechanical Engineering Science, vol. 223, no. 1, p. 103-141, D0I:10.1243/09544062JMES1168.
[21]	ISO 6336-2:2019. Calculation of load capacity of spur and helical gears - Part 2: Calculation of surface durability (pitting). International Organization for Standardization, Geneva.
[22]	ISO 6336-5:2016. Calculation of load capacity of spur and helical gears - Part 5: Strength and quality of materials. International Organization for Standardization. Geneva.
[23]	ISO/TS 6336-20:2017. Calculation of load capacity of spur and helical gears - Part 20: Calculation of scuffing load capacity (also applicable to bevel and hypoid gears) -Flash temperature method. International Organization for Standardization. Geneva.
[24]	ISO/TS 6336-21:2017. Calculation of load capacity of spur and helical gears - Part 21: Calculation of scuffing load capacity (also applicable to bevel and hypoid gears) -Integral temperature method. International Organization for Standardization. Geneva.
[25]	Blok, H. (1963). The flash temperature concept. Wear, vol. 6, no. 6 p. 483-494, D0I:10.1016/0043-1648(63)90283-7.
[26]	Theyse, F.H. (1967). Die Blitztemperaturhypotese nach Blok und ihre praktische Anwenung bei Zahnräder. Schmierungstechnick, vol. 14, p. 22-29.
[27]	Yu, L., Wang, G., Zou, S. (2017). The calculation of meshing efficiency of a new type of conical involute gear. Strojniški vestnik - Journal of Mechanical Engineering, vol. 63, no. 5, p. 320-330, D0I:10.5545/sv-jme.2016.3843.
Operating Performance of External Non-Involute Spur and Helical Gears: A Review
388
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
[28]	Xue, J., Li, W., Qin, C. (2014). The scuffing load capacity of involute spur gear systems based on dynamic loads and transient thermal elastohydrodynamic lubrication. Tribology International, vol. 79, p. 74,83, D0l:10.1016/j. triboint.2014.05.024.
[29]	Xu, X., Fan, X., Wei, P., Yang, B. (2020). Research on the lubrication characteristics of harmonic gear transmission meshing areas. Strojniški vestnik - Journal of Mechanical Engineering, vol. 66, no. 9, p. 513-522, D0I:10.5545/sv-jme.2020.6667.
[30]	Dowson, D., Higginson, G.R. (1977). Elastohydrodinamic Lubrication. Pergamon Press, Oxford.
[31]	ISO/TS 6336-22:2018. Calculation of load capacity of spur and helical gears - Part 22: Calculation of micropitting load capacity. International Organization for Standardization. Geneva.
[32]	Bergstedt, E., Lin, J., Olofssen, U. (2020). Influence of gear surface roughness on the pitting and micropitting life. Proceedings of the Institution of Mechanical Engineers,, Part C: Journal of Mechanical Engineering Science, vol. 234, no. 24, p. 4953-4961, D0I:10.1177/0954406220931541.
[33]	Clarke, A., Evans, H.P., Snidle, R.W. (2015). Understanding micropitting in gears. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 230, no. 7-8, p. 1276-1289, D0I:10.1177/0954406215606934.
[34]	Liu, H., Liu, H., Zhu, C., Zhou, Y. (2019). A review on micropitting studies of steel gears. Coatings, vol. 9, no. 1, art. ID 42, D0I:10.3390/coatings9010042.
[35]	ISO 6336-3:2019. Calculation of load capacity of spur and helical gears - Part 3: Calculation of tooth bending strength. International Organization for Standardization. Geneva.
[36]	Litvin, F.L., Fuentes, A., Gonzalez-Perez, A.I., Carnevali, L., Sep, T.M. (2000). New version of Novikov-Wildhaber helical gears: computerised design, simulation of meshing and stress analysis. Computer Methods in Applied Mechanics and Engineering, vol. 191, no. 49-50, p. 57 07-5740, D0I:10.1016/ S0045-7825(02)00482-6.
[37]	Zorko, D., Tavčar, J., Duhovnik, J. (2018). The influence of the tooth profile shape on the stress-strain in the gear. International Journal for Science, Technics and Innovations for the Industry: Machines, Tecnologies, Materials, vol. 12, no. 4, p. 153-156.
[38]	Li, S. (2015). A thermal tribo-dynamic mechanical power loss model for spur gear pairs. Tribology International, vol. 88, p. 170-178, D0I:10.1016/j.triboint.2015.03.022.
[39]	Liu, M., Xu, P., Zhang, J., Ding, H. (2019). Analysis of mechanical power loss of a helical gear pair based on the starved thermalelastohydrodynamic lubrication model. Industrial Lubrication and Tribology, vol. 72, no. 3, p. 333340, D0I:10.1108/ILT-06-2019-0225.
[40]	Zigeltrum, A., Lohner, T., Stahl, K. (2017). TEHL simulation on the influence of lubricants on load-dependent gear losses. Tribology International, vol. 113, p. 252-261, D0I:10.1016/j. triboint.2016.12.018.
[41]	Shuai, M., Shuai, M., Guoguang, J., Xiang, Y.Z. (2019). Research on dynamic load-sharing characteristics of two-stage
asymmetric star gear system. IEEE Access, vol. 7, no. 1, p. 126799-126811, D0I:10.1109/ACCESS.2019.2937940.
[42]	Li, T., Xu, H., Tian, M. (2020). A loaded analysis method for RV cycloidal-pin transmission based on the minimum energy principle. Strojniški vestnik - Journal of Mechanical Engineering, vol. 66, no. 11, p. 655-667, D0I:10.5545/sv-jme.2020.6868.
[43]	Peng, Y., Song, A., Shan, Y., Lin, X. (2017). A novel arc-tooth-trace cycloid cylindrical gear. Mechanism and Machine Theory, vol. 118, p. 180-193, D0I:10.1016/j. mechmachtheory.2017.08.009.
[44]	Wildhaber, E. (1926). Helical gearing, Patent 1601750. United States Patent and Trademark Office, Washington.
[45]	Novikov, M.L. (1956). USSR Patent 109750, USSR, Moscow.
[46]	Radzevich, S.P. (2020). High-conformal Gearing: Kinematics and Geometry. Elsevier, Amsterdam.
[47]	Radzevich, S.P. (2018). An examination of high-conformal gearing. Gear Solutions, p. 31-39.
[48]	Shotter, B.A. (1978). The Lynx transmission and conformal gearing. SAE Technical Paper, art ID 781041, D0I:10.4271/781041.
[49]	Dyson, A., Evans, H., Snidle, R. (1986). Wildhaber-Novikov circular arc gears: Geometry and kinematics. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, vol. 403, no. 1825, p. 313-340, D0I:10.1098/rspa.1986.0014.
[50]	Dyson, A., Evans, H., Snidle, R. (1989). Wildhaber-Novikov circular arc gears: Some properties of relevance to their design. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, vol. 425, no. 1869, p. 341-363, D0I:10.1098/rspa.1989.0110.
[51]	Coulbourne, J.R. (1989). The contact stresses in Novikov gears. Mechanism and Machine Theory, vol. 24, no. 3, p. 223229, D0I:10.1016/0094-114X(89)90061-X.
[52]	Niemann, G. (1961). Novikov gear system and other special gear systems for high load carrying capacity. VDI-Berichte, 47, p. 5-12.
[53]	Tsay, C.-B., Wong, Z.H., Tao, S. (1989). The mathematical model of Wildhaber-Novikov gears applicable to finite element stress analysis. Mathematical and Computer Modelling, vol. 12, no. 8 p. 939-946, D0I:10.1016/0895-7177(89)90199-4.
[54]	Evans, H.P., Snidle, R.W. (1993). Wildhaber-Novikov circular arc gears: Elastohydrodinamics. Journal of Tribology, vol. 115, no. 3, p. 487-492, D0I:10.1115/1.2921664.
[55]	Markowski, T., Batsch, M. (2015). Tooth contact analysis of Novikov convexo-concave gears. Advances in Manufacturing Science and Technology, vol. 39, no. 1, p. 53-69, D0I:10.2478/ amst-2015-0004.
[56]	Batsch, M., Markowski, T., Legutko, S., and Krolczyk, G.M. (2018). Measurement and mathematical model of convexo-concave Novikov gear mesh. Measurement, vol. 125, p. 516525, D0I:10.1016/j.measurement.2018.04.095.
[57]	Batsch, M. (2016). Comparative fatigue testing of gears with involute and convexo-concave teeth profiles. Advances in Manufacturing Science and Technology, vol. 40, no. 2, p. 5-25, D0I:10.2478/amst-2016-0007.
[58]	Batsch, M. (2020). Mathematical model and tooth contact analysis of convexo-concave helical bevel Novikov gear mesh.
389
Okorn, I. - Nagode, M. - Klemenc, J.
Strojniski vestnik - Journal of Mechanical Engineering 67(2021)5, 256-271
Mechanism and Machine Theory, vol. 149, art. ID. 103842, D0I:10.1016/j.mechmachtheory.2020.103842.
[59]	Litvin, F.L., Feng, P.H., Lagutln, S.A. (2000). Computerised generation and simulation of meshing and contact analysis of new type of Novikov-Wildhaber helical gears. NASA ARL-CR-428.
[60]	Litvin, F.L., Lu, J. (1993). Computerised simulation of generation, meshing and contact of double circular-arc helical gears. Mathematical and Computer Modelling, vol. 18, no. 5, p. 31-47, D0I:10.1016/0895-7177(93)90131-H.
[61]	Yang, S.-C. (2009). Mathematical model of a stepped triple circular-arc gear. Mechanism and Machine Theory, vol. 44, no. 5, p. 1019-1031, D0I:10.1016/j. mechmachtheory.2008.05.013.
[62]	Ariga, Y., Nagata, S. (1985). Load capacity of a new W-N gear with a basic rack of combined circular and involute profile. Journal of Mechanical Design, vol. 107, no. 4, p. 565-572, D0I:10.1115/1.3260764.
[63]	Hlebanja, G., Hlebanja, J. (2009). Uniform power transmission gears. Strojniški vestnik - Journal of Mechanical Engineering, vol. 55, no. 7-8, p. 472-483.
[64]	Liang, D., Chen, B., Tan, R., Liao, R. (2015). Geometric design and analysis of gear transmission with double circular arc-involute tooth. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 231, no. 11, p. 2100-2109, D0I:10.1177/0954406215625678.
[65]	Chen, Z., Ding, H., Li, B., Lou, L., Zhang, L., Yang, Y. (2017). Geometry and parameter design of novel circular arc helical gears for parallel-axis transmission. Advances in Mechanical Engineering, vol. 9, no. 2 p. 1-11, D0I:10.1177/1687814017690957.
[66]	Chen, Z., Zeng, M., Fuentes-Aznar A. (2020). Geometric design, meshing simulation and stress analysis of pure rolling rack and pinion mechanisms. Journal of Mechanical Design, vol. 142, no. 3, art. ID 031122, D0I:10.1115/1.4044954.
[67]	Chen, Z., Zeng, M., Fuentes-Aznar, A. (2020). Geometric design, meshing simulation and stress analysis of pure rolling cylindrical helical gear drives. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 234, no. 15 p. 3102-3115, D0I:10.1177/0954406220912265.
[68]	Berlinger, B.E., Coulbourne, J.R. (2012). Convoloid gearing technology. Gear Solutions, p. 50-57.
[69]	Hlebanja, J. (1976). Konkav-konvexe Verzahnung: Ermitlung der Zahnflanken und einige Grenzfälle. Antriebstechnik, vol. 15, no. 6, p. 324-329.
[70]	Hlebanja, J., Hlebanja, G., Umberger, M. (2020). S-gear design rules. Ventil, vol. 26, no. 4, p. 254-263.
[71]	Hlebanja, J. (1979). Betriebserfahrung an Kammwalzgetrieben mit Sonderverzahnungen. VDI-Berichte, 332, p. 283-288.
[72]	Hlebanja, J., Hlebanja, G. (2010). Spur gears with a curved path of contact for small gearing dimensions. VDI-Berichte, 2108, p. 1281-1294.
[73]	DIN 3962-1:1978-08 (1978). Tolerances for Cylindrical Gear Teeth; Tolerances for Deviations of Individual Parameters. Deutsches Institut für Normung, Berlin.
[74]	ISO 14635:2006. Gears - FZG test procedures - part 1: FZG test metod A/8,3/90 for relative scuffing load-carryng capacity of oils. International Organization for Standardization, Geneva.
[75]	Hlebanja, J., Okorn, I. (1996). Investigation of tooth surface-durability of non-involute spur gears. VDI-Berichte, 1230, p. 443-450.
[76]	Hlebanja, G., Hlebanja, J., Okorn, I. (2000). Research of gears with progressive path of contact. Proceedings of DETC'00; International Power Transmission and Gearing Conference, DETC2000/PTG-14384, p. 1-7, D0I:10.1115/detc2000/ptg-14384.
[77]	Hlebanja, J., Hlebanja, G. (2005). Anwendbarkeit der S-Verzahnung im Getriebebau. Antriebstechnick, vol. 44, p. 34-38.
[78]	Hlebanja, J., Hlebanja, G. (2005). Tooth flank durability of internal S-gears. VDI-Berichte, 1904, p. 385-393.
[79]	Hlebanja, G., Hlebanja, J. (2013). Influence of axis distance variation on rotation transmission in S-gears. Proceedings of the International Conference on Gears, p. 669-679.
[80]	Kulovec, S., Duhovnik, J. (2013). Variation of S-gear shape and the influence of the main parameters. VDI-Berichte, 2199, p. 1535-1541.
[81]	Hlebanja, G., Hlebanja, J., Čarman, M. (2009). Cylindrical wormgearings with progressively curved shape of teeth flanks. Strojniški vestnik - Journal of Mechanical Engineering, vol. 55, no. 1, p. 5-14.
[82]	Duhovnik, J., Zorko, D., Sedej, L. (2016). The effect of the teeth profile shape on polymer gear pair properties. Technical Gazette, vol. 23, no. 1, p. 199-207, D0I:10.17559/TV-20151028072528.
[83]	Zorko, D., Kulovec, S., Tavčar, J., Duhovnik, J. (2017). Different teeth profile shapes of polymer gears and comparision of their performances. Journal of Advanced Mechanical Design, System and Manufacturing, vol. 11, no. 6, p. 1-10. D0I:10.1299/jamdsm.2017jamdsm0083.
[84]	Hlebanja, G., Kulovec, S., Zorko, D., Hlebanja, J., Duhovnik, J. (2017). Influence of the tooth flank shape on thermal load of the gear. Proceedings of the International Conference on High Performance Plastic Gears, p. 1583-1592, D0I:10.51202/9783181022948-1583.
[85]	Trobentar, B., Kulovec, S., Hlebanja, G., Glodež, S. (2020). Experimental failure analysis of S-polimer gears. Engineering Failure Analysis, vol. 111, no. 4, art. ID 104496. D0I:10.1016/j.engfailanal.2020.104496.
[86]	Hlebanja, G., Kulovec, S. (2018). Thermal properties of S-gears in comparison with involute gears. IOP Conference Series: Materials Science and Engineering, vol. 393, art. ID. 012067, D0I:10.1088/1757-899X/393/1/012067.
[87]	Wadagaonkar, S., Shinde, S. (2015). Cosine gear stress analysis with experimental validation and comparison with involute gear. International Journal of Innovative Science, Engineering & Technology, vol. 2, no. 3, p. 438-443.
[88]	Liu, L., Meng, F., Ni, J. (2019). A novel non-involute gear designed based on control of relative curvature. Mechanism and Machine Theory, vol. 140, p. 144-158, D0I:10.1016/j. mechmachtheory.2019.05.022.
Operating Performance of External Non-Involute Spur and Helical Gears: A Review
390
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5
Vsebina
Vsebina
Strojniški vestnik - Journal of Mechanical Engineering letnik 67, (2021), številka 5 Ljubljana, maj 2021 ISSN 0039-2480
Izhaja mesečno
Razširjeni povzetki (extended abstracts)
Andraž Maček, Janez Urevc, Miroslav Halilovič: Prepoznava oblike planih preizkušancev z uporabo optičnih metod merjenja celotnih polj deformacij in poravnava oblike z uporabo minimizacije odstopanja preslikave	SI 29
Abdulbasit Mohammed, Hirpa G. Lemu, Belete Sirahbizu: Določitev optimalne oblike lopatic za sistem
črpanja vode na vetrno energijo na izbranih lokacijah	SI 30
Mingxing Han, Yinshui Liu, Yitao Liao, Shucai Wang: Raziskava modela in dinamičnih lastnosti
novega proporcionalnega hidravličnega ventila, ki ga poganja zvočniška tuljava	SI 31
Mohsen Soori, Mohammed Asmael: Zmanjšanje preostalih napetosti in napak odklona z virtualnim
petosnim rezkanjem turbinskih lopatic	SI 32
Mario Santos, Jaime Santos: Slabitev ultrazvočnih valov v nodularni litini z razpršitvijo: eksperimentalna
študija in simulacija	SI 33
Ivan Okorn, Marko Nagode, Jernej Klemenc: Obratovalne lastnosti zunanjih neevolventnih valjastih
zobnikov z ravnimi in poševnimi zobmi: pregled	SI 34
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, SI 29
© 2021 Strojniški vestnik. Vse pravice pridržane.
Prejeto v recenzijo: 2021-01-26
Prejeto popravljeno: 2021-04-14
Odobreno za objavo: 2021-04-23
Prepoznava oblike planih preizkušancev z uporabo optičnih metod merjenja celotnih polj deformacij in poravnava oblike z uporabo minimizacije odstopanja preslikave
Andraž Maček - Janez Urevc - Miroslav Halilovič*
Univerza v Ljubljani, Fakulteta za strojništvo, Slovenija
Raziskava predstavlja razvito metodologijo za prostorsko poravnavo oblike merjenega planega preizkušanca z obliko računskega modela. Prostorska poravnava predstavlja ključen del najnovejših postopkov inverzne identifikacije parametrov konstitutivnih modelov trdnin, ki temeljijo na meritvah celotnega polja deformacij preizkušancev s heterogenim odzivom. Poravnava je namreč nujno potrebna, saj moramo med postopkom inverzne identifikacije primerjati izmerjene in izračunane odzive preizkušancev v istih (materialnih) točkah.
Razvita metodologija obsega tako prepoznavo oblike merjenega preizkušanca iz fotografij merilnega sistema korelacije digitalnih slik (DIC) kot tudi postopek prostorske poravnave. S praktičnega vidika je to ključnega pomena, saj moramo prepoznati in poravnati obliko preizkušanca, ko je ta že na lokaciji in pripravljen za meritve deformacij. Torej ima preizkušanec v tem stanju že nanešen naključni vzorec, ki je potreben za izvedbo meritev deformacij z merilnim sistemom DIC in je vpet v trgalni stroj. Razvita metoda omogoča prostorsko poravnavo meritev, kjer je preizkušanec opazovan bodisi z eno kamero (ravninsko 2D merjenje), bodisi z več. Prepoznava oblike merjenega preizkušanca temelji na primerjavi dveh fotografij preizkušanca. Prva fotografija je zajeta na klasičen način, tako kot fotografiramo preizkušanec za potrebe postopka DIC. Pred zajemom druge fotografije pa dodatno osvetlimo preizkušanec še iz zadnje smeri. Tako zajamemo dve fotografiji, ki se razlikujeta v osvetlitvi ozadja preizkušanca. S prepoznavo slikovnih točk, katerih sivinska vrednost se med slikama ne spremeni veliko, določimo obliko merjenega preizkušanca. Prostorska poravnava te oblike z obliko preizkušanca uporabljenega v izračunu temelji na preslikavi pozicije slikovnih točk na zajetih fotografijah v prostor, kjer se nahaja oblika računalniškega modela preizkušanca. Ta preslikava je pravzaprav sestavljena iz dveh preslikav. Prva predstavlja umerjanje merilnega sistema in preslika pozicijo slikovnih točk (vrstica in stolpec na zajeti fotografiji) v pozicijo koordinatnega sistema preizkušanca (X in Y os ležita v ravnini preizkušanca in predstavljata realne mere, npr. milimeter). Druga preslikava pa nadalje preslika XY pozicije v prostor računalniškega modela preizkušanca (xy). Ta v osnovi predstavlja le transformacijo togega telesa, torej premik in zasuk točk v ravnini preizkušanca. Tako lahko preslikamo lokacije vseh slikovnih točk fotografij v prostor računalniškega modela preizkušanca. Točke, ki na fotografiji predstavljajo preizkušanec in istočasno v prostoru modela predstavljajo ozadje ali obratno, predstavljajo točke neujemanja. Večje kot je število točk, večje je neujemanje med oblikama preizkušancev. Sama poravnava je izvedena v obliki optimizacijskega problema, ki išče parametre omenjene transformacije togega telesa s ciljem minimizacije števila točk neujemanja.
Robustnost metode je prikazana na primeru preizkušanca, katerega oblika spominja na strukturo metamateriala, in primerjana s splošno znano iteracijsko metodo poravnave najbližjih točk (ICP). Izkaže se, da novo razvita metoda potrebuje manj iteracij za zagotovitev poravnave preizkušanca. Uspešno pa uspe izvesti poravnavo tudi pri različnih stopnjah začetnega neujemanja in ob prisotnosti geometrijskih napak ter šuma, kjer ICP lahko konvergira k napačni poravnavi preizkušanca.
Uspešnost metode poravnave je prikazana tudi na primeru realnih meritev heterogenega preizkušanca, kjer smo obravnavali primer merjenja z le eno kamero in primer merjenja s štirimi kamerami. V obeh primerih je razvita metoda zagotovila odlično poravnavo oblik preizkušanca.
Razvita metoda predstavlja veliko praktično uporabnost v vseh primerih inverzne identifikacije na osnovi meritev celotnega polja deformacij planih preizkušancev in predstavlja novost na omenjenem področju, saj avtomatsko prepozna obliko preizkušanca in na edinstven način poskrbi za poravnavo z upoštevanjem celotne oblike preizkušanca.
Ključne besede: meritve polja deformacij, korelacija digitalnih slik, prepoznava oblike preizkušanca, poravnava oblik, inverzna identifikacija, heterogeno polje deformacij
*Naslov avtorja za dopisovanje: Univerza v Ljubljani, Fakulteta za strojništvo, Aškerčeva 6, 1000 Ljubljana, Slovenija, miroslav.halilovic@fs.uni-lj.si
SI 29
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, SI 29
© 2021 Strojniški vestnik. Vse pravice pridržane.
Prejeto v recenzijo: 2021-01-26
Prejeto popravljeno: 2021-04-14
Odobreno za objavo: 2021-04-23
Določitev optimalne oblike lopatic za sistem črpanja vode na vetrno energijo na izbranih lokacijah
Abdulbasit Mohammed1 - Hirpa G. Lemu2 - Belete Sirahbizu1
^nanstveno-tehniška univerza v Adis Abebi, Etiopija 2Univerza v Stavangerju, Norveška
Pričujoči članek obravnava rezultate optimizacije oblike in študije zmogljivosti rotorskih lopatic, ki so primerne za delo v šibkem vetru. Vetrnice so namenjene za neposredni pogon vodnih črpalk na treh izbranih ruralnih lokacijah v Etiopiji.
Rotorske lopatice vetrnic za tri izbrane lokacije so bile oblikovane na podlagi aerodinamičnega profila SG6043 ter podatkov o hitrosti vetra na teh lokacijah. Izračunani so bili aerodinamični profili rotorskih lopatic, ki zagotavljajo maksimalen koeficient moči oziroma razmerje med dejansko močjo vetrnice in razpoložljivo močjo vetra na izbranih referenčnih območjih. Optimizacija lopatic s ciljem doseganja maksimalnega izkoristka energije vetra je bila opravljena z različnimi parametri, kot so oblika lopatic, tetive, razmerje med obodno hitrostjo in hitrostjo vetra, geometrije, nastavljeni koti itd. Optimizacija lopatic je potekala na osnovi teorije gibalne količine lopatic (BEM) v programski opremi QBlade, ki je bila razvita posebej za načrtovanje lopatic vetrnih turbin. Za iteriranje oblike je bila uporabljena koda v programskem paketu MATLAB®, ki je bila razvita za ta namen.
Na podlagi raziskav so bili določeni optimizirani profili lopatic za vetrne razmere na lokacijah v Abomsi, Metehari in Ziwayju na jugovzhodu Etiopije. Sistemi za črpanje vode na treh lokacijah lahko učinkovito delujejo z vetrnicami s tremi lopaticami premera 10,74 m, 7,34 m oz. 6,34 m. Največje vrednosti koeficienta moči Cp(r) znašajo 0,4512, 0,4587 in 0,4627 pri največjih vrednostih razmerja med obodno hitrostjo in hitrostjo vetra X(r) 3,0, 3,5 oz. 4,0.
Članek obravnava načrtovanje, analizo in optimizacijo rotorskih lopatic za določitev optimalnih parametrov zmogljivosti pri specifičnih lokalnih vetrnih razmerah, vključno z velikostjo rotorja in profili. Rezultate numerične analize bo treba še dodatno preveriti in validirati v okviru teh raziskav. Na podlagi konstrukcijskih parametrov bo razvit model celotne vetrne turbine in zgrajene bodo prototipne turbine za izvajanje preskusov na treh lokacijah.
Glavni prispevek dela je v praktičnih možnostih izvedbe. Razvite metode projektiranja rotorskih lopatic za učinkovito izkoriščanje vetrnih zmogljivosti za pogon sistemov za črpanje vode, zlasti na območjih z nizko hitrostjo vetra, imajo velik pomen za izboljšanje kakovosti življenja na ruralnih območjih v Etiopiji, ki nimajo dostopa do električne energije za oskrbo s pitno vodo. To je zelo pomembno za velik delež prebivalstva te države. Ključne besede: metoda robnih elementov, dimenzioniranje vetrnic, vetrna energija, hidravlična moč, koeficient moči, sistem za črpanje vode
SI 30
*Naslov avtorja za dopisovanje: Univerza v Stavangerju, N-4036, Stavanger, Norveška, Hirpa.g.lemu@uis.no
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, SI 29
© 2021 Strojniški vestnik. Vse pravice pridržane.
Prejeto v recenzijo: 2021-01-26
Prejeto popravljeno: 2021-04-14
Odobreno za objavo: 2021-04-23
Raziskava modela in dinamičnih lastnosti novega proporcionalnega hidravličnega ventila, ki ga poganja zvočniška tuljava
Mingxing Han1* - Yinshui Liu2 - Yitao Liao1 - Shucai Wang1 1 Kmetijska univerza Huazhong, Kolidž za inženiring, Kitajska 2 Znanstveno-tehniška univerza Huazhong, Šola za tehniške vede in strojništvo, Kitajska
Proporcionalni hidravlični ventil kot glavna krmilna komponenta v vodni hidravliki ima pomembno vlogo pri delovanju krmiljenja sistemov.
Konstruiranje visokozmogljivih vodnih proporcionalnih hidravličnih ventilov je težavna naloga zaradi nizke viskoznosti in oksidativnih lastnosti medija. V procesu konstruiranja je treba razrešiti veliko težav, kot sta npr. slabo mazanje in visoko trenje. Upor na krmilnem batu zaradi nizke viskoznosti in slabe mazalnosti vode ni le velik, ampak tudi nelinearen. Delovanje oljnih razvodnih ventilov drsniškega tipa je idealno, medtem ko se pri krmilnih batih podobnih vodnih ventilov v pogojih vodnega mazanja zaradi neoptimalne kombinacije materialov (jeklo/jeklo) pojavljajo oscilacije. Za visoko dinamično zmogljivost je zato potrebna optimizacija zgradbe hidravličnega proporcionalnega ventila in elektromehanskega aktuatorja.
Pričujoča študija predstavlja nov hitroodzivni vodni proporcionalni hidravlični ventil. Aktuator ventila je izveden kot zvočniška tuljava (VCM), katere prednosti so v hitrem odzivu, visoki natančnosti krmiljenja in majhnih izmerah. Za hitro krmiljenje ventila je bil zasnovan ročični ojačevalnik, ki zagotavlja zadostno silo za nastavljanje. Razvit je bil podroben in natančen nelinearen matematični model ventila, ki upošteva parametre konstrukcije ventila in elektromagnetne lastnosti VCM. Izvedene so bile podrobne simulacije delovanja, ki so vključevale elektromagnetne simulacije, analizo lastnosti ročičnega ojačevalnika in simulacijo dinamične zmogljivosti ventila.
Rezultati simulacij kažejo, da znaša nastavitveni čas približno 28 ms, največji prenihaj pa približno 5 %. Čas odziva na stopnico znaša približno 15 ms. Zgrajeni sta bili tudi preizkuševališči za ventil in za VCM. Rezultati preskusov prototipa kažejo, da je optimalno območje hoda VCM med 4 in 15 mm. Največji prenihaj ventila znaša okrog 10 %, nastavitveni čas pa je približno 30 ms pri odpiranju in 35 ms pri zapiranju. Rezultati preskusa dokazujejo, da ima ventil dobro statično in dinamično zmogljivost krmiljenja. Trenutno še ni primerjav med oljnimi in vodnimi hidravličnimi ventili glede vpliva pomanjkljivega mazanja na dinamični odziv na stopnico. To vprašanje je zelo zanimivo in bo navdih za podrobnejše raziskave, ki jih še načrtujejo avtorji.
V študiji je podan predlog vodnega proporcionalnega hidravličnega ventila s hitrim odzivom. Prednosti VCM, ki sta izkoriščeni pri elektromehanskem pretvorniku ventila, sta visoka hitrost in natančnost krmiljenja. Nizka viskoznost in mazalnost vodnega medija zahtevata veliko silo za premikanje krmilnega bata. Za ojačevanje sile VCM je bil uporabljen ročični ojačevalnik. Dinamična zmogljivost in odpornost ventila proti motnjam sta bili izboljšani z regulacijo položaja v zaprti zanki s povratno zvezo. Rezultati testov dokazujejo, da ventil ohranja visoko odzivnost pri različnih tlakih. Pri novem hidravličnem ventilu z VCM je bistveno izboljšana dinamična zmogljivost v primerjavi s tradicionalnimi vodnimi proporcionalnimi ventili. VCM je tako lahko dobra zamenjava za elektromagnete pri tovrstnih ventilih.
Ključne besede: proporcionalni krmilni ventil, zvočniška tuljava, optimizacija, dinamična zmogljivost
*Naslov avtorja za dopisovanje: Kmetijska univerza Huazhong, Kolidž za inženiring, Wuhan, Kitajska, 430070, hanmingxing@mail.hzau.edu.cn
SI 31
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, SI 29
© 2021 Strojniški vestnik. Vse pravice pridržane.
Prejeto v recenzijo: 2021-01-26
Prejeto popravljeno: 2021-04-14
Odobreno za objavo: 2021-04-23
Zmanjšanje preostalih napetosti in napak odklona z virtualnim petosnim rezkanjem turbinskih lopatic
Mohsen Soori - Mohammed Asmael
Univerza vzhodnega Mediterana, Oddelek za strojništvo, Turčija
Preostale napetosti v delih, izdelanih s postopki odrezavanja, so posledica trenja, oblikovanja odrezkov in nastanka toplote v območju odrezavanja. Sile in povišane temperature med odrezavanjem lahko na obdelanih turbinskih lopaticah povzročijo napake odklona, ki pa morajo biti zaradi zahtevane natančnosti izdelka čim manjše. Cilj pričujoče študije je bil izboljšanje natančnosti in zanesljivosti obdelanih turbinskih lopatic z minimizacijo preostalih napetosti in napak odklona v virtualnem obdelovalnem sistemu.
Z uporabo takega sistema je mogoče napovedati in zmanjšati preostale napetosti in napake odklona po petosnem rezkanju lopatic. Z analizo po metodi končnih elementov (MKE) so bile določene rezalne sile in temperature za vsak položaj rezalnega orodja na poti obdelave. Za optimizacijo parametrov obdelave z vidika preostalih napetosti in napak odklona je bil nato uporabljen genetski algoritem.
Virtualni obdelovalni sistem je izračunal delovno temperaturo in rezalne sile pri danih parametrih v vsaki točki na poti obdelave. Izračunane sile in temperature so bile uvožene v analitično programsko opremo Abaqus R2016X FEM za analizo preostalih napetosti in napak odklona po obdelavi. Raziskovalno delo je bilo nato validirano s primerjavo preostalih napetosti, določenih eksperimentalno po metodi rentgenske difrakcije, in rezultatov virtualnega obdelovalnega sistema. Napake odklona so bile določene z meritvijo obdelanih lopatic na koordinatnem merilnem stroju.
Preiskava fizično in virtualno obdelanih lopatic je pokazala 24,2-odstotno oz. 23,9-odstotno zmanjšanje preostalih napetosti po uporabi optimiziranih rezalnih parametrov. Zmanjšanje napake odklona je bilo 26,3-odstotno oz. 24,1-odstotno.
Predstavljeni pristop bo mogoče uporabiti pri obdelavi površin letalskih aerodinamičnih profilov na petosnih CNC-obdelovalnih strojih za zmanjšanje preostalih napetosti.
Analiza literature je pokazala, da napovedovanje in zmanjševanje preostalih napetosti in napak odklona pri rezkanju z virtualnimi obdelovalnimi sistemi ni dovolj raziskano. Razviti virtualni obdelovalni sistem bo tako lahko pripomogel k izboljšanju natančnosti in zanesljivosti lopatic.
Ključne besede: virtualna obdelava, preostale napetosti, napaka odklona, optimizacija parametrov, turbinska lopatica, petosni CNC-obdelovalni stroji
SI 32
*Naslov avtorja za dopisovanje: Univerza vzhodnega Mediterana, Oddelek za strojništvo, Famagusta, Severni Ciper, Mersin 10, Turčija, mohsen.soori@emu.edu.tr
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, SI 29
© 2021 Strojniški vestnik. Vse pravice pridržane.
Prejeto v recenzijo: 2021-01-26
Prejeto popravljeno: 2021-04-14
Odobreno za objavo: 2021-04-23
Slabitev ultrazvočnih valov v nodularni litini z razpršitvijo: eksperimentalna študija in simulacija
Mario Santos* - Jaime Santos
Univerza Coimbra, CEMMPRE, Oddelek za elektrotehniko in računalništvo, Portugalska
Nodularna železova litina je razširjena zaradi svoje livnosti, visoke toplotne prevodnosti in dobrih mehanskih lastnosti, predvsem natezne trdnosti in duktilnosti. Mehanske lastnosti kovin so v veliki meri odvisne od njihove mikrostrukture. Nodularna litina se proizvaja z dodajanjem manjših količin elementov tik pred strjevanjem, npr. magnezija in cerija. Ti dodatki povzročijo rast grafitnih vozlov, njihova oblika in porazdelitev pa sta odločilni za lastnosti kovine. Neporušne preiskave mikrostruktur so zelo pomembne za identifikacijo osnovne mikrostrukture in vozlov.
Ultrazvočna karakterizacija ponuja vrsto prednosti v primerjavi s porušnimi metalografskimi metodami. Interakcijo ultrazvočnih valov z mikrostrukturo je mogoče spremljati z akustičnimi parametri, kot sta hitrost in slabitev. Obstajata dva glavna mehanizma slabitve ultrazvočnih valov: absorpcija in razpršitev. Absorpcija je povezana s toplotnimi izgubami, histerezo in viskoznimi izgubami. Do razpršitve prihaja zaradi heterogenosti v obliki mej med zrni, praznin, vključkov, delcev drugih faz in poroznosti. Prav ta mehanizem slabitve je najbolj pomemben pri heterogenih materialih, kot je železova litina.
V pričujočem delu je bila eksperimentalno in s simulacijami ocenjena slabitev ultrazvočnih valov zaradi razpršitve na strukturah s kompleksno porazdelitvijo, ki je značilna za železovo litino. Uporabljena je bila eksperimentalna določitev slabitve po impulzni metodi. V ta namen so bili obdelani odboji ultrazvočnih impulzov na prednjih in na zadnjih ploskvah šestih preizkušancev prizmatične oblike. Preizkušanci so bili pripravljeni po standardnih metalografskih metodah in analiza optičnih mikroposnetkov je pokazala nizko prostorsko variabilnost porazdelitve velikosti vozlov. Slabitev z razpršitvijo je bila ocenjena po modelih Truella in Papadakisa, rezultati pa so bili primerjani z rezultati eksperimentov. Poudarek je bil na modelu Papadakisa, ki upošteva porazdelitev velikosti sipalcev. Rezultati obeh analitičnih modelov se močno razlikujejo od eksperimentalnih.
Za analizo vpliva velikosti in koncentracije sipalcev v preučevanih strukturah, ki sta pomembni za učinkovit razvoj modelov, so bila pri izvedbi simulacij uporabljena orodja MATLAB k-Wave. Ta simulacijska orodja uporabljajo preprost napovedovalni model na osnovi rešitve v časovni domeni po psevdospektralni metodi v k-prostoru za povezovanje akustičnih enačb prvega reda za homogene in heterogene medije v eni, dveh in treh razsežnostih. Paket k-Wave omogoča modeliranje poljubnih virov in zaznavanje površin z usmerjenimi elementi. Pri širjenju valov je mogoče upoštevati nelinearnost, heterogenosti in zakon absorpcije.
Simulator, ki je bil razvit za preizkušance, je omogočil prediktivni model, v katerem je slabitev sorazmerna s peto potenco velikosti sipalcev, formulacija pa je podobna tisti, ki jo predlagajo analitični modeli. Predhodno eksperimentalno pridobljeno znanje o slabitvi valov omogoča ustvarjanje simulacijskih modelov s spreminjanjem koncentracije sipalcev in enakovredne velikosti. Na ta način je omogočeno posnemanje struktur s kompleksnimi porazdelitvami sipalcev z vidika njihove enakomerne koncentracije in velikosti.
Ključne besede: modeliranje, anizotropija, impulzna ultrazvočna metoda, simulacija, slabitev ultrazvoka
*Naslov avtorja za dopisovanje: DEEC, Polo II, Pinhal de Marrocos, 3030-290 Coimbra, Portugalska, marioj@deec.uc.pt
SI 33
Strojniški vestnik - Journal of Mechanical Engineering 67(2021)5, SI 29
© 2021 Strojniški vestnik. Vse pravice pridržane.
Prejeto v recenzijo: 2021-01-26
Prejeto popravljeno: 2021-04-14
Odobreno za objavo: 2021-04-23
Obratovalne lastnosti zunanjih neevolventnih valjastih zobnikov z ravnimi in poševnimi zobmi: pregled
Ivan Okom - Marko Nagode - Jernej Klemenc
Univerza v Ljubljani, Fakulteta za strojništvo, Slovenija
Namen preglednega znanstvenega članka je bil podati sistematičen pregled neevolventnih zobnikov, njihovih prednosti in slabosti ter možnosti uporabe v praktičnih aplikacijah. V strokovni in znanstveni literaturi nismo našli celovitega pregleda od starih oblik, kot so cikloidno ozobje in Novikov ozobje, do ozobij, ki so bila razvita v zadnjem desetletju. Želja avtorjev je bila pripraviti kritičen pregled, ki bo uporaben tako za raziskovalce novih oblik neevolventnih zobnikov kot za uporabnike le-teh v industriji.
Večina zobnikov v industrijski praksi ima evolventno obliko bokov zob. Z vidika bočne nosilnosti, obrabe in izkoristka ta oblika ni optimalna. Ključne slabosti evolventnega ozobja so: omejeno minimalno število zob pastorka, visoka kontaktna obremenitev zaradi konveksno-konveksnega dotika bokov zob in neugodne kinematske razmere (velike drsne hitrosti) na začetku in koncu ubiranja. Številni raziskovalci so razvili in analizirali različne neevolventne oblike bokov zob, ki se niso širše uveljavile v praksi. Glavna razloga sta stroški izdelave in občutljivost na točnost izdelave in montaže.
V	preteklih desetletjih so bile raziskave večinoma usmerjene v izboljšavo lastnosti evolventnih zobnikov. Podrobno so bili raziskani vplivi geometrije, mazanja, točnosti izdelave in novih materialov na nosilnost in hrupnost ozobja. Ker v članku primerjamo evolventne zobnike z neevolventnimi, so uvodoma povzete glavne lastnosti evolventnih zobnikov. Pojasnjene so možnosti za izboljšanje njihovih obratovalnih lastnosti (uporaba profilnih premikov, profilne korekcije, nesimetrično ozobje).
Ocena lastnosti evolventnih in neevolventnih ozobij temelji na uveljavljenih kriterijih, kot so Hertzov kotalni tlak, temperatura v dotiku, debelina oljnega sloja, upogibna napetost v korenu zoba in hrupnost ozobja. V drugem poglavju so kriteriji za oceno podrobno predstavljeni. Na podlagi enačb, ki popisujejo razmere v dotiku bokov zob, so pojasnjeni trendi, ki jih raziskovalci upoštevajo pri razvoju neevolventnih zobnikov. Pri oblikovanju neevolventnih zobnikov je cilj doseči čim večji srednji krivinski radij, kar posledično pomeni zmanjšanje Herzovega tlaka, zmanjšanje temperature v dotiku in povečanje debeline oljnega sloja. Drugi ključni cilj je čim manjše relativno drsenje v dotiku, kar prispeva k manjši toplotni obremenitvi zobnikov in boljšemu izkoristku. Hrupnost zobnikov je manjša, če je stopnja prekrivanja večja. Večina neevolventnih zobnikov ima manjšo stopnjo prekrivanja, kar negativno vpliva na dinamične sile in hrupnost.
V	tretjem poglavju so opisno in s slikami predstavljene različne oblike neevolventnih zobnikov. Povzeti in komentirani so rezultati raziskav iz navedene literature. Podrobneje so predstavljeni rezultati raziskav na S-zobnikih, ki so bili razviti in preizkušani na Fakulteti za strojništvo Univerze v Ljubljani. S preizkusi S-zobnikov in primerljivih evelventnih zobnikov je bila potrjena večja bočna nosilnost in manjša toplotna obremenitev S-zobnikov. Potencial na področju S-zobnikov predstavljao zobniki s poševnimi zobmi, ki pa še niso bili eksperimentalno raziskani. V četrtem poglavju je združen pregled lastnosti predstavljenih neevolventnih zobnikov v primerjavi z evolventnimi zobniki. Uporabljeni kriteriji so bočna nosilnost, korenska nosilnost, toplotna obemenitev, hrupnost in občutljivost ozobja na napake izdelave in montaže. Večina neevolventnih zobnikov ima večjo korensko in bočno nosilnost ter manjšo toplotno obremenitev.
Za širšo uveljavitev nevolventnih zobnikov v industrijski praksi bodo potrebne nadaljnje eksperimentalne potrditve teoretičnih predpostavk. Z razvojem novih tehnologij in materialov se uporabnost neevolventnih zobnikov sicer povečuje. Ker pa je pri neevolventnih zobnikih mnogo težje doseči enak razred točnosti, so le-ti še vedno omejeni na specifične primere uporabe. Trenutno stanje razvoja in raziskav kaže, da neevolventni zobniki ne morejo izpodriniti evolventnih v večjem obsegu.
Ključne besede: neevolventni zobniki, profil boka zoba, ubirnica, S-zobniki, cikloidno ozobje, Novikov ozobje, nosilnost ozobja
SI 34
*Naslov avtorja za dopisovanje: Univerza v Ljubljani, Fakulteta za strojništvo, Aškerčeva 6, 1000 Ljubljana, Slovenija, ivan.okorn@fs.uni-lj.si
Guide for Authors
All manuscripts must be in English. Pages should be numbered sequentially. The manuscript should be composed in accordance with the Article Template given above. The maximum length of contributions is 12 pages (approx. 5000 words). Longer contributions will only be accepted if authors provide justification in a cover letter. For full instructions see the Information for Authors section on the journal's website: http://en.sv-jme.eu .
SUBMISSION:
Submission to SV-JME is made with the implicit understanding that neither the manuscript nor the essence of its content has been published previously either in whole or in part and that it is not being considered for publication elsewhere. All the listed authors should have agreed on the content and the corresponding (submitting) author is responsible for having ensured that this agreement has been reached. The acceptance of an article is based entirely on its scientific merit, as judged by peer review. Scientific articles comprising simulations only will not be accepted for publication; simulations must be accompanied by experimental results carried out to confirm or deny the accuracy of the simulation. Every manuscript submitted to the SV-JME undergoes a peer-review process.
The authors are kindly invited to submit the paper through our web site: http://ojs.sv-jme.eu. The Author is able to track the submission through the editorial process - as well as participate in the copyediting and proofreading of submissions accepted for publication - by logging in, and using the username and password provided.
SUBMISSION CONTENT:
The typical submission material consists of:
-	A manuscript (A PDF file, with title, all authors with affiliations, abstract, keywords, highlights, inserted figures and tables and references),
-	Supplementary files:
•	a manuscript in a WORD file format
•	a cover letter (please see instructions for composing the cover letter)
•	a ZIP file containing figures in high resolution in one of the graphical formats (please see instructions for preparing the figure files)
•	possible appendicies (optional), cover materials, video materials, etc. Incomplete or improperly prepared submissions will be rejected with explanatory comments provided. In this case we will kindly ask the authors to carefully read the Information for Authors and to resubmit their manuscripts taking into consideration our comments.
COVER LETTER INSTRUCTIONS:
Please add a cover letter stating the following information about the submitted paper:
1.	Paper title, list of authors and their affiliations. One corresponding author should be provided.
2.	Type of paper: original scientific paper (1.01), review scientific paper (1.02) or short scientific paper (1.03).
3.	A declaration that neither the manuscript nor the essence of its content has been published in whole or in part previously and that it is not being considered for publication elsewhere.
4.	State the value of the paper or its practical, theoretical and scientific implications. What is new in the paper with respect to the state-of-the-art in the published papers? Do not repeat the content of your abstract for this purpose.
5.	We kindly ask you to suggest at least two reviewers for your paper and give us their names, their full affiliation and contact information, and their scientific research interest. The suggested reviewers should have at least two relevant references (with an impact factor) to the scientific field concerned; they should not be from the same country as the authors and should have no close connection with the authors.
FORMAT OF THE MANUSCRIPT:
The manuscript should be composed in accordance with the Article Template. The manuscript should be written in the following format:
-	A Title that adequately describes the content of the manuscript.
-	A list of Authors and their affiliations.
-	An Abstract that should not exceed 250 words. The Abstract should state the principal objectives and the scope of the investigation, as well as the methodology employed. It should summarize the results and state the principal conclusions.
-	4 to 6 significant key words should follow the abstract to aid indexing.
-	4 to 6 highlights; a short collection of bullet points that convey the core findings and provide readers with a quick textual overview of the article. These four to six bullet points should describe the essence of the research (e.g. results or conclusions) and highlight what is distinctive about it.
-	An Introduction that should provide a review of recent literature and sufficient background information to allow the results of the article to be understood and evaluated.
-	A Methods section detailing the theoretical or experimental methods used.
-	An Experimental section that should provide details of the experimental set-up and the methods used to obtain the results.
-	A Results section that should clearly and concisely present the data, using figures and tables where appropriate.
-	A Discussion section that should describe the relationships and generalizations shown by the results and discuss the significance of the results, making comparisons with previously published work. (It may be appropriate to combine the Results and Discussion sections into a single section to improve clarity.)
-	A Conclusions section that should present one or more conclusions drawn from the results and subsequent discussion and should not duplicate the Abstract.
-	Acknowledgement (optional) of collaboration or preparation assistance may be included. Please note the source of funding for the research.
-	Nomenclature (optional). Papers with many symbols should have a nomenclature that defines all symbols with units, inserted above the references. If one is used, it must contain all the symbols used in the manuscript and the definitions should not be repeated in the text. In all cases, identify the symbols used if they are not widely recognized in the profession. Define acronyms in the text, not in the nomenclature.
-	References must be cited consecutively in the text using square brackets [1] and collected together in a reference list at the end of the manuscript.
-	Appendix(-icies) if any.
SPECIAL NOTES
Units: The SI system of units for nomenclature, symbols and abbreviations should be followed closely. Symbols for physical quantities in the text should be written in italics (e.g. v, T, n, etc.). Symbols for units that consist of letters should be in plain text (e.g. ms-1, K, min, mm, etc.). Please also see: http://physics.nist.gov/cuu/pdf/sp811.pdf .
Abbreviations should be spelt out in full on first appearance followed by the abbreviation in parentheses, e.g. variable time geometry (VTG). The meaning of symbols and units belonging to symbols should be explained in each case or cited in a nomenclature section at the end of the manuscript before the References.
Figures (figures, graphs, illustrations digital images, photographs) must be cited in consecutive numerical order in the text and referred to in both the text and the captions as Fig. 1, Fig. 2, etc. Figures should be prepared without borders and on white grounding and should be sent separately in their original formats. If a figure is composed of several parts, please mark each part with a), b), c), etc. and provide an explanation for each part in Figure caption. The caption should be self-explanatory. Letters and numbers should be readable (Arial or Times New Roman, min 6 pt with equal sizes and fonts in all figures). Graphics (submitted as supplementary files) may be exported in resolution good enough for printing (min. 300 dpi) in any common format, e.g. TIFF, BMP or JPG, PDF and should be named Fig1.jpg, Fig2.tif, etc. However, graphs and line drawings should be prepared as vector images, e.g. CDR, AI. Multi-curve graphs should have individual curves marked with a symbol or otherwise provide distinguishing differences using, for example, different thicknesses or dashing.
Tables should carry separate titles and must be numbered in consecutive numerical order in the text and referred to in both the text and the captions as Table 1, Table 2, etc. In addition to the physical quantities, such as t (in italics), the units [s] (normal text) should be added in square brackets. Tables should not duplicate data found elsewhere in the manuscript. Tables should be prepared using a table editor and not inserted as a graphic.
REFERENCES:
A reference list must be included using the following information as a guide. Only cited text references are to be included. Each reference is to be referred to in the text by a number enclosed in a square bracket (i.e. [3] or [2] to [4] for more references; do not combine more than 3 references, explain each). No reference to the author is necessary.
References must be numbered and ordered according to where they are first mentioned in the paper, not alphabetically. All references must be complete and accurate. Please add DOI code when available. Examples follow.
Journal Papers:
Surname 1, Initials, Surname 2, Initials (year). Title. Journal, volume, number, pages, DOI code.
[1]	Hackenschmidt, R., Alber-Laukant, B., Rieg, F. (2010). Simulating nonlinear materials under centrifugal forces by using intelligent cross-linked simulations. Strojniški vestnih - Journal of Mechanical Engineering, vol. 57, no. 7-8, p. 531-538, DOI:10.5545/sv-jme.2011.013.
Journal titles should not be abbreviated. Note that journal title is set in italics. Books:
Surname 1, Initials, Surname 2, Initials (year). Title. Publisher, place of publication.
[2]	Groover, M.P. (2007). Fundamentals of Modern Manufacturing. John Wiley & Sons, Hoboken.
Note that the title of the book is italicized.
Chapters in Books:
Surname 1, Initials, Surname 2, Initials (year). Chapter title. Editor(s) of book, book title. Publisher, place of publication, pages.
[3]	Carbone, G., Ceccarelli, M. (2005). Legged robotic systems. Kordič, V., Lazinica, A., Merdan, M. (Eds.), Cutting Edge Robotics. Pro literatur Verlag, Mammendorf, p. 553576.
Proceedings Papers:
Surname 1, Initials, Surname 2, Initials (year). Paper title. Proceedings title, pages.
[4]	Štefanič, N., Martinčevič-Mikič, S., Tošanovič, N. (2009). Applied lean system in process industry. MOTSP Conference Proceedings, p. 422-427.
Standards:
Standard-Code (year). Title. Organisation. Place.
[5]	ISO/DIS 16000-6.2:2002. Indoor Air - Part 6: Determination of Volatile Organic Compounds in Indoor and Chamber Air by Active Sampling on TENAX TA Sorbent, Thermal Desorption and Gas Chromatography using MSD/FID. International Organization for Standardization. Geneva.
WWW pages:
Surname, Initials or Company name. Title, from http://address, date of access.
[6]	Rockwell Automation. Arena, from http://www.arenasimulation.com, accessed on 200909-07.
EXTENDED ABSTRACT:
When the paper is accepted for publishing, the authors will be requested to send an extended abstract (approx. one A4 page or 3500 to 4000 characters or approx. 600 words). The instruction for composing the extended abstract are published on-line: http://www.sv-jme. eu/information-for-authors/ .
COPYRIGHT:
Authors submitting a manuscript do so on the understanding that the work has not been published before, is not being considered for publication elsewhere and has been read and approved by all authors. The submission of the manuscript by the authors means that the authors automatically agree to transfer copyright to SV-JME when the manuscript is accepted for publication. All accepted manuscripts must be accompanied by a Copyright Transfer Agreement, which should be sent to the editor. The work should be original work by the authors and not be published elsewhere in any language without the written consent of the publisher. The proof will be sent to the author showing the final layout of the article. Proof correction must be minimal and executed quickly. Thus it is essential that manuscripts are accurate when submitted. Authors can track the status of their accepted articles on http://en.sv-jme.eu/.
PUBLICATION FEE:
Authors will be asked to pay a publication fee for each article prior to the article appearing in the journal. However, this fee only needs to be paid after the article has been accepted for publishing. The fee is 380 EUR (for articles with maximum of 6 pages), 470 EUR (for articles with maximum of 10 pages), plus 50 EUR for each additional page. The additional cost for a color page is 90.00 EUR (only for a journal hard copy; optional upon author's request). These fees do not include tax.
Strojniški vestnik -Journal of Mechanical Engineering Aškerčeva 6, 1000 Ljubljana, Slovenia, e-mail: info@sv-jme.eu
http://www.sv-jme.eu
Contents
Papers
Andraž Maček, Janez Urevc, Miroslav HaLiLovič:
Flat Specimen Shape Recognition Based on Full-Field Optical Measurements and Registration Using Mapping Error Minimization Method
Abdulbasit Mohammed, Hirpa G. Lemu, Belete Sirahbizu:
Determining Optimum Rotary Blade Design for Wind-Powered Water-Pumping
Systems for Local Selected Sites
Mingxing Han, Yinshui Liu, Yitao Liao, Shucai Wang:
Investigation on the Modeling and Dynamic Characteristics
of a Novel Hydraulic Proportional Valve Driven by a Voice Coil Motor
Mohsen Soori, Mohammed Asmael:
Virtual Minimization of Residual Stress and Deflection Error in the Five-Axis Milling of Turbine Blades
Mario Santos, Jaime Santos:
Ultrasonic Scattering Attenuation in Nodular Cast Iron: Experimental and Simulation Studies
Ivan Okorn, Marko Nagode, Jernej Klemenc:
Operating Performance of External Non-Involute Spur and Helical Gears: A Review
203
214
223
245
9770039248001