Strojniški vestnik
Journal of Mechanical Engineering
Strojniški vestnik - Journal of Mechanical Engineering (SV-JME)
Aim and Scope
The international journal publishes original and (mini)review articles covering the concepts of materials science, mechanics, kinematics, thermodynamics, energy and environment, mechatronics and robotics, fluid mechanics, tribology, cybernetics, industrial engineering and structural analysis.
The journal follows new trends and progress proven practice in the mechanical engineering and also in the closely related sciences as are electrical, civil and process engineering, medicine, microbiology, ecology, agriculture, transport systems, aviation, and others, thus creating a unique forum for interdisciplinary or multidisciplinary dialogue.
The international conferences selected papers are welcome for publishing as a special issue of SV-JME with invited co-editor(s).
Editor in Chief
Vincenc Butala
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
Technical Editor Pika Škraba
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
Founding Editor
Bojan Kraut
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Editorial Office
University of Ljubljana, Faculty of Mechanical Engineering SV-JME, Aškerčeva 6, SI-1000 Ljubljana, Slovenia Phone: 386 (0)1 4771 137 Fax: 386 (0)1 2518 567 info@sv-jme.eu, http://www. sv-jme.eu
Print: Grafex, d.o.o., printed in 310 copies Founders and Publishers
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
University of Maribor, Faculty of Mechanical Engineering, Slovenia
Association of Mechanical Engineers of Slovenia Chamber of Commerce and Industry of Slovenia, Metal Processing Industry Association
President of Publishing Council
Branko Širok
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
International Editorial Board
Kamil Arslan, Karabuk University, Turkey
Hafiz Muhammad Ali, University of Engineering and Technology, Pakistan Josep M. Bergada, Politechnical University of Catalonia, Spain Anton Bergant, Litostroj Power, Slovenia
Miha Boltežar, UL, Faculty of Mechanical Engineering, Slovenia
Franci Čuš, UM, Faculty of Mechanical Engineering, Slovenia
Anselmo Eduardo Diniz, State University of Campinas, Brazil
Igor Emri, UL, Faculty of Mechanical Engineering, Slovenia
Imre Felde, Obuda University, Faculty of Informatics, Hungary
Janez Grum, UL, Faculty of Mechanical Engineering, Slovenia
Imre Horvath, Delft University of Technology, The Netherlands
Aleš Hribernik, UM, Faculty of Mechanical Engineering, Slovenia
Soichi Ibaraki, Kyoto University, Department of Micro Eng., Japan
Julius Kaplunov, Brunel University, West London, UK
Iyas Khader, Fraunhofer Institute for Mechanics of Materials, Germany
Jernej Klemenc, UL, Faculty of Mechanical Engineering, Slovenia
Milan Kljajin, J.J. Strossmayer University of Osijek, Croatia
Peter Krajnik, Chalmers University of Technology, Sweden
Janez Kušar, UL, Faculty of Mechanical Engineering, Slovenia
Gorazd Lojen, UM, Faculty of Mechanical Engineering, Slovenia
Thomas Lübben, University of Bremen, Germany
Janez Možina, UL, Faculty of Mechanical Engineering, Slovenia
George K. Nikas, KADMOS Engineering, UK
José L. Ocana, Technical University of Madrid, Spain
Miroslav Plančak, University of Novi Sad, Serbia
Vladimir Popović, University of Belgrade, Faculty of Mech. Eng., Serbia Franci Pušavec, UL, Faculty of Mechanical Engineering, Slovenia Bernd Sauer, University of Kaiserlautern, Germany Rudolph J. Scavuzzo, University of Akron, USA Arkady Voloshin, Lehigh University, Bethlehem, USA
Vice-President of Publishing Council
Jože Balič
University of Maribor, Faculty of Mechanical Engineering, Slovenia
Strojniški vestnik Journal of Mechanical "''i^' Engineering
Cover: The cover shows a system for adaptive robotic deburring of die-cast parts. On the upper part of the image, the reference workpiece fixed on the robot arm and laser profilometer are shown during 3D surface measuring. On the lower part, a workpiece with a sketched coordinate system is presented. Orange arrows show the burr to be removed. On the left side, the 3D measurement is shown and on the right, displacements after registration between the reference and currently processed workpieces. Image Courtesy:
Marko Leber, MKreativ, www.mkreativ.si Urban Pavlovčič, University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
ISSN 0039-2480
General information
Strojniški vestnik - Journal of Mechanical Engineering is published in 11 issues per year (July and August is a double issue).
Institutional prices include print & online access: institutional subscription price and foreign subscription €100,00 (the price of a single issue is €10,00); general public subscription and student subscription €50,00 (the price of a single issue is €5,00). Prices are exclusive of tax. Delivery is included in the price. The recipient is responsible for paying any import duties or taxes. Legal title passes to the customer on dispatch by our distributor.
Single issues from current and recent volumes are available at the current single-issue price. To order the journal, please complete the form on our website. For submissions, subscriptions and all other information please visit: http://en.sv-jme.eu/.
You can advertise on the inner and outer side of the back cover of the journal. The authors of the published papers are invited to send photos or pictures with short explanation for cover content.
We would like to thank the reviewers who have taken part in the peerreview process.
© 2015 Strojniški vestnik - Journal of Mechanical Engineering. All rights reserved. SV-JME is indexed / abstracted in: SCI-Expanded, Compendex, Inspec, ProQuest-CSA, SCOPUS, TEMA. The list of the remaining bases, in which SV-JME is indexed, is available on the website.
The journal is subsidized by Slovenian Research Agency.
Strojniški vestnik - Journal of Mechanical Engineering is available on http://www.sv-jme.eu, where you access also to papers' supplements, such as simulations, etc.
Contents
Strojniški vestnik - Journal of Mechanical Engineering volume 62, (2016), number 4 Ljubljana, April 2016 ISSN 0039-2480
Published monthly
Papers
Hubert Kosler, Urban Pavlovčič, Matija Jezeršek, Janez Možina: Adaptive Robotic Deburring of Die-
Cast Parts with Position and Orientation Measurements Using a 3D Laser-Triangulation Sensor 207 Shijun Ji, Huijuan Yu, Ji Zhao, Xiaolong Liu, Mingxu Zhao: Ultra-Precision Machining of a Large
Amplitude Sinusoidal Ring Surface Based on a Slow Tool Servo	213
Andrej Škrlec, Jernej Klemenc: Estimating the Strain-Rate-Dependent Parameters of the Cowper-
Symonds and Johnson-Cook Material Models using Taguchi Arrays	220
Dan Ni, Minguan Yang, Bo Gao, Ning Zhang, Zhong Li: Flow Unsteadiness and Pressure Pulsations in
a Nuclear Reactor Coolant Pump	231
Deniz Duran, Celalettin Karadogan: Determination of Coulomb's Friction Coefficient Directly from
Cylinder Compression Tests	243
Ahmet Selim Dalkilig, Ali Celen, Alican Cebi, Somchai Wongwises: Effect of Refrigerant Type and
Insulation Thickness on Refrigeration Systems of Land and Sea Vehicles	252
Strojniški vestnik - Journal of Mechanical Engineering 62(2016)4, 207-212 © 2016 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2015.3227	Original Scientific Paper
Adaptive Robotic Deburring of Die-Cast Parts with Position and Orientation Measurements Using a 3D Laser-Triangulation Sensor
Hubert Kosler1 - Urban Pavlovčič2* - Matija Jezeršek2 - Janez Možina2 1 Yaskawa Motoman, Slovenia 2 University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
A system for adaptive robotic deburring with a correction of the errors in workpiece positioning is presented. The correction is based on 3D measurements of the workpiece's surface and its registration to the target surface, measured on a reference workpiece. The surface measurement is performed with a laser-triangulation profilometer. The reference tool path is determined using robot teaching on the reference, already deburred, workpiece. The positioning errors of the currently processed workpiece are compensated by tool path adaptation in accordance with the registration results by means of rotation and translation. The experiments showed that the average precision of localization is 0.06 mm and the average bias between the true and measured values is 0.23 mm. The developed adaptive system is also applicable in other similar applications where it is difficult to ensure repeatable clamping of the workpiece. Keywords: adaptive robotic machining, deburring, laser triangulation, position error correction, localization
Highlights
•	A new, simplified method for workpiece position and orientation corrections is presented.
•	The 3D measurement of the workpiece is based on a laser-triangulation technique.
•	The differences in position and orientation are calculated using an iterative closest point algorithm.
•	The tool path of the current workpiece is adapted by rotating and translating the reference tool path in accordance with the measured positioning error.
•	The method is especially applicable for manufacturing lines with large batch sizes, such as the robotic deburring of die-cased parts.
0 INTRODUCTION
Burr formation (see Fig. 1) is a major concern in the surface and edge finishing of workpieces, since it may injure a human worker and reduce the quality of the workpieces, [1] and [2]. To remove the burr, another process must be introduced into a manufacturing line, i.e., deburring. Unfortunately, manual deburring is time consuming, costly and demands a very high level of skill and experience to maintain consistency [3]. However, it is still common, even in today's most fully automated factories [4]. The process itself adds little or no added value to the product, but the costs for some parts can be as high as 35 % of total cost of the part, [5] and [6]. Therefore, the need to automate the process is obvious.
Most attempts at deburring automation are based on the use of robots to manipulate the workpiece or deburring tool along a predefined path. Usually, they are programmed manually by the robot operators. This process is often referred to as robot teaching. The conventional robot-deburring approach is based on the assumptions that the workpiece has no defects and is located at a known position [7]. This is why the robot can travel along a rigid programmed path
[8]. Unfortunately, those assumptions are not always true: die-casted workpieces may vary in geometry to a certain level, and position and orientation errors occur when the workpiece is fixed in the jig or grasped by the robot. While we can realistically assume that workpieces with defects are removed from the production line in steps prior to deburring, differences in the orientation and the position between grasped workpieces are inevitable. Minimal misalignments can be neglected when active force control is introduced, [5] and [9]. However, with increasing misalignments the accuracy of the deburring decreases and, consequently, too much or too little material is removed. That is why minimizing these differences prior to deburring is helpful [10].
This issue is traditionally solved by either specially designed fixtures that ensure a repeatable workpiece position or by measuring the position and orientation of the workpiece's key features [7]. The latter approach is more suitable when the batch size is large. Several authors developed systems to correct these misalignments in applications for robot deburring or other robotized processes. Song et al. [10] presented an approach where the tool path is generated based on a CAD model. Then the taught
points are matched to a generated tool path, so that the modified tool path is acquired. Habibi and Pescaru [11] registered a patent for a system that trains a robot to recognize and localize objects using 3D vision. Biegelbauer and Vincze [12] 3D measured the surface, segmented the range image and fitted a cylinder to detect the actual position of the bore. A system that scans and localizes the workpieces in 3D for assembly and pick-and-place operations was presented by Skotheim et al. [13]. Rajaraman et al. [14] developed a laser-scanner-based localization system for automatic robot welding.
In contrast to the presented approaches, our method preserves the robot's teaching step which is made on the reference workpiece. During normal operation the workpiece misalignment relative to the reference is measured and the robot tool path is properly adapted. A custom laser-triangulation sensor for a harsh industrial environment was developed for 3D measurements of the surface of a currently processed workpiece. The positional misalignment is then determined using a registration algorithm that calculates the misalignment relative to the reference workpiece. This approach requires a minimal adjustment of the system and only one additional software module.
1 EXPERIMENTAL SETUP
The deburring of a die-cast part is presented in Fig. 1. The burr (indicated by orange arrows) is located on a contact of the upper and lower parts of the mold (the X-Z plane of the workpiece). The die-cast part is processed with the system schematically presented in Fig. 2. The robot first grasps it at the bottom side. Then it makes a 3D scan of the side that is visible in Fig. 1. For that purpose the robot moves the part under the laser-triangulation sensor that is fixed in space relative to the floor. The 3D data is then used to calculate the position and orientation of the current workpiece relative to the reference one (see the second paragraph for details). Finally, the part is deburred with the deburring tool, which is also fixed relative to the floor.
A Yaskawa Motoman MA1800 robot [15] with a DX100 controller is used for the workpiece manipulation. It is a vertical jointed-arm able to handle a payload of 15 kg within 3.2 m of vertical and 1.8 m of horizontal reach, and with a repeatability of 0.08 mm.
The 3D measurements are based on the lasertriangulation principle [16], which is used in a custom-developed laser-triangulation sensor (Yaskawa
MOTOSense profilometer). The sensor consists of a laser line projector (Flexpoint MVnano, wavelength 660 nm, power 100 mW, line thickness 0.1 mm) and a CCD camera (Basler, Ace acA645-100gm, 659x494 pix, 100 FPS), which are attached into a housing. The aluminum housing protects the projector and the camera from dust and any potential collision. Its dimensions and the measuring range are shown in Fig. 3. The precision of the measured profile is 0.1 mm and the typical scanning resolution is 0.2 mm, which means that the measurement time for a 20-mm-long surface is 1 s.
Fig. 1. Die-cast part used in the experiments. The location of the burr is indicated with orange arrows and is approximately the X-Z sectioning plane
3D measuring phase
MOTOsense
deburring tool
deburring phase
_ robot )__<_
controller
Fig. 2. Schematic of the system with the 3D measuring and deburring stages
The images acquired by the triangulation sensor are processed with a PC. In the first step the profiles
are detected from each image with sub-pixel accuracy. In the next step the series of profiles is transformed into a 3D surface according to the known positions of the camera, the laser projector and the robot arm. Finally, the workpiece misalignment is calculated and the correction is sent to the robot controller.
I_45_U
Fig. 3. Dimensions of the laser-triangulation sensor and its measuring range
2 POSITION AND ORIENTATION MEASUREMENT
The position and orientation measurement method is based on an assumption that the shape of each part does not change significantly, i.e., it is expected to be within the allowable tolerances. Thus, only the misalignment of the current workpiece with respect to the reference workpiece is determined and the transformation of the taught tool path is calculated. A flowchart of the method is presented in Fig. 4. The robot teaching is performed on one workpiece that is manually deburred prior to the teaching. Hereinafter, it is referred to as the reference workpiece, but note that it is not a special workpiece, since any workpiece from the same batch can serve as a reference. The surface of the reference workpiece is measured in 3D
after the robot teaching and, which is most important, during the same clamping.
When other workpieces are about to be deburred with the robot system, their surface is first measured in 3D. Then the relative rotation and translation with respect to the reference workpiece are calculated based on a registration using the iterative closest point (ICP) algorithm [17], implemented as an ICP function in trimesh2 SDK [18]. The algorithm minimizes the distance between the target (the surface of the reference workpiece) and the source (the surface of the currently processed workpiece) sets of points. The registration is done by translating and rotating the latter. The results of the algorithm are a transformation matrix (translation and rotation) and a goodness-of-fit estimator.
Fig. 4. Flowchart of the adaptive deburring process
The tool path taught on the reference workpiece is then transformed to match the position and orientation of the currently processed workpiece using the transformation matrix.
The central region of the workpiece was selected for 3D measurement, as it is shown in Fig. 5. The geometry of this region confines all translations and rotations, which is the required condition for reliable and stable registration [17].
3 EXPERIMENTS
The version of the ICP algorithm used in the experiments uses four parameters: the maximum distance between the corresponding points (MPD), the maximum number of iterations (MNI), and two values to determine whether the convergence has been reached. The third parameter determines the number of the last N iterations (NNE), where the error is increased in order to finish the algorithm. Value of N is determined by the fourth parameter [18].
Two experiments were conducted in order to characterize the proposed method. With the first one we characterized the precision of the proposed method. A 3D measurement of the same workpiece was repeated 10 times: the first measurement served as a target surface to which all the other surfaces (sources) were registered. The workpiece was not detached during the measurement repetitions.
With the second experiment the accuracy of the proposed method was assessed. In this case the systematic offsets along the X, Y and Z axes were incrementally added to the robot's scanning path. Offsets were executed in a 0.5 mm step in each direction separately, from 0 mm to 4 mm. The same workpiece was measured after each offset and the workpiece was not detached during the test. The measured surfaces were registered in order to determine the offsets, which should be equal to the offsets introduced by the robot.
For a statistical analysis of the differences between the results of the registration of the full and reduced point clouds an analysis of variance (ANOVA) at the 95 % level of significance was used. Except where noted differently, the results are provided in the notation Average value (Standard deviation). Geomagic Studio (Raindrop Geomagic) was used for the visualization of the results.
4 RESULTS
The laser-triangulation sensor measured the profiles for every 0.2 mm of the move, so that each measurement was composed of about 500 profiles, containing 659 points each. After the exclusion of the background points, the average number of points in each point cloud was 77582 (493). A typical 3D measurement of the workpiece is shown in Fig. 5a.
The MPD parameter of the ICP algorithm was set to 10 mm based on the maximal expected displacement. Other three parameters were set on default values: MNI = 100, NNE = 5 and N = 7, as it is advised in [18]. The influence of the geometry
and parameters required for the correct alignment is further described in [19]. Whether the registration was successful the deviation between the reference and the current surface was checked.
a)	b)
Fig. 5. a) Measured surface of the workpiece; and b) map of the deviations between the registered target and the source surface
The average deviations were 0.06 mm (0.08 mm) and a typical result is visualized in Fig. 5b. The average calculation time for the registration was 1.52 s (0.10 s). In order to speed up the calculation, we investigated the influence of the point reduction on the point cloud. The results show that if the number of points in the cloud is reduced to 6 % of the original value, the calculation time was reduced to 0.22 s (0.03 s), while the change in accuracy of the registration was not statistically significant (the ANOVA p-values were 0.25, 0.92 and 0.64 for translations in X, Y and Z directions, respectively).
The results of repeated measurements of the same workpiece without any intentionally introduced offset indicates that the calculated offset (bias) was 0.03 mm (0.06 mm), 0.01 mm (0.01 mm) and -0.07 mm (0.12 mm) in the X, Y, and Z directions, respectively.
The results of the second experiment are shown in Fig. 6. The diagram shows the Euclid distance between the measured and the robot offsets versus the offset from the origin. The mean Euclid distance of all the points is 0.23 mm (0.12 mm). The average distances in the individual directions are 0.05 mm (0.10 mm), -0.02 mm (0.03 mm), -0.03 mm (0.15 mm) in the X, Y, and Z directions, respectively. Although the workpiece was not deliberately rotated during the experiment, the registration algorithm returned minimal rotations of -0.02° (0.02°), 0.01°
(0.02°) and -0.02° axes, respectively.
(0.06°) around the X, Y and Z
0.5
'0.4
" 0.3 §
-d
X) 0.2
3 ш
0.1
i			□			□
mean + std	•			□		° +
						+
mean n	о	о		+ •	□ +	•
•	□		•		•	................"8"—
9 1 mean - std		□ •	+		0	Translation in: + X axis
	+	+	о	о		о Yaxis □ 7 axis • all 3 axes
						i Ь
Offset from origin [mm]
Fig. 6. Euclid distance between the measured and reference points to the displacement in each direction
5 DISCUSSION
The results of the registration were compared with the results acquired with commercially available software for 3D surfaces using the Global registration function in Geomagic Studio. Although the results using the ICP algorithm from our program were slightly better (approximately 0.01 mm), the differences were not statistically significant (ANOVA, p = 0.05).
We estimate that the accuracy of determining the positional misalignment of the workpiece is mainly affected by the accuracy of the robot arm. For our case it is approximately 0.1 mm, which means that about 50 % of the distances between the true and measured origins from the robot. On the other hand, since the same robot with the same accuracy is used for the deburring, the same accuracy can be attributed to the deburring system as a whole.
The second parameter that affects the precision of the whole system is the laser-triangulation sensor. Reference [20] shows that increasing the sensor precision from 1.6 mm to 0.3 mm improves the registration precision from 1.60° to 0.12°. In our case the sensor precision is 0.1 mm, which is very close to the standard deviations between the measured and true offsets.
Since the expected offsets between the true and the current workpiece position are less than 4 mm, we implemented the ICP algorithm, which found the local minimum in terms of deviation. In cases where larger displacements are expected, the introduction of course registration using feature detection and matching (e.g., SIFT [21], MSER [22], SURF [23] feature detectors) should be considered.
6	CONCLUSIONS
An adaptive robotic system was developed for the deburring of die-cast parts with position and orientation error correction. It uses a custom-developed laser-triangulation sensor to measure the 3D shape of the workpiece's surface and the ICP registration algorithm to determine the orientation and position error with respect to the reference workpiece, on which the robot was taught. The experiments show that the errors in the position and orientation in the range up to 4 mm are reduced to 0.23 mm (0.12 mm). The remaining errors can be further reduced by selecting a more accurate robot and triangulation sensor.
The developed adaptive system is also applicable in other similar applications where it is difficult to ensure the repeatable clamping of a workpiece.
7	REFERENCES
[1]	Song, H.C., Song, J.B. (2013). Precision robotic deburring based on force control for arbitrary shaped workpiece using CAD model matching. International Journal of Precision Engineering and Manufacturing, vol. 14, no. 1, p. 85-91, DOI:10.1007/s12541-013-0013-2.
[2]	Ton, T., Park, H., Ko, S. (2011). Experimental analysis of deburring process on inclined exit surface by new deburring tool. CIRP Annals Manufacturing Technology, vol. 60, no. 1, p. 129-132, D0I:10.1016/j.cirp.2011.03.124.
[3]	Valente, C.M.O., Oliveira, J.F.G. (2004). A new approach for tool path control in robotic deburring operations. ABCM Symposium Series in Mechatronics, vol. 1, p. 124-133.
[4]	Bagde, S.T. (2014). Development of combined deburring and inspection system. IOSR Journal of Mechanical and Civil Engineering, vol. 9, p. 63-68.
[5]	Kazerooni, H. (1988). Automated robotic deburring using impedance control. IEEE Control System Magazine, vol. 8, no. 1, p. 21-25, DOI:10.1109/37.464.
[6]	Asakawa, N., Toda, K., Takeuchi, Y. (2002). Automation of chamfering by an industrial robot; for the case of hole on free-curved surface. Robotics and Computer Integrated Manufacturing, vol. 18, no. 5-6, p. 379-385, DOI:10.1016/ s0736-5845(02)00006-6.
[7]	Srinivasan, H., Harrysson, O.L.A., Wysk, R.A. (2015). Automatic part localization in a CNC machine coordinate system by means of 3D scans. The International Journal of Advanced Manufacturing Technology, vol. 81, no. 5, p. 1127-1138, DOI:10.1007/s00170-015-7178-z.
[8]	Jaweera, N., Webb, P. (2010). Measurement assisted robotic edge deburring of aero engine components. WSEAS Transactions on Systems and Control, vol. 5, no. 3, p. 174-183.
[9]	Nagata, F., Kusumoto, Y., Fujimoto, Y., Watanabe, K. (2007). Robotic sanding system for new designed furniture with freeformed surface. Robotics and ComputerIntegrated
Manufacturing, vol. 23, no. 4, p. 371-379, D0l:10.1016/j. rcim.2006.04.004.
[10]	Song, H.C., Kim, B.S., Song, J.B. (2012). Tool path generation based on matching between teaching point and CAD model for robotic deburring. IEEE/ASME International Conference on Advanced Intelligent Mechatronics, p. 890-895, D0I:10.1109/ AIM.2012.6265921.
[11]	Habibi, B., Pescaru, S. (2004). Method and apparatus for single camera 3D vision guided robotics. US Patent 6,816,755, US Patent & Trademark Office, New York.
[12]	Biegelbauer, G., Vincze, M. (2006). 3D vision-guided bore inspection system. ICVS IEEE International Conference on Computer Vision Systems, p. 1-22, DOI:10.1109/ICVS.2006.1.
[13]	Skotheim, 0., Lind, M., Ystgaard, P., Fjerdingen, S.A. (2012). A flexible 3D object localization system for industrial part handling. IEEE/RSJ International Conference on Intelligent Robots and Systems, p. 3326-3333, D0I:10.1109/ IR0S.2012.6385508.
[14]	Rajaraman, M., DawsonHaggerty, M., Shimada, K., Bourne, D. (2013). Automated workpiece localization for robotic welding. IEEE International Conference on Automation Science and Engineering, p. 681-686, D0I:10.1109/CoASE.2013.6654062.
[15]	Yaskawa Motoman Robotics. MA1800 Arc Welding, from www. motoman.com/datasheet/MA1800.pdf accessed on 201510-22.
[16]	Jezeršek, M., Možina, J. (2003). A laser anamorph profilometer. Strojniški vestnik - Journal of Mechanical Engineering, vol. 49, no. 2, p. 76-89.
[17]	Besl, P.J., McKay, N.D. (1992). Method for registration of 3D shapes. Proceedings of SPIE 1611, Sensor Fusion IV: Control Paradigms and Data Structures, D0I:10.1117/12.57955.
[18]	Rusinkiewicz, S. (2015). trimesh, from gfx.cs.princeton.edu/ proj/trimesh2/ accessed on 2015-10-22.
[19]	Rusinkiewicz, S., Levoy, M. (2001). Efficient variants of the ICP algorithm. Proceedings of 3rd International Conference on 3-D Digital Imaging and Modeling, p. 145-152, D0I:10.1109/ IM.2001.924423.
[20]	Pavlovčič, U, Diaci, J., Možina, J, Jezeršek, M. (2013). Characterization of the headtotrunk orientation with handheld optical 3D apparatus based on the fringe projection technique. BioMedical Engineering OnLine, vol. 12, art. num. 96, D0I:10.1186/1475-925X-12-96.
[21]	Lowe, D.G. (1999). Object recognition from local scale-invariant features. The Proceedings of the 7th International Conference on Computer Vision, vol. 2, p. 1150-1157, D0I:10.1109/ICCV.1999.790410.
[22]	Matas, J., Chum, O., Urban, M., Pajdla, T. (2002). Robust wide-baseline stereo from maximally stable extremal regions. Image and Vision Computing, vol. 22, no. 10, p. 761-767, D0I:10.1016/j.imavis.2004.02.006.
[23]	Bay, H., Ess, A., Tuytelaars, T., Van Gool, L. (2008). Speeded-Up Robust Features (SURF). Computer Vision and Image Understanding, vol. 110, no. 3, p. 346-359, D0I:10.1016/j. cviu.2007.09.014.
Strojniški vestnik - Journal of Mechanical Engineering 62(2016)4, 213-219 © 2016 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2015.2528	Original Scientific Paper
Ultra-Precision Machining of a Large Amplitude Sinusoidal Ring Surface Based on a Slow Tool Servo
Shijun Ji - Huijuan Yu - Ji Zhao* - Xiaolong Liu - Mingxu Zhao Jilin University, School of Mechanical Science and Engineering, China
A sinusoidal ring surface with a large size, amplitudes in the macro-level, sub-micrometre form accuracy and nanometre surface roughness plays an important role in modern optical systems and precision calibration, but it is difficult to fabricate with traditional cutting methods. A sinusoidal ring surface with a submillimetre amplitude, machined using a slow tool servo (STS) assisted by single-point diamond turning (SPDT) is presented in this paper. The kinematic characteristic of machine axes is analysed to assess the feasibility of the turning process. The tool path generation, tool geometry optimization, and tool radius compensation are investigated for fabrication of the desired surface. A fabricating experiment of a 0.4 mm amplitude sinusoidal ring surface is performed, and the surface profile is measured with a Talysurf PGI 1240, a measurement system for the small to medium sized aspheric optics. After dealing with original measured curves, a form accuracy 0.274 pm in peak-to-valley error (PV) and surface roughness 7.5 nm in Ra are obtained for the machined surface, and it can be seen that STS can be used for machining of a large amplitude sinusoidal ring surface. Keywords: sinusoidal ring surface, STS, kinematic characteristic, ultra-precision machining
Highlights
•	A sinusoidal ring surface with a submillimetre amplitude is machined using a slow tool servo (STS) assisted by single-point diamond turning (SPDT).
•	The kinematic characteristic of machine axes is analysed to assess the feasibility of the turning process.
•	Analysis about tool path generation and tool radius compensation, and optimization of the cutting tool parameters are necessary before the actual cutting experiment.
•	A form accuracy 0.274 pm in peak-to-valley error (PV) and surface roughness 7.5 nm in Ra are obtained for the machined surface.
0 INTRODUCTION
Freeform surfaces can make a smaller, lighter, higherperformance optical system and have been the subject of increased attention in recent years [1] to [3]. The sinusoidal ring surface is a typical shape of freeform surfaces, which can be used in positioning multi-axis precision motion stages and as aerospace components. The surfaces are required to be of large size, and have sub-micrometre form accuracy and a nanometre surface finish in those applications [4] and [5].
Although the sinusoidal ring surface is a rotationally symmetric surface, the continuously changing position and moving in the Z direction still make the processing difficult to perform while the workpiece is turned with a spindle. Several methods can be used to manufacture the sinusoidal ring surface, such as grinding, polishing and lapping. However, these methods use multiple processes, and it is not easy for them to fabricate complex surfaces, such as sinusoidal surfaces with high form accuracy and highly smooth surfaces [6]. Presently, SPDT with fast tool servo (FTS) and STS are becoming pre-eminent due to their ability to achieve optical quality in a one-stage process without subsequent processing [7] and [8]. As for FTS, the cutting tool can be rapidly and accurately
positioned in the manufacturing process; therefore, a freeform surface with micro/nano-structures can be fabricated. However, the amplitude of the swing is usually small, which cannot be used in large sagittal surface machining. For example, Lu et al. [9] made a sinusoidal grid surface using a fast tool servo system with a profile wavelength is 350 дт, and a peak-to-valley amplitude of approximately 16 ^m. Rasmussen et al. [10] designed a piezoelectric-driven cutting tool system, which can dynamically control the depth of cut; it can be used to machine slightly non-circular workpieces. With an amplification mechanism, the tool was able to produce a 52 ^m travel range with a frequency of 200 Hz. As for STS, all axes can be under full coordinate position control. Moreover, the form accuracy and the quality of a machined surface depend on the proper selection of the tool path, cutting strategies and cutting parameters [11] to [13]. Several papers have examined sinusoidal surfaces in recent years. For instance, Zhang et al. [1] presented studies on the sinusoidal surface machining with a proposed cylindrical coordinate method using the STS. With the proposed method, a 40 ^m amplitude sinusoidal surface with a nanometric finish of 5.54 nm in Ra is achieved. Gao et al. [6] manufactured a sinusoidal angle grid surface with a wavelength of 100 ^m and
an amplitude of 100 nm over a large surface area with a diameter of 150 mm. Zhu et al. [14] described the manufacturing process of a sinusoidal grid surface using FTS on a Zr-based bulk metallic glass (BMG); an 80 дт amplitude sinusoidal surface is obtained using diamond cutting tools.
The challenge of manufacturing a submillimetre amplitude sinusoidal ring surface using a three-axis (X, Z, C) turning machine is to achieve the ultra-precision form accuracy [15] and [16]. Thus, this paper presents a theoretical research on the machining of sinusoidal ring surface using the STS, including tool path generation, kinematic characteristic of two linear axes and the selection of cutting tool geometry. Then, a typical sinusoidal ring surface is machined with STS to demonstrate the feasibility of the theoretical analyses.
1 TOOL PATH GENERATION
In single-point diamond turning, a workpiece is mounted on the spindle, and the schematic diagram of cutting the sinusoidal ring surface using STS is shown in Fig. 1, where a and h are the wavelength and the double amplitude of the sinusoidal surface, respectively. The sinusoidal ring surface in cylindrical coordinate system (X, в, Z) can be expressed as:
h . 2n
z = — sm I — x 2 I a

2
(1)
where D is the maximum diameter of the sinusoidal ring surface.
Diamond tool
Fig. 1. Schematic diagram of cutting a sinusoidal ring surface
The diamond tool is moved along the Z direction in coordination with the movement of the workpiece
along the X direction in the single-point diamond turning process. To obtain the continuous and smooth surface, the tool path of the diamond tool must be accurately generated in advance, and a planar spiral is used to generate the cutting tool path. The path of tool contact points is obtained subsequently by projecting the planar spiral on the desired surface in the Z direction. The planar spiral curve is related to the cutting parameters [14], such as feed rate f [mm/ min] and spindle speed v [rev/min]. The expression in polar coordinates can be shown as follows:
f 2nv
0,
x * D 2
(2)
According to Eqs. (1) and (2), the calculated path of tool contact points for the sinusoidal ring surface can be illustrated by Fig. 2.
Fig. 2. Path of tool contact points for sinusoidal ring surface
To avoid the influence of tool radius and to obtain a surface with superior form accuracy, the geometry of the diamond tool nose must be considered, since it is not an ideal point but a half-cylinder with radius R [17] and [18]. In fact, by projecting points of the planar spiral on the desired surface, cutting contact points can then be obtained. In order to generate the motion commanders of CNC machine tools, the cutting contact points must be converted into cutter location points; that is to say, the tool radius compensation must be conducted. Since every cutting contact point has a unique curvature along the cutting trajectory on the machining surface, the tool compensation can be achieved by offsetting the tool contact points with a value of tool radius R in the normal vector direction rather than only by offsetting the tool contact points with a value of tool radius R in the Z direction. The
a)
x' = x + Rsint y' = y + Rcost
b)	x[mm]
Fig. 3. Tool radius compensation a) the schematic diagram of tool radius compensation b) tool radius compensation for sinusoidal ring surface
schematic diagram of tool radius compensation is illustrated in Fig. 3, where point P (x,y) is an arbitrary point on the tool cutting contact path. The normal vector n at point P on the tool cutting contact path determined the compensation direction, and the offset value is the tool nose radius R; then, P'(x',y') is the corresponding compensation tool location point. All of the compensation tool locational points constitute the tool location path. The tool location point can be calculated by:
where в is the angle between normal vector n and perpendicular line. After tool radius compensation with Eq. (3), Fig. 4 shows the final tool location path of sinusoidal ring surface, which drives the diamond tool accurately in a ultra-precision machine.
2 ANALYSIS OF KINEMATIC CHARACTERISTICS
Since the machining efficiency is greatly influenced by the velocity and accelerated speed limitation of every axis, it is necessary to analyse the kinematic characteristic of every axis before cutting. The movement of cutting point in X axis and Z axis along with the turning of C axis is illustrated in Fig. 5 under the Х-в-Z coordinate system.
To analyse the kinematic characteristic of the fabricating process, the cutting surface has been stretched out as shown by Fig. 5. Obviously, the expanded shape of the sinusoidal ring surface is a straight sinusoidal wave surface. The point Pt is an arbitrary cutting point on the cutter path, and the displacement of X axis to the cutting point is a linear relation with the angular displacement в in Х-в-Z coordinate system. While, the relation between displacement of Z axis and angular в can be indicated by:
= z (e) =—sin \f( W 2 I av
+ z0.
(4)
Fig. 4. Tool path for sinusoidal ring surface in 3D
Fig. 5. The kinematic analysis of diamond tool
As indicated by Eq. (4), the displacement of Z axis of the cutting point is a sinusoidal function of angular displacement в, and the amplitude of the sinusoidal function is constant h/2; therefore, the wavelength is changed to 2nav/f. Moreover, it also can be seen that the position and moving direction in the Z direction are continuously changeable with the
change of angular displacement в, so the STS of Z axis is inevitably used to provide accurate displacement control. The speed and acceleration of the X and Z axes can be obtained by calculating the first derivative and second derivative of displacement function X=Х(в) and Z=Z(e), respectively. The maximum velocity and maximum acceleration of linear axes should not exceed the corresponding maximum value of the machine tool, which are 4000 mm/min and 1980 mm/ s2 for the Nanoform250 machine tool.
3 PARAMETER SELECTION OF DIAMOND TOOL
In the ultra-precision turning process, the proper selection of the cutting tool plays an important role in the processing planning, since the quality of the final product is greatly influenced by the diamond tool geometry [19] to [21]. Thus, cutting tool geometry must be optimized according to practical cutting conditions. [22].
The cutting tool radius selection of the diamond tool is presented in Fig. 6, which is based on the stretch-out view of cutting surface. As shown in Fig. 6a, point Pi is an arbitrary cutting point on the tool path, and vector nti is the tangent vector for path curve at point Pi. Point P'i and n'ti are the projection point and vector of Pi and nti on the X-O-в plan, respectively. Plane у is given by passing P'i and vertical to n'ti, and line lr is the intersection line of plane у and the cutting surface. The enlarged partial view of the cutting condition is presented in Fig. 6b, which illustrates that the cutting interference can occur when curvature radius and arc angle of cutting point are both smaller than the corresponding value R and в of the diamond tool, such as cutting point Pi, and the opposite condition is applicable, such as point Pk. Therefore, the optimal cutting tool radius and included angle should be less than the minimum value of curvature radius and arc angle at all cutting points.
Clearance angle is another indispensable parameter for the selection of tool geometry. As shown in Fig. 7a, plane ö is established by passing the point P'i and vector n'ti, and line lc is the intersection line of plane ö and the cutting surface. The selection of clearance angle at one cutting point is determined by the inclination angle of the tangential line at this point on the cutting curve, and the cutting interference occurs in the case that the clearance angle less than the inclination angle of the tangential line, such as point Qi, while the opposite condition point Qj is feasible, as demonstrated in Fig. 7b. With the purpose of avoiding interference at any cutting points, the optimal clearance angle should be selected for machining, and
R
a)
b) 11	Pk
Fig. 6. Optimal tool radius selection; a) schematic diagram of radius selection, and b) partial enlarged detail view
a)
b)
Fig. 7. Optimal clearance angle selection; a) schematic diagram of clearance angle selection, and b) partial enlarged detail view
the complementary angle of the optimal clearance angle must be larger than the maximum value of the inclination angle in all cutting point positions.
4 EXPERIMENTS AND DISCUSSION
A sinusoidal ring surface expressed by Eq. (1) with h = 0.4 mm, a = 2 mm and D = 20 mm was simulated and fabricated in the Nanoform250 ultra-precision machining system. Both the X and Z axes of the Nanoform250 are constructed with non-contact oil hydrostatic bearing sideways, non-contact linear motors, and non-contact encoder feedback, liquid cooled systems for slides. In this manner, there is zero stick-slip or friction and thermal stability. The material of the workpiece is Al-6061. The diamond tool with a nose radius of 0.506 mm, an included angle of 120°, and a clearance angle of 10° was applied, which satisfies the analyses in Section 4. The machining quality on computerized numerically controlled (CNC) machine tools is sensitive to the machining parameters [23]. In this paper, the cutting conditions we adopted are listed in Table 1.
Table 1. Cutting conditions
Parameters	Values
X increment per rev	0.01 mm
C increment	1°
Max speed of X axis	0.011 mm/s
Max speed of Z axis	0.0063 mm/s
Frequency of Z axis	0.0055 Hz
Acceleration of X and Z axes	0.2 g
Amplitude of the surface	0.4 mm
were conducted to ensure high machining precision. The aluminium workpiece was fabricated with STS at a spindle speed of 60 rpm. The actual fabrication structure is displayed in Fig. 8. A rotationally symmetrical fixture was applied to connect the vacuum chuck and workpiece. The designed cutting path, as shown in Fig. 4, was employed to perform the cutting experiment.
Before ultra-precision manufacturing the workpiece, centring the tool with an LVDT video tool set station and dynamic balancing of the spindle
Fig. 8. Photos of in-situ machining
To reveal the surface quality of the machined part, the surface roughness and form accuracy of the machined sinusoidal ring surface were measured and evaluated with a high magnification microscope (Contour GT) and a precision form measurement system (Taylor Hobson), respectively. Firstly, Contour GT, a white-light interferometric microscope, is used to measure surface roughness. The measuring area is a rectangular shape of 60 цт x 45 цт, which is about 3 mm to the outside boundary of the workpiece. In the process of measurement, the
a)
b)
Fig. 9. Measurement result of surface roughness for machined sinusoidal surface: a) linear roughness measur. b) areal roughness measur.
Measure length [mm]
Fig. 10. The rotated measured profile curve and theoretical sinusoidal curve
machined surface was placed under the optical lens of the Contour GT. Then, the adjustment button was rotated forward and backward repeatedly until the interference fringe appeared. After that, the instrument was set to automatic measurement mode. Through automatic measuring the surface of the part and analysis with its corresponding software, the testing result is given in Fig. 9. In Fig. 9a, cutting direction and measuring direction are illustrated. Linear roughness measurement is used to assess the surface quality, the value of surface roughness (Ra) along the measuring direction line is about 7.5 nm. In Fig. 9b, areal roughness measurement is used to evaluate the surface quality, the value of surface roughness (Sa) is 4.6 nm. It can be seen that the surface roughness of the machined part meets the requirement of the ultra-precision machining.
The measurement profile of the sinusoidal ring surface was performed by Taylor Hobson. To obtain
X 10""
		
*V A. л-v JЛ г— J	A 1 Л К, и,	La_A_A./
Q	A*.	г /III" ^ г 1 ..... t гЯ\ ^J г	г	г \ . У r_^
0 2 4 6 8 10 12 14 16 18 20
Diameter [mm]
Fig. 11. The result of form accuracy for machined sinusoidal ring surface
form accuracy of the machined surface, further analysis for the profile is necessary. Fig. 10 represents the comparison between the measured sinusoidal profile which has been obtained by rotating the original measurement sinusoidal curve and theoretical sinusoidal curve. After post-processing, it can be seen they have a good match. The deviation along the normal direction between the rotated sinusoidal profile and the theoretical one is considered to be the form accuracy, and the result is demonstrated in Fig. 11. It can be observed that the form accuracy PV value 0.274 цт is obtained for the machined sinusoidal ring surface with STS.
5 CONCLUSIONS
In this paper, the STS assisted diamond turning is applied to machine a sub-millimetre sinusoidal ring surface in ultra-precision. For the purpose of improving the machining precision, several theoretical analyses are implemented before the actual cutting experiment, including tool path generation, tool radius compensation, kinematic analyses of axes, and optimization of the cutting tool parameters. To evaluate the form accuracy of the machined surface, the measured sinusoidal profile is rotated to match the theoretical sinusoidal profile, and the deviation value of two profile curves can be obtained. To verify the validity of theoretical analyses, a sinusoidal ring surface is manufactured with the developed theories. The results of surface roughness and form accuracy of the machined surface are 7.5 nm in Ra and 0.274 цт in PV, respectively. The experiment results reveal that the theoretical analyses are effective to improve the machining precision.
6 ACKNOWLEDGEMENTS
This work is supported by the National Key Basic
Research and Development Program (973 program)
of China (grant no. 2011CB706702), Natural Science
Foundation of China (grant nos. 51305161 and
51135006), Jilin Province Science and Technology
Development Plan, item (20130101042JC).
7 REFERENCES
[1]	Zhang, X.D., Fang, F.Z., Wang, H.B., Wei, G.S., Hu, X.T. (2009). Ultra-precision machining of sinusoidal surfaces using the cylindrical coordinate method. Journal of Micromechanics and Micro-engineering, vol. 19, no. 5, p. 054004, DOI:10.1088/0960-1317/19/5/054004.
[2]	Wang, Y., Lu, Z., Liu, H., Zhang, H. (2007). Application of freeform surface prism. Infrared and Laser Engineering, vol. 36, no. 3, p. 319-321, D0I:10.3969/j.issn.1007-2276.2007.03.009. (in Chinese)
[3]	Liu Y. (1998). Designing the integrated helmet-mounted display for pilots. Proceedings of SPIE 3560, Display Devices and Systems, p. 175-181, DOI:10.1117/12.319671.
[4]	Gao, W., Dejima, S., Yanai, H., Katakura, K., Kiyono, S., Tomita, Y. (2004). A surface motor-driven planar stage integrated with a XY0Z surface encoder for precision positioning. Precision Engineering, vol. 28, no. 3, p. 329-337, DOI:10.1016/j. precisioneng.2003.12.003.
[5]	Jiang, X., Scott, P., Whitehouse, D. (2007). Freeform surface characterization - a fresh strategy. CIRP Annals - Manufacturing Technology, no. 56, no. 1, p. 553-556, DOI:10.1016/j.cirp.2007.05.132.
[6]	Gao, W., Araki, T., Kiyono, S., Okazaki, Y., Yamanaka, M. (2003). Precision nano-fabrication and evaluation of a large area sinusoidal grid surface for a surface encoder. Precision Engineering, vol. 27, no. 3, p. 289-298, DOI:10.1016/S0141-6359(03)00028-X.
[7]	Zhang, X.D., Fang, F.Z., Wu, Q.Q., Liu, X.L., Gao, H.M. (2013). Coordinate transformation machining of off-axis aspheric mirrors. The International Journal of Advanced Manufacturing Technology vol. 67, no. 9-12, p. 2217-2224, DOI:10.1007/ s00170-012-4642-x.
[8]	Gong, H., Fang, F.Z., Hu, X.T. (2012). Accurate spiral tool path generation of ultra-precision three-axis turning for non-zero rake angle using symbolic computation. The International Journal of Advanced Manufacturing Technology, vol. 58, no. 9-12, p. 841-847, DOI:10.1007/s00170-011-3433-0.
[9]	Lu, H., Lee, D.-W., Lee, S.-M., Park, J.-W. (2012). Diamond machining of sinusoidal grid surface using fast tool servo system for fabrication of hydrophobic surface. Transactions of Nonferrous Metals Society of China, no. 22, suppl. 3, p. 787792, DOI:10.1016/S1003-6326(12)61805-6.
[10]	Rasmussen, J.D., Tsao, T.C., Hanson, R.D., Kapoor, S.G. (1994). Dynamic variable depth of cut machining using piezoelectric actuators. Journal of Machine Tool and Manufacture, vol. 34, no. 3, p. 379-392, DOI:10.1016/0890-6955(94)90007-8.
[11]	Yin, Z.Q., Dai, Y.F., Li, S.Y., Guan, C.L., Tie, G.P. (2011). Fabrication of off-axis aspheric surfaces using a slow tool servo. International Journal of Machine Tools and Manufacture, vol. 51, no. 5, p. 404-410, DOI:10.1016/j. ijmachtools.2011.01.008.
[12]	Yi, A.Y., Li, L. (2005). Design and fabrication of a micro-lens array by use of a slow tool servo. Optics Letters, vol. 30, no. 13, p. 1707-1709, DOI:10.1364/OL.30.001707.
[13]	Li, L., Allen, Y.Y., Huang, C., Guan, C.L., Tie, G.P. (2006). Fabrication of diffractive optics by use of slow tool servo diamond turning process. Optical Engineering, vol. 45, no. 11, p. 113401, DOI:10.1117/1.2387142.
[14]	Zhu, Z., Zhou, X., Liu, Q., Lin, J.; Zhao, S. (2012). Fabrication of micro-structured surfaces on bulk metallic glasses based on fast tool servo assisted diamond turning. Science of Advanced Materials, vol. 4, no. 9, p. 906-911, DOI:10.1166/ sam.2012.1374.
[15]	Wang, Y., Zhao, Q., Shang, Y., Lv, P., Guo, B., Zhao, L. (2011). Ultra-precision machining of Fresnel microstructure on die steel using single crystal diamond tool. Journal of Materials Processing Technology, vol. 211, no. 12, p. 2152-2159, DOI:10.1016/j.jmatprotec.2011.07.018.
[16]	Ding, X., Rahman, M. (2012). A study of the performance of cutting polycrystalline Al 6061 T6 with single crystalline diamond micro-tools. Precision Engineering, vol. 36, no. 4, p. 593-603, DOI:10.1016/j.precisioneng.2012.04.009.
[17]	Fang, F.Z., Venkatesh, V.C. (1998). Diamond cutting of silicon with nanometric finish. CIRP Annals - Manufacturing Technology, vol. 47, no. 1, p. 45-49, DOI:10.1016/S0007-8506(07)62782-6.
[18]	Rakuff, S., Cuttino, J.F. (2009). Design and testing of a longrange, precision fast tool servo system for diamond turning. Precision Engineering, vol. 33, no. 1, p. 18-25, DOI:10.1016/j. precisioneng.2008.03.001.
[19]	Yu, X., Wang, L. (1999). Effect of various parameters on the surface roughness of an aluminum alloy burnished with a spherical surfaced polycrystalline diamond tool. International Journal of Machine Tools and Manufacture, vol. 39, no. 3, p. 459-469, DOI:10.1016/S0890-6955(98)00033-9.
[20]	Liu, K., Melkote, S.N. (2007). Finite element analysis of the influence of tool edge radius on size effect in orthogonal micro-cutting process. International Journal of Mechanical Sciences, vol. 49, no. 5, p. 650-660, DOI:10.1016/j. ijmecsci.2006.09.012.
[21]	Zong. W.J., Li, D., Cheng, K., Sun, K., Liang, Y.C. (2007). Finite element optimization of diamond tool geometry and cutting-process parameters based on surface residual stresses. The International Journal of Advanced Manufacturing Technology, vol. 32, no. 7-8, p. 666-674, DOI:10.1007/s00170-005-0388-z.
[22]	Wang, L.,Wang, D., Gao, Y. (2015). Investigations on the effects of different tool edge geometries in the finite element simulation of machining. Strojniški vestnik - Journal of Mechanical Engineering, vol. 61, no. 3, p. 157-166, DOI:10.5545/sv-jme.2014.2051.
[23]	Čuš, F., Župerl, U. (2015). Surface roughness control simulation of turning processes. Strojniški vestnik - Journal of Mechanical Engineering, vol. 61, no. 4, p. 245-253, DOI:10.5545/sv-jme.2014.2345.
Strojniški vestnik - Journal of Mechanical Engineering 62(2016)4, 220-230 © 2016 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2015.3266	Original Scientific Paper
Estimating the Strain-Rate-Dependent Parameters of the Cowper-Symonds and Johnson-Cook Material Models
using Taguchi Arrays
Andrej Škrlec* - Jernej Klemene
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
In order to reduce R&D costs, a product's behaviour during use is predicted with numerical simulations in the early phases of R&D. If a structure is subjected to high-strain-rate loading this effect should be considered in the material models that are used for the numerical simulations. This article shows how the strain-rate-dependent material parameters can be determined by combining the experimental data with the numerical simulations. The presented methodology is based on the application of Taguchi arrays to find the most appropriate values for the strain-rate-dependent parameters. The presented methodology is applied in a practical case, for which the parameters of the Cowper-Symonds and Johnson-Cook material models are estimated.
Keywords: strain-rate-dependent material behaviour, finite-element method, explicit dynamics, Taguchi arrays Highlights
•	A procedure for modelling the strain-rate-dependent behaviour of materials is presented.
•	The parameters of the Cowper-Symonds and Johnson-Cook material models are estimated on the basis of an impact experiment.
•	A simulation plan for estimating the material parameters is developed with the help of Taguchi orthogonal arrays.
0 INTRODUCTION
To reduce the costs of a research and development (R&D) process and optimise the product's design, while ensuring the necessary safety, effectiveness and reliability of the newly developed product, finite-element (FE) simulations of the product's behaviour under real operating conditions are widely applied. To obtain reasonable simulation results the operating conditions as well as the product's geometry and its material properties must be known. Usually, in the R&D process the structural loads are assumed on the basis of similar products, are obtained by numerical simulations or are defined by the customer. If we assume that the structural loads resulting from the operating conditions are known, it is the material properties that have the greatest influence on the product's behaviour for a given geometry. If the structure is subjected to extreme mechanical loading conditions (e.g., impacts during crash tests or different burst tests) it is of tremendous importance to consider the strain-rate-dependent material properties when performing a FE simulation.
It is known from the literature [1] that quasi-static loading does not have a significant influence on the material's yield stress and the stress-strain relationship, but this changes if the strain rate increases. The increased values of the strain rate cause an increase in the material's yield stress and change the material's stress-strain behaviour in the
plastic domain. Investigations of this effect have been performed by many researchers over the last century, like Hopkinson, Charpy, Taylor [1], Zerilli, Armstrong, Johnson [2], [3], etc. In the last 20 years, strain rate influence on material behaviour is still interesting for researchers like El-Magd [4], Zhao and Gary [5], Huh et al. [6] to [8], etc.
When using the explicit dynamic FE code to simulate extreme loading conditions, such as impact phenomena, the material models that consider the strain-rate dependency of the material's plastic curve are commonly applied. The three most commonly applied material models in researches [9] to [15] that consider the strain-rate effects are: Cowper-Symonds, Johnson-Cook, and Zerilli-Armstrong. Since the Cowper-Symonds and Johnson-Cook material models are simpler than the Zerilli-Armstrong material model, we considered only the former two models in our research. The main difference between them is how they account for the strain-rate effects. Consequently, the number of material parameters that describe the plastic stress-strain relationship with the strain-rate effects is different ([16] and [17]): • Yield stress according to the Cowper-Symonds material model:
^y =
a0 +ß-
E, ■ E E - E
-s
eff

(1)
Flow stress according material model [3]:
to the Johnson-Cook
G flow ~

B (f )
1 + c ■ In
и p
T - T
T - T
melt room
(2)
where a0 is the reference yield stress, E is the material's elastic modulus, Et is its tangent modulus, ß is the hardening coefficient, ePpf is the effective plastic strain, è is the strain rate and è0 is the reference strain rate. C and P are the strain-rate parameters of the Cowper-Symonds material model; B, n, c and m are the strain- and strain-rate-dependent parameters of the Johnson-Cook material model. The Johnson-Cook material model considers the influence of temperature on the stress-strain behaviour, but not the Cowper-Symonds material model. However, the temperature effects can be omitted from the Johnson-Cook material model if its parameter m is set to zero. The temperature influences can be omitted from the Johnson-Cook material models if they do not influence the material behaviour. The Cowper-Symonds and Johnson-Cook models are the most commonly applied material models when performing explicit dynamic simulations up to moderate strain rates, (i.e., up to è = 104) ([14] and [15]), at which the changing temperature does not have a significant influence.
The main problem, linked to the material parameters that consider the strain-rate effects, is that they cannot be simply measured and determined, thus they are empirically determined through special experimental and optimisation processes ([9] to [13]). For the above-mentioned material models the parameters that consider the strain-rate dependence were investigated for many different materials. In the literature, their typical values for mild steel ([9], [10] and [11]), high-yield-strength steel [12], aluminium alloys [13], titanium alloys [13], etc. can be found. Since, on the other hand, no link between the strain-rate parameters and the chemical composition of the material was found, these parameters should be identified individually for each material under consideration.
In this article we will compare the performance of Cowper-Symonds and Johnson-Cook material models when applied for simulating the impact of a ball on a thin steel plate. The parameters of the two material models, which cannot be estimated from the tensile test, were determined from the impact experiments using the LS-DYNA explicit FE code that was combined with a Taguchi array to reduce the numerical processing effort. The article is structured
as follows. After the introductory section and the theoretical background, the experimental arrangement and the experimental results are presented. The article continues with a presentation of the results and their discussion and ends with a concluding section and a list of references.
1 THEORETICAL BACKGROUND 1.1 Applied Material Models
In this article two material models (i.e., the Cowper-Symonds model and a simplified Johnson-Cook model) that consider the strain-rate effects on the material's behaviour are compared. For each material model three parameters were estimated on the basis of the impact experiment by using the explicit dynamic simulations combined with a parameter factorization according to the Taguchi array.
The Cowper-Symonds (C-S) material model with a bilinear characteristic is determined with the following parameters ([16] and [17]): elastic modulus E, Poisson's number v, tangent modulus Et, hardening coefficient ß, material density p, and the parameters C and P that describe the dependency of the yield stress ay on the strain rate è , see Eq. (1). Therefore, the flow stress is calculated as follows:
-eff

E, ■ E E - E,
eff
i
'eff■
(3)
For the C-S material model the strain rate influences only the yield stress ay. This means that the plastic curves (the flow stress as a function of strain) are parallel. The larger the strain rate, the higher the flow-stress curve, see Fig. 1a. For high-strain-rate applications the parameters C and P in Eqs. (1) and (3) are usually not estimated from the tensile test, due to the limitations of the existing tensile-test equipment (maximum strain rates are up to a few hundreds of s-1). Together with the tangent Et they were determined on the basis of the impact experiment.
Since the temperature effects were neglected in our case, the simplified Johnson-Cook (J-C) material model was applied with the parameter m from Eq. (2) being equal to 0. This material model is determined with the following parameters ([16] and [17]): elastic modulus E, Poisson's number v, reference yield stress s0, exponent of the flow-stress curve n, scale factor B for the effective plastic strain, sensitivity c to the logarithm of the strain rate and material density p. The flow stress is then given by:
m
v s0 J
Effective plastic strain - eeff [-]	Effective plastic strain - eeff [-]
Fig. 1. Flow stress for a) the Cowper-Symonds material model and b) the Johnson-Cook material model
	( p \n		1 + c - In	fs p 1
Gflow	CTo + B )			
	L J			V so )
(4)
For the J-C material model the strain rate influences the flow stress Oßow in its whole range and the flow-stress curves are not parallel. Their non-linearity depends on the exponent n. The larger the strain rate, the higher the flow-stress curve, see Fig. 1b. The three parameters B, n and c in Eq. (4) cannot be estimated from the tensile test and were determined from the impact experiment.
1.2 Simulation Plan for Material-Parameter Estimation
For each of the two strain-rate-dependent material models, three parameters (i.e., Et, C and P for the C-S material model and B, n and c for the J-C material model) need to be determined. These parameters are usually determined using a reversed engineering approach with the help of numerical simulations that reproduce the actual experiment ([11], [12] and [18]). This means that a series of numerical simulations with different combinations of material parameters are carried out to establish which combination of material parameters best fits the experimental results. The problem is that the strain-rate-dependent parameters can be selected from a domain with a range that spans over many orders of magnitude.
For this reason it is almost impossible to run a full-factorial simulation plan to determine the optimal strain-rate-dependent parameters, because the processing time would be prohibitive even in the case when simulations are run on a supercomputer. To shorten the processing time for estimating the material parameters it was decided to apply orthogonal Taguchi arrays for a simulation-plan set-up. The
reason for this decision was that the use of Taguchi arrays results in a significant reduction in the number of studied parameter combinations. Their advantage is that the assigned combinations of parameters are approximately equally distributed over the search space, which is not always the case with, e.g., Latin hyper-cubes.
We applied a L81(340) Taguchi array and transformed it into the L81(910) orthogonal array using the linear graph in Fig. 2 ([19] and [20]). This was done because we needed as many levels per parameter as possible and L81(910) best suited this requirement.
Fig. 2. Linear graph for the transformation of a L81(340) Taguchi array into a L81(910) array
Using the L81(910) orthogonal array the three material parameters (Et, C and P for the C-S material model and B, n and c for the J-C material model) are attributed to the first three columns of the L81(910) array. The levels were chosen from a very wide domain to account for all the possible parameter values for different kinds of steels. This makes the proposed approach general, although it was tested for the case of E185 steel. Since the ranges of the individual parameters span over more than one order
of magnitude, the logarithms of these parameters represent the nine-level factors in the Taguchi array. In this manner 81 combinations of the parameter triples (Et, C, P) and (B, n, c) were obtained that cover the whole domains of these parameters, see Fig. 3 for the parameters of the C-S material model and Fig. 4 for the parameters of the J-C material model.
Table 1. Parameters levels for the C-S material model
Original parameter levels of parameter C [ms-i]		
0.2154	0.4642 1.0000 2.1544	4.6416
10.0000	21.5443 46.4159 100.0000	
Original parameter levels of parameter P [-]		
1.0000	1.7783 3.1623 5.6234	10.0000
17.7828	31.6228 56.2341 100.0000	
Original parameter levels of parameter Et [GPa]		
0.1000	0.1778 0.3162 0.5623	1.0000
1.7783	3.1623 5.6234 1 0.0000	
Fig. 3. Distribution of C-S material-parameter combinations over their domains
Fig. 4. Distribution of J-C material-parameter combinations over their domains
By applying the Taguchi arrays for the simulation plan the number of material-parameter combinations that need to be simulated was reduced by a factor of nine when compared to the full-factorial simulation plan.
Table 2. Parameter levels for the J-C material model
Original parameter levels of parameter B [GPa]	
0.1000	0.1778 0.3162 0.5623 1.0000
1.7783	3.1623 5.6234 10.0000
Original parameter levels of parameter n [-]	
0.001	0.0024 0.0056 0.0133 0.0316
0.0750	0.1778 0.4217 1.0000
Original parameter levels of parameter c [-]	
0.001	0.0024 0.0056 0.0133 0.0316
0.0750	0.1778 0.4217 1.0000
For each combination of the material parameters six numerical simulations were carried out to account for the six different boundary and initial conditions, see Table 3. Therefore, when applying the simulation plan according to the L81(910) Taguchi array 816 = 486 numerical simulations were performed for each of the two material models.
A cost function that measures the deviations of the experimental and simulation results was defined so that it measured the difference between the experimentally determined and simulated data for the indentation depth H and the position of the indentation centre Z for the specimen. The averages of the experimentally determined values for these two geometrical parameters for two specimen thicknesses and three different velocities of the ball are listed in Table 3. The cost function that was used to assess the goodness-of-fit between the experiments and the simulations is defined as follows for the C-S and J-C material models:
( k w,
f ( E,, P, C ) = -k
t(( -H*š) +
i=1
+ (1 - Wh )-t ( - Zsš )
V
(5a)
f ( B, n, c) =-k
f k WH „
i=1
ZK - Hsiš ) + i=1
+(1 -wH)-Z( -Zsš)
(5b)
where к is the number of experimental results with different boundary and initial conditions (according to Table 3, к = 6). Hexp and Zexp, are the averaged
measured maximum indentation depth and its z coordinate for the specimen. Hsim and Zsim are the maximum indentation depth and its z coordinate for the specimen obtained by simulations. wH was a weighting factor in the two cost functions and was equal to 0.5 in our case.
2 EXPERIMENTAL ARRANGEMENT
2.1 Measurement of a Static Stress-Strain Curve
The methodology that was described in Section 1 was applied to characterise the strain-rate-dependent material behaviour of a mild steel E185. Its static material properties (elastic modulus, yield stress, ultimate tensile stress) were measured according to the ASTM E8/E8M standard [21] on Zwick/Roell Z050 testing equipment.
Fig. 5. Zwick/Roell Z050 test stand and the specimen
The test stand and the specimen geometry are presented in Fig. 5. A total of 21 specimens were tested. The average yield stress and the ultimate tensile strength were 185 MPa and 350 MPa, respectively.
450
400
350
£"300 a
£^250 b
200 150 100 50
0 _i_i_i_i_i_i_i_i_i_
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
The measured engineering stress-strain curves are presented in Fig. 6a and the resulting average true-stress-true-strain curve is presented in Fig. 6b. This true-stress-true-strain curve was determined with Eqs. (6) and (7), see Dowling [22].
ё = In (1 + ё),	(6)
j = j(l + s).	(7)
e and a were the corresponding average engineering strain and stress, respectively.
2.2 Experimental Determination of the Material Behaviour at High Strain Rates
The main objective was to determine the strain-rate-dependent material parameters for simulating the behaviour of a mild-steel sheet metal that is used as a shield during a turbine burst test. To identify the corresponding material parameters we designed experimental apparatus for shooting a steel ball at a flat specimen. The experimental arrangement was based on the ASTM D5420 standard [23], which describes a test method for measuring the impact resistance of a flat rigid plastic specimen by means of a striker impacted by a falling weight. During the impact between the ball and the flat specimen strains can be measured on the left-hand side of the specimen with strain gauges. After the impact test the gross geometric data, i.e., the indentation depth H and the position of the indentation centre Z, were measured, see Fig. 7.
The experimental apparatus was mounted on a testing machine that was originally built for burst tests on supercharger structures, see Fig. 7.
50 ■
0 -1-1-1-1-1-1-1-1-1-
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
m
Fig. 6. a) Measured engineering stress-strain curves; and b) the resulting true-stress-true-strain curve
Fig. 7. a) An experimental device; b) front view and c) upper view of the experimental specimen
To study the different dynamic behaviours of steel plates with different thicknesses, balls with a diameter of 12 mm and a weight of 7 g were shot at the centre of the specimens at an angle of 20° with different velocities that depended on the specimen thickness, see Table 3. The initial conditions and the plate thicknesses where chosen with regard to the machine's limitations and the expected application of the results (thin-shelled supercharger burst shield). The specimens (see Fig. 7) were metal plates with dimensions of 98 mm x 60 mm. Two different sheet-metal thicknesses were tested: 1 mm and 1.5 mm. The specimens were fixed along the shorter side. The free area of the specimen was 60 mm x 60 mm. The impact velocities were measured on the testing machine just before the impact point using photo-sensors. For each combination of the steel-plate thickness and the impact velocity, three test repetitions were usually performed. The average and the standard deviation of the indentation depth and the position of the indentation centre are presented in Table 3. Some of the plate specimens were also equipped with strain gauges for measuring the strains to obtain the strain rates. Such specimen preparation was time consuming and there were only a few specimens for which the strain gauges did not break off during the measurements. Nevertheless, some strain measurements were successful and an agreement between the measured and simulated strain rates for the 1-mm-thick plate and the impact velocity of 109 m/s was as follows: the measured peak strain rate was approximately 160 s-1, whereas the simulated values at that spot were 200 s-1 to 250 s-1. The major causes for the difference of 30 % were the idealisation of the FE model and the FE mesh resolution, since the strain-rate decay from the impact
point to the border is progressive and it is difficult to calculate the true strain-rate value at the spot of the strain gauges. The simulated peak strain rate at the impact point was 5000 s-1.
Table 3. Combinations of boundary conditions and results from experiments
Experiment Specimen condition thickness number [mm]
Average
ball velocity [m/s]
Measured max.
indentation depth, average, Hexp [mm]
Position of the max. indentation depth, average,
1	1	103	11.37	34.83
2	1	109	12.12	34.89
3	1	121	13.07	35.11
4	1.5	121	10.38	30.19
5	1.5	131	11.53	29.60
6	1.5	139	12.65	30.49
From the results in Table 3 we can see that the scatter of the experimental results is relatively small, which means that the experimental arrangement was appropriate for our study. Furthermore, we can see that the indentation depth increases with the increasing velocity. The indentation depth at a thickness of 1.5 mm is smaller than for the 1-mm-thick steel plate. Besides, it is clear from Table 3 that for the 1.5-mm-thick specimens, the point with the deepest indentation is approximately 30 mm from the lower side of the specimen, just in the centre of the impact. This is not the case for the 1-mm-thick specimens, for which the point with the largest indentation is 5 mm from the centre of the impact. The lack of displacement of the deepest imprint for the thicker plate is a consequence of the fact that after the plastification almost all of the kinetic energy of the ball was consumed. This was not
the case for the thinner plate and the ball proceeded in its direction of travel, though causing an extension of the imprint in the vertical direction. This implies different impact dynamics for the specimens with different thicknesses and this should be replicated by the numerical simulations if the strain-rate-dependent material parameters are properly identified.
3 RESULTS AND DISCUSSION
3.1 FE Model for Identification of the Material Parameters
The LS-DYNA explicit dynamic FE code was used to identify the material parameters of the C-S and J-C material models [16]. The FE model that was used to simulate the ball-impact experiment from Section 2 is presented in Fig. 8. The steel sheet model had 5436 four-node and three-node shell finite elements. The mesh density around the impact area was larger than in the wider region of the specimen model, to accurately simulate the indentation. The mesh density was chosen to optimise the processing time for a reasonable accuracy of the deformation. The rigid ball was modelled with 448 solid finite elements. In the finite-element model, the nodes on both sides of the thin sheet-metal plate were fixed (Fig. 8). A rigid ball was shot into the centre of the sheet at an angle of 20° with different impact velocities, which are presented in Table 3.
Fig. 8. FE model for a ball-plate Impact simulation
Between the flat specimen and the rigid ball was an AUTO_SURFACE_TO_SURFACE contact with the friction coefficient ц = 0.2. This was an approximate average value from different references [24] and [25]. During the simulation, a strain at the left-hand side of the sheet and the gross geometric
dimensions of the specimen were recorded for further processing, like in the experiment (Fig. 7).
The C-S material model was applied using the MAT_PLASTIC_KINEMATIC (MAT3) material model from LS-DYNA ([16] and [17]). This material model is defined with the following parameters, see also Eqs. (1) and (3): material density, elastic modulus, Poisson's ratio, yield stress o-0, tangent modulus Et and the C-S parameters P and C. The parameters Et, P and C were estimated using the procedure from Section 1 and the above-described FE model.
The J-C material model was applied using the MAT_SIMPLIFIED_ JOHNSON_COOK (MAT98) material model from LS-DYNA ([16] and [17]). This material model is defined with the following parameters, see also Eq. (2): material density, elastic modulus, Poisson's ratio, yield stress o-0, and the J-C parameters B, n and c. The latter three parameters were estimated using the procedure from Section 1 and the above-described FE model.
The values of the other material parameters were fixed in our simulation: the material density was 7850 kg/m3, the elastic modulus was 2.1105 N/mm2, the yield stress was o0 = 185 MPa (see Fig. 6) and the Poisson's ratio was v = 0.3. The reference strain rate £0 for the J-C model was 1 s-1 and was taken as a default value from LS-DYNA. This value usually follows from the tensile-test arrangement. However, its choice does not phenomenologically influence the results, since the variation of the parameter £0 monotonically influences the changes of the estimated parameter c:
' "" (1)
C F(1 ) — = ln %
C F( }
2 b0
3.2 Results and Discussion
Table 1 and Fig. 3 represent factor levels for the original domains of the C-S parameters Et, P and C. Table 2 and Fig. 4 represent factor levels for the original domains of the J-C parameters B, n and c. For each of the 81 combinations of the three material parameters for the C-S and J-C material models six simulations were carried out, i.e. three different impact velocities combined with two different specimen thicknesses. The FE simulations were carried out on a numerical server with two Intel Xeon X5670 2.93-GHz processors, 48 GB of RAM and a Linux operating system. The time spent for one numerical simulation on one processor's core was about 2 hours.
The values of the cost function from Eq. (5a) for the C-S material model are presented in Fig. 9. The values of the cost function from Eq. (5b) for
the J-C material model are presented in Fig. 10. In both diagrams the logarithms of the cost functions from Eqs. (5a) and (5b) are presented. The best three combinations of the C-S parameters Et, P and C and the J-C parameters B, n and c are listed in Table 4.
We can see from Fig. 9 and Table 4 that the best combinations of the C-S parameters Et, P and C are in the middle at the top of the original domain. The best combinations of the J-C parameters B, n and c are on the right of the original domains, see Fig. 10 and Table 4. The parameter combinations with the worst cost-function values failed to reproduce, in particular the position of the maximum indentation depth during simulations.
If we look at the results for both material models, we can conclude that the optimal values of the individual parameters can be up to two orders of magnitude distant from each other. This means that the ranges for the original domains of both the C-S and J-C parameter triples were too wide for a reliable estimation of these parameters. For this reason we decided to narrow the ranges of the C-S and J-C parameter domains around their most promising values from Table 4. This was followed by a new simulation plan that was again composed with the help of the L81(910) Taguchi array, but with the Et, P, C and B, n, c parameter levels taken for the narrowed domains in the same manner as was the case for the original domain. The narrowed domains were as follows: (i) for the C-S material model: C = 10 ms-1 to 46.4159 ms-1 P = 3.1623 to 10 and Et = 0.5623 GPa to 1.7783 GPa; (ii) for the J-C material model: B = 0.1778 GPa to 3.1623 GPa, n = 0.1778 to 1 and c = 0.005623 to 0.1778.
Fig. 9. Cost-function values for the original domains of the Cowper-Symonds parameters
Fig. 10. Cost-function values for the original domains of the Johnson-Cook parameters
Fig. 11. Cost-function values for the narrowed domains of the Cowper-Symonds parameters
Fig. 12. Cost-function values for the narrowed domains of the Johnson-Cook parameters
We actually performed a "nested" Taguchi simulation plan in the second phase of the simulations.
The resulting distributions of cost-function values for the narrowed domains are presented in Fig. 11 for the C-S parameters and in Fig. 12 for the J-C parameters. The best three combinations of the C-S parameters Et, P and C and the J-C parameters B, n and c for the narrowed domains are much closer together when compared to the original domains for those parameters. The average values from the three best solutions are listed in Table 5. We can conclude from Table 5 that our estimations of the parameters Et, P, C and B, n, c are comparable to the values reported for mild steels in the literature ([9] to [11]), despite different experimental arrangements and the fact that the strain rates during the experiments were 10 or more times higher in our case. We can conclude that we obtained reasonable estimates of the parameters Et, P, C and B, n, c for our case.
Plastic flow curves a - speff for the two material models are presented in Fig. 13. They were calculated for different strain rates with the averaged parameters E, P, C and B, n, c from Table 5. We can see from this figure that the estimates of the parameters for both material models were consistent, because the flow curves span a similar domain of the a - speff space for the two material models.
From the results we can conclude that the described methodology, which combines the nested design with the FE simulations, can be very time
efficient for estimating the parameters of material models that govern the material's behaviour at high strain rates. The added value of our approach is meaningful, especially in the cases when the number of parameters that need to be identified is relatively high, with a wide range of potential parameter values. With the Taguchi orthogonal array, a reasonable estimate of the material parameters can be found with a relatively small computing effort and a short time. For example, if we applied a methodology that is based on genetic algorithms and was originally developed for estimating the foam-material-model parameters [18], the processing times would be approximately 100-times longer. This would be appreciated if the numbers of parameters to be identified and the wide ranges were to be increased.
4 CONCLUSIONS
The article presents a general approach to estimating the parameters that govern a material's behaviour at high strain rates. In our approach the Taguchi experimental design was combined with the FE code LS-DYNA to estimate the material parameters using the results of the impact test between the ball and thin sheet metal. The presented approach was applied to the realistic case of a material-parameter estimation for two different material models, i.e., the C-S and the J-C material models.
Table 4. The best three combinations for the original domains of the C-S and J-C material parameters
		Parameter C [ms-1]	Parameter P [/]	Parameter Et [GPa]	Cost-function value [/]
Cowper-Symonds	Combination 1	21.5443	5.6234	1.0000	1.524
material model	Combination 2	46.4159	10.000	1.0000	1.588
	Combination 3	2.1544	3.1623	1.0000	1.826
		Parameter B [GPa]	Parameter n [/]	Parameter c [/]	Cost-function value [/]
Johnson-Cook	Combination 1	0.1778	0.1778	1.0000	2.020
material model	Combination 2	0.5623	0.4217	0.0056	2.259
Combination 3	1.7783	1.0000	0.0056	2.129
Table 5. The average of the best three combinations for the narrowed domains of the C-S and J-C material parameters
		Parameter C [ms-1]	Parameter P [/]	Parameter Et [GPa]	Cost-function value [/]
Cowper-Symonds material model	Our average	41.0133	6.2000	0.9550	1.522
	Belingardi et al. [12]	3.006-4.987	1.329-1.619	-	-
	Marais et al. [10]	2.000	5.000	-	-
	Markiewicz et al. [11]	1.150	7.750	-	-
		Parameter B [GPa]	Parameter n [/]	Parameter c [/]	Cost-function value [/]
Johnson-Cook	Our average	1.9250	0.8183	0.0972	2.138
material model	Singh et al. [9]	0.779-2.692	0.743-0.928	0.0144-0.021	-
	Marais et al. [10]	0.292	0.310	0.025	-
It turned out that it is possible to obtain reliable estimates of the C-S parameters (Et, P, C) and J-C parameters (B, n, c) with a nested design-of-simulation approach using only two iteration runs. The materialparameter estimates for the two models are consistent and comparable to the results from the literature.
Fig. 13. Plastic flow curves for different strain rates (in s-1) for a) the Cowper-Symonds, and b) the Johnson-Cook material models
5 REFERENCES
[1]	Meyers, M.A. (1994). Dynamic Behavior of Materials. John Wiley & Sons, New York, D0I:10.1002/9780470172278.
[2]	Armstrong, R.W., Walley, S.M. (2008). High strain rate properties of metals and alloys. International Materials Reviews, vol. 53, no. 3, p. 105-128, D0I:10.1179/174328008X277795.
[3]	Johnson, G.R., Cook, W.H. (1983). A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures. Proceedings of the 7th International Symposium on Ballistics.
[4]	El-Magd, E. (1994). Mechanical properties at high strain rates. Le Journal de Physique IV, vol. 4, no. C8, p. 149-170, D0I:10.1051/jp4:1994823.
[5]	Zhao, H., Gary, G. (1996). The testing and behaviour modelling of sheet metals at strain rates from 10-4 to 104 s-1. Materials Science and Engineering: A, vol. 207, no. 1, p. 46-50, D0I:10.1016/0921-5093(95)10017-2.
[6]	Huh, H., Kim, S.B., Song, J.H., Yoon, J.H., Lim, J.H., Park, S.H. (2009). Investigation of elongation at fracture in a high speed sheet metal forming process. Steel Research International, vol. 80, no. 5, p. 316-322.
[7]	Huh, H., Lim, J.H., Park, S.H. (2009). High speed tensile test of steel sheets for the stress-strain curve at the intermediate strain rate. International Journal of Automotive Technology, vol. 10, no. 2, p. 195-204, D0I:10.1007/s12239-009-0023-3.
[8]	Kim, S.B., Hoon, H. (2008). Evaluation of the failure elongation of steel sheets for an auto-body at the intermediate strain rate. Key Engineering Materials, vol. 385-387, p. 749-752, DOI:10.4028/www.scientific.net/kem.385-387.749.
[9]	Singh, N.K., Cadoni, E., Singha, M.K., Gupta, N.K. (2013). Mechanical behavior of a structural steel at different rates of loadings. Proceedings of the International Symposium on Engineering under Uncertainty: Safety Assessment and Management, p. 859-868, D0I:10.1007/978-81-322-0757-3_57.
[10]	Marais, S.T., Tait, R.B., Cloete, T.J., Nurick, G.N. (2004). Material testing at high strain rate using the split Hopkinson pressure bar. Latin American Journal of Solids and Structures, vol. 1, no. 3, p. 219-339.
[11]	Markiewicz, E., Ducrocq, P., Drazetic, P. (1998). An inverse approach to determine the constitutive model parameters from axial crushing of thin-walled square tubes. International Journal of Impact Engineering, vol. 21, no. 6, p. 433-449, D0I:10.1016/S0734-743X(98)00004-9.
[12]	Belingardi, G., Chiandussi, G., Ibba, A. (2006). Identification of strain-rate sensitivity parameters of steel sheet by genetic algorithm optimisation. Brebbia, C.A. (ed), High Performance Structures and Materials III. WIT Press, Wessex Institute of Technology, p. 201-210, DOI:10.2495/HPSM06021.
[13]	Rule, W.K. (1997). A numerical scheme for extracting strength model coefficients from Taylor test data. International Journal of Impact Engineering, vol. 19, no. 9-10, p. 797-810, D0I:10.1016/S0734-743X(97)00015-8.
[14]	Kurtaran H., Buyuk, M., Eskandarian, A., (2003). Ballistic impact simulation of GT model vehicle door using finite element method. Theoretical and Applied Fracture Mechanics, vol. 40, no. 2, p. 113-121, D0I:10.1016/S0167-8442(03)00039-9.
[15]	Schwer, L.E., Hacker, K., Poe, K. (2006). Perforation of metal plates: laboratory experiments and numerical simulations. Proceedings to the 9th Annual LS DYNA Users Conference.
[16]	Hallquist, J.O. (1998). LS-DYNA Theoretical Manual, Livermore Software Technology Corporation, Livermore.
[17]	Hallquist, J.O. (2007). LS-DYNA Keyword User's Manual -Version 971, Livermore Software Technology Corporation, Livermore.
[18]	Škrlec, A., Klemenc, J., Fajdiga, M. (2014). Parameter identification for a low-density-foam material model using numerical optimisation procedures. Engineering Computations, vol. 31, no. 7, p. 1532-1549, DOI:10.1108/EC-03-2013-0100.
[19]	Taguchi, G. (1987). System of Experimental Design: Engineering Methods to Optimize Quality and Minimize Costs, Vol. 1. American Supplier Institute, Dearborn.
[20]	Taguchi, G. (1987). System of Experimental Design: Engineering Methods to Optimize Quality and Minimize Costs, Vol. 2. American Supplier Institute, Dearborn.
[21]	ASTM E8/E8M-09. 2010. Standard test methods for tension testing of metallic materials. ASTM International. West Conshohocken, D0l:10.1520/E0008_E0008M-09.
[22]	Dowling, N.E. (1993). Mechanical behavior of materials: Engineering Methods for Deformation, Fracture, and Fatigue. Prentice-Hall Inc., Englewood Cliffs.
[23]	ASTM D5420-10. 2010. Standard Test Method for Impact Resistance of Flat, Rigid Plastic Specimen by Means of a Striker Impacted by a Falling Weight (Gardner Impact). ASTM International. West Conshohocken, D0I:10.1520/D5420-10.
[24]	Wittel, H. (2009). Roloff/Matek Maschinenelemente: Tabellenbuch, Vieweg&Teubner, Wiesbaden.
[25]	Klemenc, J., Wagner, A., Grubisic, V., Fajdiga, M. (2011). Schienenzustand und Gleitreibungszahl zwischen Rad und Schiene. ETR, vol. 60, no. 11, p. 34-38.
Strojniški vestnik - Journal of Mechanical Engineering 62(2016)4, 231-242 © 2016 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2015.3192	Original Scientific Paper
Flow Unsteadiness and Pressure Pulsations in a Nuclear Reactor Coolant Pump
Dan Ni* - Minguan Yang - Bo Gao - Ning Zhang - Zhong Li
Jiangsu University, School of Energy and Power Engineering, China
Unsteady flow induced by rotor-stator interaction is detrimental to the safe operation of the nuclear reactor coolant pump, so it is essential to clarify flow structures and pressure pulsations in such pumps, especially within the spherical casing. In this paper, unsteady flow characteristics in a mixed-flow nuclear reactor coolant model pump were investigated using large-eddy simulation (LES) method. Results show that at the nominal flow rate, in two particular diffuser channels near the spherical casing discharge nozzle, the flow structures are uneven compared with that in the other flow channels. The reason is associated with the position of the flow channel with respect to the spherical casing nozzle. Large- scale flow separation and backflow structures easily occur at the regions near these two channels. In the right and the middle region of the casing nozzle, due to the large-scale separate flow and high vorticity magnitude, unsteady flow structures are more complicated in comparison with the other regions. It has been found that the vorticity spectra and the pressure spectra almost have the same main excitation frequencies. Therefore, it has been confirmed that for particular regions pressure pulsations are determined by the shedding vortex wake from the diffuser blade trailing edge of the nuclear reactor coolant model pump.
Keywords: nuclear reactor coolant pump, flow unsteadiness, pressure pulsation, vorticity magnitude Highlights
•	Numerical investigation of flow unsteadiness and pressure pulsation in a nuclear reactor coolant pump.
•	Flow separations and backflows easily occur at the diffuser channels near the discharge nozzle region.
•	At the left region of the casing, the discharge nozzle is affected significantly by an intense rotor-stator interaction effect.
•	In the right and middle of the casing nozzle, unsteady flow structures are more complicated.
•	An unsteady vortex-shedding effect would motivate evident component in pressure spectrum.
0 INTRODUCTION
A nuclear reactor coolant pump (RCP) [1], the device in the reactor core and the steam generator that transfers the heat energy [2], is the "heart" of a nuclear reactor driving the circulation flow of the coolant in the main loop [3]. The nuclear reactor coolant pump is the only rotating part of the nuclear island, so it belongs in a nuclear safety grade one facility. Moreover, it is also the main energy-consuming equipment in the nuclear power plants, which must be guaranteed to operate continuously over a long term and without trouble. An investigation on the internal flow in the mixed-flow nuclear reactor coolant pump is one of the key problems for the development of reactor coolant pumps in large pressurized water reactors. Some investigations have been conducted to study the design of the mixed-flow pump [4] and [5] and flows instabilities [6] and [7]. However, most of them only focus on the unsteady pressure pulsation in impeller [8] and [9] or the performance of the nuclear reactor coolant pump [10] and [11]. The unsteady flow structure of a mixed-flow nuclear reactor coolant pump, especially in certain specific regions, is very important to the safety analysis of a nuclear reactor, as it could generate serious flow-induced vibration, threatening the pump integrity [8]. Therefore, a
comprehensive analysis and prediction of pressure pulsation caused by intense rotor-stator interaction are essential for the design of the mixed-flow nuclear reactor coolant pump [12] and [13]. At present, in order to improve the efficiency and stable operation of the nuclear reactor coolant pump, the complicated internal unsteady flow structures should be thoroughly illustrated.
The nuclear reactor coolant pump has a special structure equipped with a spherical casing, which determines a typical and complex flow pattern within the pump. Knierim et al. [14] designed a new type of reactor coolant pump for a 1400 MW nuclear power plant. The impeller and diffuser were gradually optimized based on the computational fluid dynamics (CFD). The volute casing is a spherical shape with a discharge nozzle facing the impeller, and the flow structures in the region around the discharge nozzle are uniform. The region of the discharge nozzle itself is characterized by the fact that the flow below the outlet port is divided into two parts. One portion flows out of the casing discharge nozzle, whereas the other portion circulates around the casing once more prior to exiting. Kato et al. [15] and [16] described internal flows of a high-specific-speed mixed-flow pump at low flow rates using large-eddy simulation (LES), and it showed that the head-flow curve exhibited weak
instability characteristics. Posa et al. [17] reported the LES method used in a mixed-flow pump, where a structured cylindrical coordinate solver with optimal conservation properties was utilized in conjunction with an immersed-boundary method. Moreover, the overall agreement with the experimental results was excellent, demonstrating the robustness and feasibility of the approach in rotating machinery applications. It demonstrated that the LES method was adequate to predict complicated flow patterns of a high-specific-speed mixed-flow pump. With the development of non-contact measuring techniques, unsteady particle image velocimetry (PIV) [18] and laser Doppler velocimetry (LDV), measuring techniques are often applied to investigate complex unsteady internal flow in pumps, in such a way that no external disturbance is imposed on the flow field. Miyabe et al. [6] and [7] used PIV and pressure fluctuation measurements to investigate the propagation mechanism of a rotating stall in a mixed-flow pump. It was found that unstable performance was caused by periodic large-scale abrupt backflow generated from the diffuser to the impeller outlet. However, most research focuses on the design and the instability flow of the impeller and diffuser in a mixed-flow pump [19] to [22], and unsteady flows in a mixed-flow nuclear RCP with the specific spherical casing are rarely conducted. Consequently, true internal flow structures have not been thoroughly revealed.
In this study, based on the LES method, the internal unsteady flow in a mixed-flow RCP model pump is studied. Monitoring points are mounted on the impeller inlet, diffuser channels, and spherical casing to achieve pressure pulsation characteristics. Finally, pressure pulsations and complicated flow structures are combined together to clarify flow unsteadiness in the RCP model pump.
1 NUMERICAL METHOD 1.1 Governing Equations
By neglecting the incompressible effect of the fluid, LES governing equations [15] to [17] can be described in the following form.
u / \ —P ) = 0,
dx
(1)
д д д — (piii ) +-(Puiuj ) =-
St
dx
dx
да i—. дх
Л
дР-дт, (2)
дх дх
where ut (i = 1, 2, 3) is the grid-scale velocity component, p is the grid-scale static pressure, p is
the density and ß is the kinematic viscosity, and at] is the stress tensor having the form of Eq. (3).
f dU. ды.^ —- + —-
dx
v J
dx
2	ды.
-V — 5.
3	dx '
(3)
where Si]- is Kronecker delta. Ti]- is the sub-grid-scale (SGS) stress tensor, which is defined in Eq.(4).
t = pu - puuj.	(4)
In this study, the Smagorinsky-Lilly SGS fixed-coefficient model [23] is applied to close the equations. In this model, the eddy viscosity is defined in Eq. (5).
= PLs S :
(5)
where LS is the mixed length of the grid and
2
|S| = \l2SijSij , where Si]- is the deformation tensor.
The function of LS is defined in Eq. (6).
LS = min (kd, CSV "'),
(6)
where к is the von-Karman constant, d is the distance to the nearest wall, V is the volume of the computational cells and Cs is the Smagorinsky constant. Cs is usually a constant 0.1.
1.2 Mesh Generation
A mixed-flow nuclear reactor coolant model pump was first designed. The main parameters of the RCP model pump are shown in Table 1. The RCP mainly consists of an impeller, diffuser, and spherical casing. The structure of the RCP model pump is illustrated in Fig. 1.
The quality of the mesh has a critical influence on the numerical simulation accuracy. The Ansys-ICEM mesh generation tool was adopted to generate structured mesh cells of the model pump in order to obtain precise unsteady flow structures and pressure fluctuation characteristics. To guarantee the computational accuracy, the fine grids were necessarily required by the LES model. Thus, grid cells near the solid wall were encrypted, especially on the surface of the blades, where significant pressure gradients and flow separation may easily occur [12] and [13]. Meanwhile, a detailed grid sensitivity analysis with four different mesh numbers has been carried out, as summarized in Table 2. After mesh sensitivity checks, it is noted that the pump head under nominal flow rate decrease 0.2 % with the mesh number increasing and, finally, the value of the pump head is constant. Finally, Case 3 of a total mesh number nearly 7*106 is selected for numerical simulation considering the
current computational ability. The numerical results demonstrate that at the near wall region of the impeller and the diffuser, the y+ value is about 1. In the other region of the impeller and the diffuser, the y+ value is lower than 5. Finally, the averaged y+ value of the whole computational domain of the impeller and the diffuser is approximately 4.5, which could provide adequate resolution in the critical regions of the computation domain. Fig. 2 shows the mesh of the RCP model pump impeller.
Table 1. Main design parameters
Parameters			Values	
Nominal flow rate QN			848 m3/h	= 0.236 m3/s
Designed head HN			12.7 m	
Rotating speed n			1480 r/min	
Specific speed ns	= (3.65«VQN )/ HN 075			390
Impeller outlet diameter d2			268 mm	
Impeller outlet width b2			84 mm	
Casing diameter D			637.5 mm	
Impeller blade number Z				4
Diffuser blade number Zd				12
Table 2. Mesh sensitivity check				
Parameters	Case 1	Case 2	Case 3	Case 4
Inlet section	136,567	257,863	336,835	461,823
Impeller	1,237,651	2,210,789	3,231,400	4,247,913
Diffuser	884,132	1,007,541	1,230,624	1,545,004
Spherical casing	1,060,434	1,735,866	2,137,419	3,845,103
Outlet section	184,471	360,547	502,281	651,079
Total elements	3,503,255	5,572,606	7,438,559	10,750,922
Pump head in design point [m]	12.405	12.378	12.369	12.365
1.3 Solution Parameters
The LES method [26] and [27] has a better ability to predict unsteady flow than the RANS method in general [28] and [29], especially for the flow separation [30] to [32]. Before the calculation of transient LES simulation, the steady calculation should be achieved. In this study, steady calculation using the standard ks model is first conducted, and the results are set as the initial boundary for transient LES simulation. The pump inlet is defined as the velocity inlet boundary condition and a constant pressure p = 101325 Pa is imposed at the pump outlet. The coupling of the rotational domain (impeller) and the stationary domain (diffuser) is solved by moving mesh.
Fig. 2. The mesh of the impeller and its partially refined structured grids
Fig. 1. Diagram of the mixed-flow RCP model pump
In order to illustrate the flow unsteadiness characteristics of the RCP model pump, especially within the spherical casing, time step size is set as At = 112.6 ^s in order to have enough resolution of unsteady flow results during calculation. Thus, each impeller revolution will be calculated in a time sequence of 360 time steps corresponding to 1° of the impeller rotation speed. The numerical residual convergence criterion is set as 10-5 in order to ensure the result to be converged. The unsteady convergence also has a great influence on the achieved results, so, at least, fifty revolutions were achieved to guarantee the numerical accuracy.
a)
c)
Fig. 3. Details of the monitoring points on the RCP model pump; a) impeller inlet, b) diffuser channels, and c) discharge nozzle of spherical casing
1.4 Monitoring Points
Monitoring points are mounted on the principal districts to obtain the comprehensive understanding of flow unsteadiness characteristics. Many researchers have already investigated the unsteady characteristics within impeller channels [9], so the emphasis is laid upon pressure pulsation within the diffuser channels, the spherical casing, and the pump inlet. Therefore, five monitoring points are mounted on the impeller inlet as seen in Fig. 3a. Fig. 3b shows monitoring points in six diffuser channels out of a total of 12 flow channels. All the points are locating at the mid-span plane of the diffuser. The diffuser channel on three monitoring points (such as D-A1, D-A2, and
D-A3) is defined as Channel A. The Channel B and Channel C have similar definitions. It is well known that unsteady flow distribution in the spherical casing of the RCP is quite complicated, especially at the area near the discharge nozzle of casing [14]. Thus, monitoring points are mounted on the mid-span plane of the spherical casing nearing the discharge nozzle zone as shown in Fig. 3c. As discussed in the preceding, the complicated flow pattern is expected in this region. Therefore, a total of 18 monitoring points are selected for investigation in this region. These points are divided into three groups for later analysis as shown in Fig. 3c.
2 EXPERIMENTAL SETUP
In order to validate the accuracy of the current numerical method, experimental investigation of the RCP model pump is carried out in a closed-loop test rig to guarantee the measuring accuracy. Details of the test loop are shown in Fig. 4. The water temperature is about 25 °C during experiments. Flow rates of the RCP at various conditions are measured with an electronic flow meter with an absolute accuracy of
Fig. 4. Closed test loop
Fig. 5. Test used impeller
±0.2% of the measured value. Two pressure gauges at the pump inlet and outlet pipings are used to obtain the pump head with an accuracy of ±1%. An auxiliary pump is used to overcome the pressure drop in the system. The rotating speed of the RCP model pump was ensured to be at the design value 1480 r/min by adopting a frequency inverter during the experimental process. Fig. 5 depicts the test of the used impeller.
3 RESULTS AND DISCUSSIONS
3.1 Validation of the Numerical Method
In order to validate the reliability of CFD analysis, CFD results are compared with experimental results of the RCP model pump. Fig. 6 shows the comparisons between experiments and the numerical results obtained from transient LES models. As observed, the maximum difference from the comparison is about 3.5 %. At nominal flow rates, the difference between the predicted and measured values is about 2.4 %. The discrepancy between CFD and experimental results may be attributed to the grid resolution around the blade surface, and the current mesh resolution may be insufficient to resolve the turbulent boundary layer developing on the blade surface. Furthermore, effects of the leakage flow through the clearance are not considered in the present numerical simulation [15]. However, it can be found that the tendency of CFD results has a good agreement with the test results. From performance comparison, it can be concluded that the current numerical method could capture the mainly unsteady flow structures of the RCP model pump.
3.2 Unsteady Pressure Pulsation Characteristics at the Impeller Inlet
In order to analyse the unsteady pressure pulsations, pressure coefficient Cp is defined in Eq. (7).
Cp =AP/ (0.5pu2),
nnD2
u2 =-2,
2 60
AC = C - C
p	pmax pi
(7)
(8) (9)
where AP is pressure amplitude, u2 is tangential velocity at the impeller outlet, and D2 is the impeller outlet diameter. ACp is the difference value, which is the maximum of Cp minus minimum of Cp.
cp 0.010
a)
0.006 0.004 0.01)2 0.0U0 0.006 0.004 0.002 0.000
O. 0.004 0.002 0.000 0.006 0.004 0.002 0.000 0.0012
./врг
Jfm-
100 г , 1	»1 200	4 , i 300 400	i . .r . i 500 600 700 -A2
		It . i	
100	200	300 400 it i	500 600 700 -A3
100 1	200	300 400 Л , J	500 600 700 -A4
100	200 i	300 400 1 . ,	500 600 700 -A5 ■i.l,.,
Fig. 6. Performance comparison between experimental and numerical results
100	200	300	400	500	600	700
tjj	Frequency [Hz]
Fig. 7. a) Pressure coefficients of the monitoring points at the impeller inlet, and b) frequency domain of each monitoring point at impeller inlet
The rotating speed of the impeller keeps a constant value of 1480 r/min, so the impeller-rotating frequency (fR) is 24.7 Hz and the blade-passing frequency
Obpf) is 99 Hz. Fig. 7a shows pressure coefficients of the monitoring points in the impeller inlet during one impeller-rotating period. It is observed that the pressure signal of the points fluctuates significantly with the impeller rotating, and pressure coefficients of the monitoring points (A1~A4) have an approximately symmetrical distribution in one impeller rotating cycle. Differences of pressure coefficient among these four monitoring points are quite small. It indicates that pressure signals in the periphery of the impeller inlet are nearly consistent with the impeller rotating. However, the point A5 at impeller inlet centre shows different pressure coefficient distribution compared with the other four measuring points, and the pressure coefficient value is lower than others. It does not show any obvious periodicity in one rotating cycle. The fast Fourier transform (FFT) method was adopted to change the time domain signals into the frequency domain signals to analyse pressure spectra. The hanning window was adopted to reduce the energy leakage during signal process [33]. Fig. 7b shows pressure spectra of the monitoring points at the impeller inlet. The circumferential monitoring points (A1~A4) have similar pressure spectra, and the dominant peaks occur at the blade-passing frequency 0BPF). Moreover, high harmonics at 297 Hz (3 fBPF) could also be identified with an amplitude about 1/6 of that at fBPF. However, the predominant component of the point in the impeller inlet centre (A5) occurs at three times the blade-passing frequency (3 fBPF).
3.3 Unsteady Pressure Pulsation Characteristics within the Diffuser Channels
Fig. 8 presents pressure coefficient of monitoring points within three diffuser channels near the discharge nozzle in one rotating cycle. It is evident that pressure signals change periodically, and four crests and troughs exist during one impeller revolution due to rotor-stator interaction. From the diffuser inlet to the outlet, when the pressure coefficients present wave crests, the pressure amplitudes gradually decrease and when they present wave troughs, the amplitude increases gradually in one rotating cycle. Obviously, near the diffuser channel inlet, the pressure coefficients fluctuate more greatly than the diffuser outlet. It is indicated that rotor-stator interaction effect is obvious, and it is weakened gradually along with the diffuser channels. Meanwhile, pressure coefficients of these three diffuser channels have different distributions.
In order to quantitatively analyse this phenomenon, the difference values of pressure
a)
b)
c)
Fig. 8. Pressure coefficient in diffuser: a) Channel A, b) Channel B, and c) Channel C
coefficients are analysed on the diffuser monitoring points in a rotating cycle, as shown in Fig. 9. It is
indicated that diffusers Channel C and Channel B have relatively large difference values of the pressure coefficient, especially for Channel C. These two flow channels are near the casing discharge nozzle. Thus, it is concluded that in a rotating cycle in these two flow channels, the flow structure is a more uneven compared with other flow channels. The reason is associated with the location of the flow channel, and larger flow separations or backflows occur easily near the two flow channels.
Fig. 9. Pressure coefficients difference values on diffuser monitoring point
Fig. 10. Pressure pulsation amplitudes at blade-passing frequency (fBPF) and higher harmonic
Fig. 10 shows the pressure pulsation amplitude at the blade-passing frequency (fBPF) and its high harmonic of monitoring points in diffuser channels
(Channels A~C). Obviously, along with the flow channels, the blade-passing frequency and higher harmonic become increasingly weak, and the main excitation frequency in the diffuser channel is the blade-passing frequency fePF). Pressure pulsation amplitude values at the blade-passing frequency (fBPF) and its second harmonic (2 fBPF) in Channel B are larger than the other two channels. Channel C has the second largest amplitude value, and the pressure pulsation amplitude value is minimum in Channel A. However, at higher harmonics (3 fBPF and 4 fBPF), there are some different features, especially in Channel B. At 3 fBPF, the amplitude values of Point D-B3 and Point D-B2 in diffuser Channel B are lower than the other channels. Meanwhile, at 4 fBPF, the amplitude of Point D-B1 is less than Point D-A1 and Point D-C1. That is because some other complex unsteadiness flows in flow Channel B and the rotor-stator interaction effects become weaker. Thus, it is indicated that the blade-passing frequency due to the rotor-stator interaction effect is obviously different in every diffuser channel. Furthermore, internal flows in different diffuser channels, especially near the spherical casing discharge nozzle region, are complicated.
3.4 Unsteady Pressure Pulsation Characteristics within Spherical Casing
Due to the particular spherical casing, flow structures show a greater difference in comparison with the conventional spiral volute. From this perspective, it is important to analyse unsteady flow characteristics within the spherical casing.
In order to analyse the unsteady pressure spectra in the spherical casing, Figs. 11 to 13 show the pressure pulsation characteristics around the discharge nozzle region. Monitoring points at the discharge nozzle region are divided into three groups symmetrically. Region 1 is located on the left of the discharge nozzle. Region 2 is located in the right of the discharge nozzle, and Region 3 is located in the middle of the discharge nozzle. The vertical scales in Figs. 11 to 13 are different, due to the increasing role of the low frequencies components.
Fig. 11 presents the pressure spectra in Region 1. From the left side of the casing to the middle of the casing, pressure amplitudes at the blade-passing frequency fBPF) and its high harmonic (2fBPF) decrease gradually. Near the wall surface, these monitoring points (C1, C2, C8, and C13) have the main excitation frequency, which is the blade-passing frequency fBPF). It is indicated that the rotor-stator interaction effect is
o.ooo 0.002
50	100	150	200	250	300	350	400
^10 Hz							C9
UVwVvv	vjlv		- I\				
50	100	150	200	250	300	350	400
iy,\3 Hz	1						C13
"^AWMA/1	lAjJflA		Л				
50	100	150	200	250	300	350	400
,/fio Hz f\							C14
	jJW	\ллЛ~-.—Л _					
150 200 250
Frequency [Hz]
Fig. 11. Pressure spectra of monitoring points in Region 1
Fig. 12. Pressure spectra of monitoring points in Region 2
significant. However, the monitoring point C3 has a complex pressure spectrum as indicated by many peaks appearing. As observed, the dominant excitation frequency is 21 Hz, and the second larger peak is 46 Hz. The interval of frequencies is about 20 Hz to 50 Hz in the low-frequency band. These low-frequency signals are related to the complicated internal flow structure developing at the casing discharge nozzle region. For C9, C13, and C14, significant excitation frequencies occur within the low-frequency band, i.e. 10 Hz, 13 Hz, and 10 Hz respectively. The interval of the main frequencies is approximately 10 Hz to 15 Hz. Therefore, from analysis of the pressure spectra characteristics in Region 1, it is concluded that pressure pulsations of the points at the left of the casing discharge nozzle are affected significantly by an intense rotor-stator interaction effect. However, for the points near the middle of the casing discharge nozzle, the rotor-stator interaction effect is less significant, and the dominant component mainly appears in the low-frequency band.
0.006 0.004 0.002 0.000 0.006 0.004 0.002 0.000 0.006 (1.004 0.002 0.000 0.006 0.004 0,002 0.000

-C4

2 fn pi
50	100	150	200	250	300	350	400
• у, 18 Hz						-СЮ	
50 Г /,24 Hz	100	150	200	250	300	350 -C15	400
yltyw							
50 Г 24 Hz	100	150	200	250	300	350 -CIS	400
fVv	■Л/Wj	VUVA-					
50	100 150 200 250 300 350 400
Frequency [Hz]
Fig. 13. Pressure spectra of monitoring points in Region 3
Pressure spectra of monitoring points in Region 2 are presented in Fig. 12. It is noted that pressure spectrum characteristics of the points in Region 2 are more complex compared with that in Region 1. Components at the blade-passing frequency and its second harmonic are weaker. As observed, the main
excitation frequencies of different monitoring points differ obviously, namely 21 Hz for C5 and C6, 19 Hz for C7 and C17, 18 Hz for C11, 28 Hz for C12 and 13 Hz for C16. The interval of the main frequencies can be defined between 10 Hz and 50 Hz. Therefore, it can be concluded from pressure spectrum analysis, for the monitoring points at the right region of the casing discharge nozzle, that the rotor-stator interaction is not the dominant source of pressure fluctuations. In this region, many components in the low-frequency band are captured, which are related to the complex flow structures developing in this area.
Fig. 13 presents pressure spectra of the monitoring points in Region 3. As noted, components at the blade-passing frequency almost could not be observed, so the rotor-stator interaction effect in this region is also not significant. In the low-frequency band, the main excitation frequency of the monitoring point of C4 is 13 Hz. For C10, evident peaks occur at 18 Hz and 24 Hz for C15, C18. In this region, the interval of the main frequencies is 10 Hz to 50 Hz.
Table 3. Main excitation frequency of each monitoring point
Monitoring points	Main excitation frequency [Hz]	Main frequencies interval [Hz]
C1, C2, C8, C13	/bpf = 99	fBPF
C3, C5, C6	21, 46	20 to 50
C10, C11	18	10 to 20
C9, C14	10	10 to 15
C15, C18	24	20 to 50
C7, C17	19	10 to 20
C4, C16	13	10 to 20
C12	28	20 to 50
Fig. 14. Distributions of static pressure at the centre plane marked with the main excitation frequencies
The main excitation frequency of each monitoring point is concluded in Table 3. Fig. 14 gives the details of the static pressure distribution at the centre plane of the casing. From Table 3 and Fig. 14, it can be concluded that the pressure spectra of the points near the left of the casing nozzle are dominated by rotor-stator interaction. And in this region, flow separation does not occur. However, near the middle and right of the spherical casing discharge nozzle, several dominant peaks of low frequencies occur in pressure spectra. Although these frequencies are very complicated, there are some connections between these low frequencies. Near the outlet of the diffuser Channels A and B, the same frequency at about 21 Hz occurs. It is associated with complex flow structures at this region, including the vortex shedding from the diffuser blade exit and even backflow near the outlet of Channels A and B. However, for the points at the diffuser blade outlet between Channels B and C, a peak at 13 Hz dominates the pressure spectrum. This excitation frequency is caused by the shedding wake effect from the blade outlet. Near the right side of the spherical casing discharge nozzle, some different frequencies are found. Because of the highvelocity flow, unsteadiness fluid from Channels B and C impinges on the wall, and some complicated flow structures including large-scale flow separation and backflow are generated in this region, causing the main excitation frequencies.
From the above analysis, it is inferred that unsteady flow structures in diffuser channels and spherical casing discharge nozzle have a huge impact on pressure pulsations of the RCp model pump. Fig. 15 shows typical distributions of velocity for one moment at the centre plane of the casing. It is found that two typical backflows occur in Region 2. One can be defined as the run-through flow, which is located at the diffuser Channel B outlet, and it flows out of the casing discharge nozzle. And the other backflow can be defined as the circulating flow, which is located at the diffuser Channel C outlet, and it flows around the casing once more prior to exiting. Fig. 16 gives the typical instantaneous distributions of vorticity magnitude at the centre plane of the RCP. It can be clearly seen that a strong vorticity magnitude region near the diffuser Blade B and diffuser Blade C occur within the casing. Thus, in Region 2, due to the large-scale flow separation (the run-through flow and the circulating flow) and high vorticity magnitude, unsteady flow structures are more complicated than that in the other region.
In order to find the reason behind the complicated pressure pulsations in Region 2, Fig. 17 gives the
Fig. 15. Typical velocity distribution at the centre plane
t = T+At	t = T+50At
Fig. 16. Typical instantaneous distributions of vorticity magnitude at the mid-span plane of the casing
vorticity spectra of some monitoring points. Combined with Figs. 11 to 13 and Table 3, it is found that the vorticity magnitude and pressure coefficient have the same main excitation: about 13 Hz for C4, 21 Hz and 46 Hz for C5, 46 Hz for C6. Thus, it indicates that an unsteady vortex shedding effect would motivate evident component in pressure spectrum. In summary, in the right and middle of the casing discharge nozzle (Region 2 and Region 3), the internal flow structures are more complicated, as manifested by the large-scale flow separation, the backflow and the high magnitude vortex shedding from the diffuser exit in comparison with that in Region 1.
-C4
■e a
"S.60
50 100	150	Till- -Ir 1 200	250	300	350	400
13 Hz						
f 21 Hz						
It .^f 46 Hz					-C5	
						
50 100	150	200	250	300	350	400
13 Hz						
:r| bj46Hz					-C6	
50 100	150	200	250	300	350	400
г я 24 Hz						
L t					--C7	
-Vyvt/VuVWVA,						
SO	100	150 200 250 300 350 400
Frequency [Hz]
Fig. 17. Vorticity spectra of particular monitoring points
4 CONCLUSIONS
Based on the LES method, the internal unsteady flow in a mixed-flow nuclear RCP at the rated condition is studied in this paper. Some conclusions are made.
Due to the particularity of the spherical casing in the RCP, the internal flow structure is especially ambiguous and complicated. Research results show that at the nominal flow rate, flow structures differ significantly in different diffuser channels, and it is closely associated with the position of the diffuser channel related to the casing nozzle. At the diffuser channels near the nozzle region, flow separations and backflows easily occur.
At the left region of the casing discharge nozzle is affected significantly by an intense rotor-stator interaction effect. However, in the right and middle of the casing nozzle, due to the large-scale separated flow (the run-through flow and the circulating flow) and high vorticity magnitude, unsteady flow structures are more complicated in comparison with the other regions. Moreover, vorticity and pressure spectra almost have the same main excitation frequency. It is believed that some peaks in pressure spectra are caused by the shedding vortex wake from the diffuser blade outlet for some specific regions.
In further study, the investigation would be concentrated on the relationship between the unsteady flow structures and pressure pulsation characteristics using PIV and pressure pulsation measuring techniques in the mixed-flow nuclear reactor coolant model pump for all the concerned conditions.
5 ACKNOWLEDGMENT
The authors gratefully acknowledge the financial support of National Natural Science Foundation of China (51476070, 51206063, 51576090), a Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).
6 REFERENCES
[1]	Gao, H., Gao, F., Zhao, X., Chen, J., Cao, X. (2011). Transient flow analysis in reactor coolant pump systems during flow coastdown period. Nuclear Engineering and Design, vol. 241, no. 2, p. 509-514, D0l:10.1016/j.nucengdes.2010.09.033.
[2]	Poullikkas, A. (2000). Two phase flow performance of nuclear reactor cooling pumps. Progress in Nuclear Energy, vol. 36, no. 2, p. 123-130, DOI:10.1016/S0149-1970(00)00007-X.
[3]	Cho, Y.-J., Kim, Y.-S., Cho, S., Kim, S., Bae, B.-U., Chung, H.-J., Youn, Y.-J., Park, J.-K., Choi, H.-S., Jeon, W.-J. Kim, B.-D., Kwon, T.-S., Song, C.-H. (2014). Advancement of reactor coolant pump
(RCP) performance verification test in KAERI. Proceedings of the 22nd International Conference on Nuclear Engineering, vol. 2B, p. 1-6, D0l:10.1016/S0149-1970(00)00007-X.
[4]	Bing, H., Cao, S.L. (2013). Multi-parameter optimization design, numerical simulation and performance test of mixed-flow pump impeller. Science China Technological Sciences, vol. 56, no. 9, p. 2194-2206, D0l:10.1007/s11431-013-5308-0.
[5]	Kim, S., Lee, K.-Y., Kim, J.-H., Kim, J.-H., Jung, U.-H., Choi, Y.-S. (2015). High performance hydraulic design techniques of mixed-flow pump impeller and diffuser. Journal of Mechanical Science and Technology, vol. 29, no. 1, p. 227-240, D0l:10.1007/s12206-014-1229-5.
[6]	Miyabe, M., Furukawa, A., Maeda, H., Umeki, I., Jittani, Y. (2008). On improvement of characteristic instability and internal flow in mixed flow pumps. Journal of Fluid Science and Technology, vol. 3, no. 6, p. 732-743, D0l:10.1299/jfst.3.732.
[7]	Miyabe, M., Maeda, H., Umeki, I., Jittani, Y., (2006). Unstable head-flow characteristic generation mechanism of a low specific speed mixed flow pump. Journal of Thermal Science, vol. 15, no. 2, p. 115-120, D0l:10.1007/s11630-006-0115-6.
[8]	Cheong, J.-S. (2000). An analytical prediction on the pump-induced pressure pulsation in a pressurized water reactor. Annals of Nuclear Energy, vol. 27, no. 15, p. 1373-1383, D0l:10.1016/S0306-4549(99)00132-2.
[9]	Wang, C., Yi, T.,Wu, Z., Liu, H., Liang, J. (2009). Analysis on pressure fluctuations of unsteady flow field in mixed-flow main coolant pump. Journal of Power Engineering, vol. 29, no. 11, p. 1036-1040.
[10]	Bing, H., Cao, S., Tan, L., Zhu, B. (2013). Effects of meridional flow passage shape on hydraulic performance of mixed-flow pump impellers. Chinese Journal of Mechanical Engineering, vol. 26, no. 3, p. 469-475, D0l:10.3901/cjme.2013.03.469.
[11]	Chan, A.M.C., Kawaji, M., Nakamura, H., Kukita, Y. (1999). Experimental study of two-phase pump performance using a full size nuclear reactor pump. Nuclear Engineering and Design, vol. 193, no. 1-2, p. 159-172, D0I:10.1016/S0029-5493(99)00150-8.
[12]	Xie, R., Shen, F., Wang, X., Shan, Y. (2011). A new CFD-based method study on modelling design of nuclear main pump impeller. Proceedings of the Asia-Pacific Power and Energy Engineering Conference, p. 1-4, D0l:10.1109/ appeec.2011.5748945.
[13]	Zhang, N., Yang, M., Gao, B., Li, Z., Ni, D., (2015). Experimental investigation on unsteady pressure pulsation in a centrifugal pump with special slope volute. Journal of Fluids Engineering, vol. 137, no. 6, p. 061103, D0l:10.1115/1.4029574.
[14]	Knierim, C., Baumgarten, S., Fritz, J., Coon, M.T. (2005). Design process for an advanced reactor coolant pump for a 1400 MW nuclear power plant. ASME Fluids Engineering Division Summer Meeting and Exhibition, p. 1-7, D0I:10.1115/ FEDSM2005-77340.
[15]	Kato, C., Mukai, H., Manabe, A. (2003). Large-Eddy simulation of unsteady flow in a mixed-flow pump. International Journal of Rotating Machinery, vol. 9, no. 5, p. 345-351, D0I:10.1155/ S1023621X03000320.
[16]	Kato, C., Mukai, H., Manabe, A. (2002). LES of internal flows in a mixed-flow pump with performance instability. Proceedings
of ASME Joint U.S.-European Fluids Engineering Division Conference, vol. 2, p. 1-8, D0l:10.1115/FEDSM2002-31205.
[17]	Posa, A., Lippolis, A., Verzicco, R., Balaras, E. (2011). Large-eddy simulations in mixed-flow pumps using an immersed-boundary method. Computers & Fluids, vol. 47, no. 1, p. 3343, DOI:10.1016/j.compfluid.2011.02.004.
[18]	Inoue, Y., Nagahara, T. (2012). Application of PIV for the flow field measurement in a mixed-flow pump. IOP Conference Series: Earth and Environmental Science, vol. 15, no. 2, p. 022011, DOI:10.1088/1755-1315/15/2/022011.
[19]	Fernandez, J., Blanco, E., Santolaria, C., Scanlon, T.J., Stickland, M.T. (2002). A numerical analysis of a mixed flow pump. ASME Joint U.S.-European Fluids Engineering Division Summer Conference, p. p. 791-798, DOI:10.1115/ FEDSM2002-31185.
[20]	Goto, A., Takemura, T., Zangeneh, M., (1996). Suppression of secondary flows in a mixed-flow pump impeller by application of three-dimensional inverse design method: part 1-design and numerical validation. Journal of Turbomachinery, vol. 118, no. 3, p. 536-543, DOI:10.1115/1.2836701.
[21]	Srivastava, S., Roy, A.K., Kumar, K. (2014). Design of a mixed flow pump impeller blade and its validation using stress analysis. Procedia Materials Science, vol 6, p. 417-424, DOI:10.1016/j.mspro.2014.07.053.
[22]	Takayama, Y., Watanabe, H. (2009). Multi-objective design optimization of a mixed-flow pump. Proceedings of the ASME Fluids Engineering Division Summer Meeting, vol. 1, p. 1-9, DOI:10.1115/FEDSM2009-78348.
[23]	Voke, P.R. (1996). Subgrid-scale modelling at low mesh Reynolds number. Theoretical and Computational Fluid Dynamics, vol. 8, no. 2, p. 131-143, DOI:10.1007/ BF00312367.
[24]	Yang, Z., Wang, F., Zhou, P. (2012). Evaluation of subgrid-scale models in large-eddy simulations of turbulent flow in a centrifugal pump impeller. Chinese Journal of Mechanical Engineering, vol. 25, no. 5, p. 911-918, DOI:10.3901/ CJME.2012.05.911.
[25]	Yamanishi, N., Fukao, S., Qiao, X., Kato, C., Tsujimoto, Y., (2006). LES simulation of backflow vortex structure at the
inlet of an inducer. Journal of Fluids Engineering, vol. 129, no. 5, p. 587-594, DOI:10.1115/1.2717613.
[26]	Piomelli, U. (1999). Large-eddy simulation: achievements and challenges. Progress in Aerospace Sciences, vol. 35, no. 4, p. 335-362, DOI:10.1016/S0376-0421(98)00014-1.
[27]	Eastwood, S.J., Tucker, P.G., Xia, H., Klostermeier, C. (2009). Developing large eddy simulation for turbomachinery applications. Philosophical transactions. Series A, Mathematical, Physical, and Engineering Sciences, vol. 367, no. 1899, p. 2999-3013, DOI:10.1098/rsta.2008.0281.
[28]	Ducros, F., Ferrand, V., Nicoud, F., Weber, C., Darracq, D., Gacherieu, C., Poinsot, T., (1999). Large-Eddy simulation of the shock/turbulence interaction. Journal of Computational Physics, vol. 152, no. 2, p. 517-549, DOI:10.1006/ jcph.1999.6238.
[29]	Kato, C., Kaiho, M., Manabe, A. (2003). An overset finite-element large-eddy simulation method with applications to turbomachinery and aeroacoustics. Journal of Applied Mechanics, vol. 70, no. 1, p. 32-43, DOI:10.1115/1.1530637.
[30]	Posa, A., Lippolis, A., Balaras, E. (2015). Large-Eddy simulation of a mixed-flow pump at off-design conditions. Journal of Fluids Engineering, vol. 137, no. 10, p. 101302, DOI:10.1115/1.4030489.
[31]	Byskov, R.K., Jacobsen, C.B., Pedersen, N., (2003). Flow in a Centrifugal Pump Impeller at Design and Off-Design Conditions-Part II: Large Eddy Simulations. Journal of Fluids Engineering, vol. 125, no. 1, p. 73-83, DOI:10.1115/1.1524586.
[32]	Zhang, N., Yang, M., Gao, B., Li, Z., Ni, D. (2016). Investigation of rotor-stator interaction and flow unsteadiness in a low specific speed centrifugal pump. Strojniški vestnik - Journal of Mechanical Engineering, vol. 62, no. 1, p. 21-31, DOI:10.5545/sv-jme.2015.2859.
[33]	Sang-hyuk, Lee, S.H., Wang, Y.-Q., Song, J.-L., (2010). Fourier and wavelet transformations application to fault detection of induction motor with stator current. Journal of Central South University of Technology, vol. 17, no. 1, p. 93-101, DOI:10.1007/s11771-010-0016-4.
Strojniški vestnik - Journal of Mechanical Engineering 62(2016)4, 243-251 © 2016 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2015.3141	Original Scientific Paper
Determination of Coulomb's Friction Coefficient Directly from Cylinder Compression Tests
Deniz Duran* - Celalettin Karadogan Atilim University Metal Forming Centre of Excellence, Turkey
In this paper, a new method is proposed for the determination of Coulomb's friction coefficient directly from cylinder compression tests. It is based on measuring the immigrated contact area (ICA), which is defined as the lateral surface that comes into contact with the platens after deformation. Preliminary sensitivity analyses showed that ICA is only a function of friction and the strain-hardening exponent at room temperature when a power law relation between the true stress and the true plastic strain is considered. Through intensive numerical simulations by using a code-driven simulation environment, an inverse calculation is done which makes the determination of the friction coefficient possible by using ICA and the strain-hardening exponent of the investigated material. At the end of compression, ICA is usually clearly visible without any precautions; thus, a simplified script is supplied which calculates ICA through digital image analysis made on the end faces of the compressed specimens. This paper includes the complete procedure to determine Coulomb's friction coefficient and a script in which the proposed method is embedded entirely. In addition, practical case studies demonstrated the ease of applicability of the proposed method by employing different materials.
Keywords: Coulomb's friction coefficient, cylinder compression test, immigrated contact area, digital image analysis Highlights
•	Immigrated contact area (ICA) is clearly visible after a cylinder compression test.
•	ICA is found to depend only on the friction and the strain-hardening exponent.
•	A directly applicable procedure is developed to evaluate friction coefficients.
•	Associated digital image analysis-based source code is provided.
0 INTRODUCTION
One of the most important factors for metal forming simulations is the selection of correct friction coefficients for the tool-workpiece interface since it considerably affects the material flow. Friction coefficients may depend on various parameters, such as contact pressure, material flow stress, relative sliding velocity, and other factors [1]. It is also well known that the friction coefficient may change during the process [2]; however, evaluation of this complicated phenomenon is subject to devoted tribological tests [3] and [4]. This challenge is usually overcome by a simple friction model with single and easy-to-determine friction coefficients.
Practically, there are three methodologies to determine friction coefficients. The first one is based on a direct computation of the friction coefficient from the experimental loads. The pin-on-disc test is a very conventional example of this category. Although this methodology is easy and elegant, the created physical conditions do not satisfactorily represent the tribology of metal-forming processes. Furthermore, the issue of virtual correctness is ignored in this approach. For a virtually-correct simulation, the numerical value of the virtual friction coefficient (or parameters) should yield a simulation result comparable to the conducted
experiment. This issue is clearly shown in the studies of Tan [5], and Hatzenbichler et al. [6]. The second methodology is an inverse approach, in which the simulation is replaced by an analytical computation, such as an upper-bound analysis, as demonstrated in [7] and [8]. Most of the practical problems are, however, difficult to solve with analytical methods due to their complicated nature. Therefore, this second methodology again lacks virtual correctness because the obtained friction coefficient is to be used in a numerical computation instead. The third methodology makes use of inverse approaches, in which various numerical simulations corresponding to a defined experiment are conducted to determine the friction coefficient, usually by comparing either the geometry of the part or the process forces. For example, ring compression tests and double-cup extrusion tests are designed to assess the geometry of the parts [9] and [10]; whereas in process tests, such as forward rod extrusion or backward can extrusion, the force-displacement history is used for the evaluation [11] and [12]. The use of strain distribution for comparison is also not uncommon [13] and [14]. This methodology, in general, requires the flow curves of the material covering the investigated process range to be known prior to the friction evaluation [2], [15] and [16]. For instance, the flow curve obtained from
tensile tests is usually not sufficient in terms of strain when analysing bulk metal-forming processes. Therefore, compression tests are mostly preferred. On the other hand, compression tests suffer from an inhomogeneous distribution of field variables due to friction in the platen-specimen interface which is typically disregarded or evaluated simultaneously with the flow curve by finite element-based inverse analysis methods [17].
Laborious inverse numerical analyses are, therefore, necessary to obtain virtually correct friction coefficients. This can be overcome if a forging-like friction experiment in which friction-sensitive behaviour is independent of material behaviour exists. Such an approach allows the determination of the functional relations between the experimental results and the friction coefficients through inverse numerical analyses only once. These functional relations later enable a quick and direct use, still being based on intensive inverse analyses. In cylinder compression tests, the term immigrated contact area (ICA) is introduced by Tan [18] and defined as the lateral surface portion that comes into contact with the platens as deformation proceeds. The preliminary sensitivity analyses show that, at room temperature, ICA is practically independent of material flow curve to a certain extent. Considering the power law formula cf= Ksn as the material flow curve, ICA is independent of strength coefficient K and only depends on the strain-hardening exponent n. Based on this material independence, this study aims to build an environment where one can directly determine the friction coefficient in compression tests using ICA and n as the inputs of a set of equations derived from intensive numerical inverse analyses and f as the ultimate output.
Accurate measurement of the ICA is considered to be a crucial matter and, therefore, should be addressed. In the literature, many successful implementations of digital image correlation (DIC) techniques can be found. For instance, Skozrit et al. [19] investigated the damage behaviour of aluminium alloys under different load cases and at different strain rates. In order to precisely calibrate the parameters of the constitutive models, the displacement and temperature distribution of the specimens have been measured by DIC and infrared thermography (IT), respectively. Min et al. [20] developed a procedure to accurately obtain principal stresses and strains, radius and sheet thickness at the pole of bulge specimens by means of DIC. Another DIC-based method has been introduced in [21] to detect forming limit strains correctly in case of multiple local necks occurrence
while testing sheet metal formability. Knysh and Korkolis [22] measured the inelastic heat fraction during the plastic deformation of different metals at different strain rates. Local kinematics have been measured using a non-contact DIC system; whereas temperature gradients have been measured with an infrared camera. DIC techniques take part in many more applications in various fields [23] to [29] and thus became a trusted method for experimental measurements. For this reason, the authors decided to use DIC to appropriately and objectively measure the ICA.
The rest of the paper is structured as follows: In the second section on Material and Method, the numerical investigation of the cylinder compression test for a specified broad range of materials and process parameters are introduced. Then, this numerical investigation is used to analyse the dependency of ICA on the material parameters and the friction coefficient. Finally, a complete procedure for the evaluation of friction coefficients based on the cylinder compression test results is developed. The third section, on Experimental Study, provides cases in which the developed procedure is utilized, and ICA is measured by DIC. The necessary tools for the application of the procedure are explained, and a complete MATLAB script is provided in the Appendix. The paper concludes in the fourth section on Discussions and Conclusion, where a critique of the developed methodology is supplied.
1 MATERIAL AND METHOD
A numerical study based on an inverse analysis is performed in order to investigate the relationship between the ICA, material parameters, and friction coefficient. Covering a practical range of material parameters and friction coefficients relying on literature, various cylinder compression test simulations are executed and evaluated in an automated manner. Based on the simulation results, the relationship between the ICA, friction coefficient, and strain-hardening exponent n, is formulated. Utilizing the developed formulas, a procedure for the evaluation of the friction coefficient is provided.
1.1 Finite Element Model
Based on the suggested specimen dimensions in [30], a cylinder of 10 mm diameter and 15 mm height is considered as the compression test specimen. Cylindrical geometry and the loading condition are both rotationally and axially symmetrical. Therefore,
a quarter section of the cylinder is modelled and discretized by four-noded full integration quadrilateral elements. Based on the conducted convergence studies for the finite element model, a special mesh structure is used, with which a precise calculation of ICA is sought, see Fig. 1. The mesh is gradually and consecutively refined four times in the region where immigration of the lateral surface to the end face takes place. This refinement corresponds to an element size decrease from 0.25 mm to 0.0156 mm. With the used mesh, simulations run free of problems up to 50 % compression, i.e. 7.5 mm total axial displacement.
Following the completion of simulations, the width of ICA, which is illustrated with u in Fig. 2, is calculated based on a simple coordinate-checking code executing at the exterior nodes of the final deformed mesh.
Fig. 1. Finite element mesh used in the calculations 1.2 Parameter Ranges
The range of parameters and the modelling approach for the material behaviour is significant with regards to the validity of the numerical study. Most of the metals exhibit strain-hardening behaviour around room temperature which is well reflected by a simple power law expression as the flow curve (f Ken). During cylinder compression around room temperature, moderate contact pressures occur. In such case, a linear relation between frictional forces and normal forces as set forth by Coulomb's friction model, Tf=, achieves non-objectionable results. Based on the aforementioned observations, the parameters to be investigated are chosen as strength coefficient K, strain-hardening exponent n, and friction coefficient f. Available models and values in the literature, such as [30], are considered for a comprehensive range of these parameters. The selected ranges and specific values for K, n and f are as follows:
K = {250, 500, 750, 1250} [ MPa ],
n = {0.05, 0.125, 0.25, 0.375, 0.5} [-],
V = {0.05, 0.08, 0.1, 0.12, 0.15, 0.2, 0.25, 0.3} [ -] .
A full factorial combination of this parameter sets results in a total of 160 simulations to be performed. These simulations are executed and evaluated automatically by a driving MATLAB script, which produces the input file, submits the simulation to MSC-Marc, extracts the final geometry, calculates u and stores the result in accordance with the used parameters.
The possible influence of elastic modulus on ICA is investigated by considering three different values: 70 GPa, 210 GPa for elastoplasticity and an infinite value for rigid plastic behaviour. The simulation results revealed that the influence of elastic modulus on ICA is minimal, less than 3 % in the worst case and that it does not show a significant trend. As a reference, 210 GPa is taken; however, the proposed method can be safely used for any type of metal regardless of its elastic modulus.
axial symmetry
platen movement
rotational symmetry
Fig. 2. Finite element mesh after 50 % compression and calculation of u
1.3 Mathematical Basis
For each combination of K, n and f ; u values are calculated through simulations. As explained in the previous sections, u is the width of the lateral (initially cylindrical) surface portion that comes into contact with the platens as deformation proceeds. Since a relation for f with regards to other parameters is sought, u values are considered as input, whereas f values are output. Fig. 3 shows f values with respect to u where simulations are grouped in terms of separate n values. A glance at these figures reveals that
u
any change in the value of K causes only negligible variation in u; therefore, u is practically independent of K and the dependency between u and / changes with respect to n. It is the general practice to use a single representative flow curve for cold forming analysis. Therefore, the cylinder compression test and the evaluation of n have to be carried out under the same conditions as the tests to determine the representative flow curve.
Bearing the strain-hardening exponent n of the investigated material in mind, one can use Fig. 3 to evaluate / by employing the measured width u. However, for the sake of convenience, interpolation equations are obtained by using a best-fit cubic polynomial equation in the following form:
№ = fn (u ) = anu3 + bnu2 + cnu + dn
(2)
where the coefficients of fn are only a function of n. Eq. (2), a third-degree polynomial equation that is chosen for curve fitting and has no physical meaning. Taking into account the goodness of fits as shown in Fig. 3, Eq. (2) represents well the considered comprehensive physical parameter range for all possible combinations. The computed coefficients, namely an, bn, cn and dn are given in Fig. 3 and summarized in Table 1.
Table 1. Coefficients of Eq. (2) computed for each n value
n	an	bn	Cn	dn
0.05	0.42708	-0.61813	0.37612	0.00752
0.125	0.2229	-0.2716	0.2518	0.0249
0.25	0.11747	-0.11819	0.23088	0.03609
0.375	0.08748	-0.1016	0.26468	0.04036
0.5	0.10083	-0.15695	0.3218	0.04166
These dependencies can be interpolated further using cubic polynomial functions, see also Fig. 4:
an = -9.78142« +11.32378«2 - 4.23226« + 0.60661, bn = 18.01546« - 20.67774« + 7.38683« - 0.92918, cn = -8.26771«3 + 9.20005« - 2.88130« + 0.49425, dn = 0.78652« - 0.91056« + 0.35799« + 0.00788. (3)
By means of the derived equation sets, friction coefficient for a platen-workpiece material pair can be computed by following the steps given below:
•	Perform a tension test to obtain the strength coefficient K and strain-hardening exponent n. If K and n fall within the range given in Eq. (1);
•	perform a compression test up to 50 % height reduction by using cylinders of 10 mm diameter and 15 mm height;
a)
b)
у = 0.08748X3- 0.10160x2 + 0.26468X + 0.04036 - Rz = 0.9
di
0.50 0.75 и [mm]
у = 0.10083X3- 0.15695хг + 0.32180x + 0.04166 - R! = 0.99961
ei
0.00 -
0.00
0.50 0.75 и [mm]
Fig. 3. Simulation results showing K-independence and u-dependence of /; for all K and / at a) n = 0.05, b) n = 0.125, c) n = 0.25, d) n = 0.375, e) n = 0.5
c)
•	measure u on one of the end faces of the compressed cylinders (elaborated further in Section 3);
•	using Eq. (3), calculate the coefficients an, bn, cn and dn;
•	using Eq. (2), calculate the friction coefficient. Neither tensile nor compression test speeds are
specified in the test procedure. Since the application takes place at room temperature, the effect of the strain rate is ignored. This approach is reasonable since the flow behaviour of most of the metals shows only a little sensitivity to strain rate at room temperature. This fact supports why a single flow curve is chosen to represent the material behaviour. However, it should be emphasized that relatively low test speeds should be preferred in order to avoid deformation heating of the specimens.
The numerical simulations and computations above are nothing more than an inverse numerical analysis that has to be performed when the inverse approach is to be used to determine the friction coefficient. As distinct from the conventional inverse analysis, results of the proposed approach can be used directly for all materials within the considered range (Eq. (1)). Noting that the proposed procedure substitutes only the numerical simulations of the inverse approach, the proof of its practicality is now complete.
0.0 0.1 0.2 0.3 0.4 0.5 0.6
n [-]
Fig. 4. Dependency of cubic function coefficients on the strain-hardening exponentn
2 EXPERIMENTAL STUDY
Following the proof and the construction of the form of relation between the strain-hardening exponent n, ICA, and the friction coefficient; a practical application of the proposed method will
be demonstrated in this section. Several examples using various material types are demonstrated, and a digital image analysis technique is introduced for the measurement of the width u of ICA on an image taken from a compressed specimen. This tool is based on the digital image analysis technique and ensures more accurate measurements as compared with available direct measurement techniques from the specimens. The complete MATLAB script corresponding to this digital image analysis technique is given in the Appendix.
2.1 Digital Image Analysis
After the compression tests, in general, ICA can be observed effortlessly by a distinct contrast difference. To improve the accuracy of the measurement, a digital image analysis technique is preferred based on the image of the compressed face. Shot from a perpendicular direction, a digital camera image of any of the end faces mostly yields a clear view of ICA. Fig. 5 shows examples of ICAs of steel (a), aluminum (b), and copper (c) specimens. The top image in Fig. 5 shows the result of the evaluation method demonstrating the fitted spline with an enlarged view. The bottom three images in Fig. 5 demonstrate the distinct contrast difference that is obtained from various materials.
The inner and outer circles generated on the end face of the compressed cylinders are used to compute ICA. However, due to material and contact interface inhomogeneities, these circles are not perfect. Therefore, equivalent radii should be calculated corresponding to these areas. For this purpose, spline curves are fitted to the exterior of the regions and the respective equivalent circles' radii are determined. The difference between these two equivalent radii yields the equivalent u. The only necessity to compute the areas in terms of real distances from pixel distances is to use a calibration factor that converts pixel distance to real distance in the plane where the observed end face lies. This can be achieved either by using a pre-specified distance between two points located on the specimen face, or by means of a precision ruler lying in the same plane as the observed contact face. The proposed quasi-automated digital image analysis approach is simple and reliable and provides values within at most 3.5 % mean deviation for around 20 evaluations performed by different users. The accuracy of the measured u and the corresponding friction coefficient is, in fact, directly related to the camera resolution. If the specimen face is covered only by 2 megapixels, a 3-pixel error in positioning
a marker corresponds to 0.03 mm variation in u and a variation of 0.0063 in the value of the friction coefficient regarding Fig. 3. This error can be reduced further by using a camera with higher resolution.
The complete MATLAB script for the quasi-automated computation of u is given in the Appendix. This script imports the picture of the compressed specimen and uses the millimetric calibration distance, which will be used later to compute the calibration factor giving the ratio of the pixel distance to the real distance. Pixels corresponding to the end points of the calibration distance are located. Following that, a sufficient number of markers are placed on the exterior of the inner circle (i.e. the area corresponding to the end face of the initial cylinder) and on the exterior of the outer circle (i.e. the final total contact area), respectively. The cases shown in Fig. 5 can be considered as good examples of a sufficient number of markers. The coordinates of the placed markers are used periodically to construct the closed cubic splines which are then used to compute the areas within each of these closed curves as shown in Fig. 5. In case deviations are observed between the interpolating spline and the borders of the regions, the procedure should be repeated for better accuracy. Finally, the circumscribed areas are used to compute the equivalent radii, the difference of which yields the equivalent u.
2.2 Case Studies
To determine the friction coefficient between a die-workpiece material pair using the proposed method, two platens out of the die material with the same surface characteristics (roughness, finish, and lubrication) have to be used. Furthermore, cylindrical specimens of 10 mm diameter and 15 mm height are needed. Compression down to half-height of 7.5 mm generates the necessary contact surface, namely ICA, for the evaluation of friction coefficient.
To demonstrate the applicability of the proposed method, cylindrical specimens are produced out of three different materials: As drawn low-alloy steel (20MnCr5), 6000 series aluminium (AlMgSi1) and cold-worked oxygen-free copper. Tension tests at room temperature are performed to extract the flow curve of the materials and are approximated by the power law function f Ksn as provided in Table 2. The used platen material is a variant of hot-work tool steel (X40CrMoV5-1), and surfaces are lightly exposed to PTFE-based spray lubricant. The authors would like to express that it is misleading to consider Table 2 to be a list of friction coefficients for readers'
access. These friction coefficients depend on surface roughness values and the lubrication condition for these specific cases.
Table 2. Material properties and computed friction coefficients
Material	K [MPa]	n [-]	u [mm]	V [-]
20MnCr5 (As drawn)	867	0.135	0.74	0.15
AlMgSil (Cold worked)	537	0.21	1.1	0.29
Copper (99.9 %)	410	0.106	0.89	0.18
Based on the processed images in Fig. 5, the equivalent u values and respective friction coefficients are evaluated for these three materials using the complete MATLAB script given in the Appendix. These values are given in Table 2 together with the materials' properties. The common and straightforward practice is to conduct a lengthy inverse analysis of the performed compression test to evaluate the friction coefficient. The proposed approach, however, bears all these inverse analyses, such as a meta-model, and allows the evaluation of the same friction coefficient in a direct manner, effortlessly.
3 CONCLUSION
In this work, a complete procedure for the direct determination of Coulomb's friction coefficients from cylinder compression tests is demonstrated. The study utilizes cylinder compression test as the friction experiment; immigrated contact area (ICA) as the indirect measure of the friction coefficient; power law cf= Ksn as the material flow curve at room temperature; Coulomb's formulation Tf=as the friction model; bilinear-displacement implementation as the numerical algorithm and MSC Marc Mentat as the simulation environment.
Numerical simulations are used to prove the dependence of ICA only on the friction coefficient and the strain-hardening exponent of the material for a broad range of material parameters. Equations are obtained to compute friction coefficients directly from the equivalent width u of ICA and the strain-hardening exponent n. Finally, an accurate measurement method based on digital image analysis is developed for the measurement of u from the compressed specimens. The MATLAB script for the complete measurement method is developed and provided in the Appendix. The usability of the method is demonstrated through compression tests performed on various materials.
The proposed method is purely based on the broadly utilized inverse approach benefiting from all of its aspects. Its practice, however, is totally
direct and straightforward and eliminates the burden of performing the inverse numerical analysis for the broad material range Eq. (1). The maximum amount of overall error coming from Eqs. (2) and (3) can be seen in Figs. 3 and 4. For the given range Eq. (1), the maximum error made using the interpolation functions Eqs. (2) and (3) in reproducing the provided base data is 10 % where the overall average error is only 2 %. Considering Fig. 3, one can conclude that if the range for the friction coefficient is confined between 0.1 < ju<0.25, the errors coming from the interpolation functions Eqs. (2) and (3) decline further. The maximum possible error obtained in this case is only 5 %, and the average error drops to 1.5 % in reproducing the provided base data Eq. (1).
A major benefit of the presented method is that it is highly practical, and hence does not require any special testing equipment for friction characterization. The compression tests for flow curve evaluation can also be used for friction characterization. Employing the proposed method, the friction coefficient becomes readily available which is necessary for the inverse evaluation of the compression test for correction purposes. The necessity of performing a tensile test for the determination of the strain-hardening exponent
can also be eliminated by using the early stages of the compression test data instead.
The applicability of the method is demonstrated for various materials using ordinary photography techniques. The results of the experimental study indicated that the proposed method replaces the lengthy inverse analysis with acceptable accuracy.
4 REFERENCES
[1]	Bowden, F.P., Tabor, D. (1982). Friction - An Introduction to Tribology (Reprint Edition), Robert E. Krieger Publishing Company, Florida.
[2]	Petersen, S.B., Martins, P.A.F., Bay, N. (1997). Friction in bulk metal forming: A general friction model vs. the law of constant friction. Journal of Materials Processing Technology, vol. 66, no. 1-3, p. 186-194, D0I:10.1016/S0924-0136(96)02518-6.
[3]	Groche, P., Müller, C., Stahlmann, J., Zang, S. (2013). Mechanical conditions in bulk metal forming tribometers-Part one. Tribology International, vol. 62, p. 223-231, DOI:10.1016/j.triboint.2012.12.008.
[4]	Groche, P., Stahlmann, J., Müller, C. (2013). Mechanical conditions in bulk metal forming tribometers-Part two. Tribology International, vol. 66, p. 345-351, DOI:10.1016/j. triboint.2012.11.028.
[5]	Tan, X. (2002). Comparisons of friction models in bulk metal forming. Tribology International, vol. 35, no. 6, p. 385-393, DOI:10.1016/S0301-679X(02)00020-8.
[6]	Hatzenbichler, T., Harrer, O., Wallner, S., Planitzer, F., Kuss, M., Pschera, R., Buchmayr, B. (2012). Deviation of the results obtained from different commercial finite element solvers due to friction formulation. Tribology International, vol. 49, p. 7579, DOI:10.1016/j.triboint.2011.12.020.
[7]	Ebrahimi, R., Najafizadeh, A. (2004). A new method for evaluation of friction in bulk metal forming. Journal of Materials Processing Technology, vol. 152, no. 2, p. 136-143, DOI:10.1016/j.jmatprotec.2004.03.029.
[8]	Martin, F., Sevilla, L., Camacho, A., Sebastian, M.A. (2013). Upper bound solutions of ring compression test. Procedia Engineering, vol. 63, p. 413-420, DOI:10.1016/j. proeng.2013.08.211.
[9]	Robinson, T., Ou, H., Armstrong, C.G. (2004). Study on ring compression test using physical modelling and FE simulation. Journal of Materials Processing Technology, vol. 153-154, p. 54-59, DOI:10.1016/j.jmatprotec.2004.04.045.
[10]	Schrader, T., Shirgaokar, M., Altan, T. (2007). A critical evaluation of the double cup extrusion test for selection of cold forging lubricants. Journal of Materials Processing Technology, vol. 189, no. 1-3, p. 36-44, DOI:10.1016/j. jmatprotec.2006.11.229.
[11]	Hsu, T.C., Huang, C.C. (2003). The friction modeling of different tribological interfaces in extrusion process. Journal of Materials Processing Technology, vol. 140, no. 1-3, p. 4953, DOI:10.1016/S0924-0136(03)00724-6.
[12]	Kuzman, K., Pfeifer, E., Bay, N. Hunding, J. (1996). Control of material flow in a combined backward can-forward rod extrusion. Journal of Materials Processing Technology, vol. 60, no. 1-4, p. 141-147, DOI:10.1016/0924-0136(96)02319-9.
[13]	Giuliano, G. (2015). Evaluation of the Coulomb Friction Coefficient in DC05 Sheet Metal Forming. Strojniški vestnik
-	Journal of Mechanical Engineering, vol. 61, no. 12, p. 709713, DOI:10.5545/sv-jme.2015.2733.
[14]	Hatipoglu, H.A., Karadogan, C. (2014). A Methodology to determine friction coefficient in stretch forming of sheet aircraft panels. Conference Proceedings of International Deep-Drawing Research Group.
[15]	Hu, C.L., Ou, H.A., Zhao, Z. (2015). An alternative evaluation method for friction condition in cold forging by ring with boss compression test. Journal of Materials Processing Technology, vol. 224, p. 18-25, DOI:10.1016/j.jmatprotec.2015.04.010.
[16]	Bay, N., Gerved, G. (1987). Tool/workpiece interface stresses in simple upsetting. Journal of Mechanical Working Technology vol. 14, no. 3, p. 263-282, DOI:10.1016/0378-3804(87)90013-1.
[17]	Cho, H., Ngaile, G., Altan, T. (2003). Simultaneous determination of flow stress and interface friction by finite element based inverse analysis technique. CIRP Annals
-	Manufacturing Technology, vol. 52, no. 1, p. 221-224, DOI:10.1016/S0007-8506(07)60570-8.
[18]	Tan, X. (2001). Friction reducing contact area expansion in upsetting. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 215, no. 2, p. 189-200, DOI:10.1243/1350650011541828.
[19]	Skozrit, I., Frančeski, J., Tonković, Z., Surjak, M., Krstulović-Opara, L., Vesenjak, M., Kodvanj, J., Gunjević, B., Lončarić, D. (2015). Validation of numerical model by means of digital image correlation and thermography. Procedia Engineering, vol. 101, p. 450-458, DOI:10.1016/j.proeng.2015.02.054.
[20]	Min, J., Stoughton, T.B., Carsley, J.E., Carlson, B.E, Lin, J., Gao, X. (2016). Accurate characterization of biaxial stress-strain response of sheet metal from bulge testing. International Journal of Plasticity, in Press, DOI:10.1016/j. ijplas.2016.02.005.
[21]	Vysochinskiy, D., Coudert, T., Hopperstad, O.S., Lademo, O.G., Reyes, A. (2016). Experimental detection of forming limit strains on samples with multiple local necks. Journal of Materials Processing Technology, vol. 227, p. 216-226, DOI:10.1016/j.jmatprotec.2015.08.019.
[22]	Knysh, P., Korkolis, Y.P. (2015). Determination of the fraction of plastic work converted into heat in metals. Mechanics of Materials, vol. 86, 71-80, DOI:10.1016/j. mechmat.2015.03.006.
[23]	Sutton, M.A., Schreier, H.W, Orteu, J-J.(2009). Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications, Springer, New York.
[24]	Milosevic, M., Milosevic, N., Sedmak, S., Tatic, U., Mitrovic, N., Hloch, S., Jovicic, R. (2016). Digital image correlation in analysis of stiffness in local zones of welded joints. Technical Gazette - Tehnički vjesnik, vol. 23, no. 1, p. 19-24, DOI:10.17559/TV-20140123151546.
[25]	Majnaric, I., Hladnik, A., Muck, T., Bolanca Mirkovic, I. (2015). The influence of ink concentration and layer thickness on yellow colour reproduction in liquid electrophotography toner. Technical Gazette - Tehnički vjesnik, vol. 22, no. 1, p. 145149, DOI:10.17559/TV-20140321230455.
[26]	Foldyna J., Zelenak M., Klich J., Hlavacek P., Sitek L., Riha Z. (2015). The measurement of the velocity of abrasive particles at the suction part of cutting head. Technical Gazette -Tehnicki vjesnik, vol. 22, no. 6, p. 1441-1446, DOI:10.17559/ TV-20140214152925.
[27]	Klancnik, S., Ficko, M., Balic, J., Pahole, I. (2015). Computer Vision-Based Approach to End Mill Tool Monitoring. International Journal of Simulation Modelling, vol. 14, no. 4, p. 571-583, DOI:10.2507/IJSIMM14(4)1.301.
[28]	Cajal, C., Santolaria, J., Samper, D., Garrido, A. (2015). Simulation of Laser Triangulation Sensors Scanning for Design and Evaluation Purposes. International Journal of Simulation Modelling, vol. 14, no. 2, p. 250-264, DOI:10.2507/ IJSIMM14(2)6.296.
[29]	Chrysochoos, A., Huon, V., Jourdan, F., Muracciole, J.M., Peyroux, R., Wattrisse, B. (2010). Use of full-field digital image correlation and infrared thermography measurements for the thermomechanical analysis of material behaviour, Strain, vol. 46, no. 1, p. 117-130, DOI:10.1111/j.1475-1305.2009.00635.x.
[30]	Banabic, D., Bunge, H.-J., Pöhlandt, K., Tekkaya, A.E. (2000). Formability of Metallic Materials, Springer-Verlag Berlin/ Heidelberg/New York/Tokyo, DOI:10.1007/978-3-662-04013-3.
5 APPENDIX
MATLAB script for the evaluation of the friction coefficient:
SpecimenImage = imread('CompSpecimen(1).jpg'); CalibrationDistance = 13.0; % 13 mm SHE = 0.21; % Strain Hardening Exponent
imshow(SpecimenImage); hold on
disp ('Click on the two locations for the length calibration:'); [xCp,yCp,button] = ginput(2);
pixelDistance = sqrt((xCp(2)-xCp(1))A2+(yCp(2)-yCp(1))A2); CalibrationFactor = pixelDistance/CalibrationDistance; plot(xCp,yCp,'y-'); hold on
disp('Pick up sufficient number of ordered points on the inner circle')
disp ('Left mouse button picks points.') %LMB
disp ('Right mouse button picks last point.') %RMB
xi = []; yi = []; n = 0; % Initially empty list
button = 1; % button 1 is the left mouse button
while button == 1
[x,y,button] = ginput(1); plot(x,y,'go') n = n+1; xi(n,1) = x;	yi(n,1) = y;
end
n = n+1; xi(n,1) = xi(1,1); yi(n,1) = yi(1,1); % closed curve % Interpolate with two splines and finer spacing. ti = 1:n; % breakpoint parameters are defined by t tis = 1: 0.01: n;
Fxis = csape(ti,xi,'periodic'); xis = fnval(Fxis,tis); Fyis = csape(ti,yi,'periodic'); yis = fnval(Fyis,tis); % Plot the interpolated curve. plot(xis,yis,'g-'); hold on
% Loop, picking up the points for the outer circle
disp ('Pick up sufficient number of ordered points on the outer circle')
disp ('Left mouse button picks points.') %LMB
disp ('Right mouse button picks last point.') %RMB
xo = []; yo = []; n = 0; % Initially, empty point list
button = 1; % button 1 is the left mouse button
while button == 1
[x,y,button] = ginput(1);	plot(x,y,'ro')
n = n+1; xo(n,1) = x;	yo(n,1) = y;
end
n = n+1; xo(n,1) = xo(1,1); yo(n,1) = yo(1,1); % closed curve
% Interpolate with two splines and finer spacing.
to = 1:n;
tos = 1: 0.01: n;
Fxos = csape(to,xo,'periodic'); xos = fnval(Fxos,tos); Fyos = csape(to,yo,'periodic'); yos = fnval(Fyos,tos); % Plot the interpolated curve. plot(xos,yos,'r-'); hold on; % Calculate areas and Immigrated Surface Width areaInsideInnerSpline = polyarea(xis,yis)/CalibrationFactorA2; areaInsideOuterSpline = polyarea(xos,yos)/CalibrationFactorA2; FinalRadiusOfInitialFace = sqrt(areaInsideInnerSpline/pi()); FinalRadiusOfContactFace = sqrt(areaInsideOuterSpline/pi());
% u = Equivalent width of the Immigrated Contact Area u = FinalRadiusOfContactFace - FinalRadiusOfInitialFace; a = -9.78142*SHEA3 + 11.32378*SHEA2 - 4.23226*SHE + 0.60661; b = 18.01546*SHEA3 - 20.67774*SHEA2 + 7.38683*SHE - 0.92918; c = -8.26771*SHEA3 + 9.200050*SHEA2 - 2.88130*SHE + 0.49425; d = 0.78652*SHEA3 - 0.910560*SHEA2 + 0.35799*SHE - 0.00788;
FrictionCoefficient = a*uA3 + b*uA2 + c*u + d
Strojniški vestnik - Journal of Mechanical Engineering 62(2016)4, 252-259 © 2016 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2015.2969	Original Scientific Paper
Effect of Refrigerant Type and Insulation Thickness on Refrigeration Systems of Land and Sea Vehicles
Ahmet Selim Dalkilig1* - Ali Celen1 - Alican Qebi1 - Somchai Wongwises2 1 Yildiz Technical University Faculty of Mechanical Engineering, Turkey 2 King Mongkut's University of Technology Thonburi, Faculty of Engineering, Thailand
Vehicles are used in the transportation industry to carry temperature-sensitive goods. In fact, they are designed to bear away perishable cartage at specific temperatures. In this study, the cold rooms of a frigoship, railroad car and truck were designed in order to store foods. In addition, the effect of refrigerants and insulation thickness of the cold rooms' refrigeration systems were taken into consideration. In the analysis, R12, R22 and their alternatives of R134a and R410A were selected as refrigerants flowing in the cycle. The evaporator capacity, the condenser capacity, the compressor work, the refrigerant flow rate and the coefficient of performance (COP) of the designed refrigeration systems for each vehicle were determined. It was observed that R134a and R410A had a slightly lower COP and required higher compressor work than R12 and R22 for a condensation temperature. Frigoship, truck and railroad car cold rooms working with R12 had COP values of 2.24, 2.63 and 3.17, respectively. Moreover, the proposed refrigerant R134a can be used in applications with its COP values of 2.16, 2.51 and 3.15 for frigoship, truck and railroad cold rooms, respectively. The influence of the insulation thickness of the wall on evaporator capacity was also investigated. It was observed that the cooling load of the frigoship, railroad car and truck cold rooms decreased 69 %, 68 % and 43 %, respectively, with increase of insulation thickness of the wall from 0.03 m to 0.3 m. Keywords: refrigeration, alternative refrigerants, COP, cold room, insulation, condenser, evaporator
Highlights
•	Frigoship, truck and railroad car cold rooms working with different refrigerants were investigated.
•	The influences of refrigerant type and insulation thickness on refrigeration systems were determined.
•	A case study was performed for anyone wants to design and investigate usual cold room issues with the fundamental equations.
•	R134a was selected as suitable refrigerant by considering COP results and environmental effects.
0 INTRODUCTION
Vehicles are commonly used in today's world in order to keep up the quality and increase the shelf life of fresh, perishable and frozen products during transportation. In this refrigeration process, the most important thing is to keep the temperature at the point where metabolic and microbial decompositions are minimized. In other words, designers design for freezing temperatures, which hinder the growth of bacteria. Moreover, designers should pay attention to the details that are necessary during the refrigeration process's design. For instance, which compressor and condenser are correct and convenient for the cold room system should be known. Additionally, efficient cold-storage space, its location with regard to transportation time from the cold room to the market area, product types stored together at the same time, building materials and seasonal temperature changes of the location should be taken into consideration.
Moureh et al. [1] aimed to reduce the temperature differences throughout palletized cargos of perishable products kept in refrigerated vehicles. For this purpose, temperature and ventilation homogeneities were characterized by using air duct systems. Numerical
and experimental studies were conducted in order to investigate airflow patterns and temperature levels within a typical loaded vehicle with and without a duct air supply system. The results revealed the importance of air ducts in decreasing temperature differences throughout the cargo.
Tso et al. [2] investigated the usage of an air curtain instead of plasticone in the refrigerated truck application. The usage of a plastic strip curtain can be dangerous because of its design for workers during loading and unloading of products. Experimental results showed that the air curtain, having available volume flow rate and velocity, can block the heat transfer through the door. Further, the infiltration heat load can be reduced 11 % by using an air curtain instead of a plastic strip one.
Hoang et al. [3] modeled heat transfer in a refrigerated vehicle. In the study, they investigated the heterogeneity of product temperature in different loads by considering generally near the rear and front. As a result of their study, CPU computing time was found to be less than 1 min instead of 60h for the same configuration (Moureh et al. [4]). Thus, input parameters could be evaluated rapidly, such as external temperature, door openings, wall insulation,
252
*Corr. Author's Address: Yildiz Technical University Faculty of Mechanical Engineering, Turkey, dalkilic@yildiz.edu.tr
thermostat setting, heat of respiration etc., on the load temperatures. Beyond this, the presence of air distribution ducts could be considered, and maximal and minimal load temperature could be predicted. Also, quality and microbiological evolution could be coupled with this approach.
Glouannec et al. [5] experimentally and numerically investigated the thermal properties of an insulation wall by considering the impact of the external environment and durability for refrigerated vans. They aimed to reduce the environmental effects of refrigerated vehicles by improving the insulation design and decreasing energy consumption. Also, the dynamic characterization of the thermal transfer insulation walls was examined. As a result of the study, aerogel and multi-foil insulation layers showed better performance during the daytime period.
Ahmed et al. [6] used phase change materials (PCMs) for refrigerated truck trailer insulation application. Heat and moisture were managed throughout the shipment period to maintain the inside of the truck at a relative humidity and constant temperature. As an analysis of this study, incentive results in lowering peak heat transfer rates and total heat flows into the refrigerated truck trailer were ensured as a heat transfer technology. Accordingly, potentially energy could be saved and pollution from diesel-driven refrigeration equipment was reduced.
Tan et al. [7] analyzed the solidification of phase change material (PCM) and the heat transfer of a cold storage system using tubes to meet the design requirement for the cryogenic energy storage of Liquefied Natural Gas (LNG) refrigerated vehicles. They found as a result of this investigation that the main thermal resistance occurred in the gaseous heat transfer fluid (HTF) inside the tube and the ice layer not only increased radial direction but also propagated it in an axial direction. Further, this study is a good reference for anyone who wants to design and optimize the cold energy storage unit of LNG-refrigerated vehicles.
Cleland [8] demonstrated the significance of refrigerated packaging system design for transportation over important distances. In this study, combining both refrigeration and packaging were taken into consideration and packaging designs for refrigerated food was frequently dictated by the use of low temperature and the high relative humidity conditions. Also, the team that played a role during packaging of perishable food had to balance the needs of marketing staff.
Morganti and Gonzalez-Feliu [9] examined deliveries of perishable goods from the urban
distribution center to local food outlets. For this investigation, Parma, Italy, was chosen since it implemented one of the city logistics projects in Europe in order to improve efficiency and reduce the negative effects of urban food distribution within the city. As a result of this study, the delivery schemes for fresh food products adopted by the corporate chain retailers, independent retailers and food services were analyzed. The policy process, the role of the Parma's public authority and the importance of the wholesale produce market as a Food Hub was searched.
Smale et al. [10] studied patterns of airflow using Computational Fluid Dynamics (CFD) in refrigeration systems so that temperature homogeneity could be checked. The conservation of food quality in chilled and frozen food applications was based on controlling of temperature during cooling, storage, transportation and display. The results indicated that the application of airflow modelling techniques enabled ventilation and temperature to be reduced and the efficiency and effectiveness of refrigeration systems to be increased.
Duret et al. [11] focused on temperature and moisture heterogeneity in cold room systems to ensure quality and safety of food products. In their study, a cold room having four apple pallets was observed. A simplified airflow pattern was recommended according to be experimented air velocity field. They tried to clarify heat exchange, water exchanges and airflow mechanisms by measuring velocity of air, temperatures of product and air, convective heat transfer coefficient and weight loss.
Kou et al. [12] wanted to examine temperature variations in refrigerated systems in order to determine their effect on the quality and microbial growth of packaged products. Owing to this research, two commercial open-refrigerated display cases were investigated under different operating conditions using packaged baby spinach products. The results demonstrated that the temperature variation in the cases was dependent on spatial location, thermostat setting, and defrost cycle period and duration of defrost.
In this investigation, the frigoship with 3391 m3 storage volume, the railroad car with 93 m3 storage volume and the truck with 113 m3 storage volume were designed to keep 57,000 kg of frozen fish, 4,080 kg of frozen meat and 23,000 kg of frozen meat, respectively. Main dimensions of the investigated vehicles were given in Table 1. Frozen fish and meat were kept at -18 °C and -20 °C, respectively. Refrigerant flow rate, condenser capacity, compressor work and coefficient of performance (COP) were figured out for the refrigeration cycle with R12,
R22 and their alternative refrigerants of R134a and R410A. In addition, effect of insulation thickness on evaporator capacity was detected.
Table 1. Main dimensions of the investigated vehicles
Specification
Dimension [m]
Length over all
80.284
Frigoship
Railroad car
Truck
Length of waterline	78.492
Rule length	76.137
Moulded breadth	12.500
Depth	61.00
Design draft	48.00
Summer draft	47.47
Total block coefficient	8.50
Height of railroad car	2.15
Length of railroad car	12.78
Width	2.6
Length overall	14.320
Height overall	3.918
Door size	2.25x2.59
Length	13
Height	3.20
Width	2.70
Door size
1.35x3.20
1 CALCULATION PROCEDURE 1.1 Heat Loss Calculations
The cooling load of a cold room is determined by summing the heat gained from walls, product, air alteration and heat sources, and it is expressed in Eq. (1). The unknown heat gains were taken into consideration as shown in Eq. (2). A 10 % increase was added to the total value of the system's cooling load, which was equal to the evaporator's heat rejection rate from the cold room.
Qt = Qw + Qp + Qa + Qh QCL = Qvp = 11- Qt.
(1) (2)
The total heat loss through each wall (such as the ceiling, windows, wall 1, wall 2, wall 3, wall 4 and the floor) can be calculated by means of Eq. (3) as follows:
Qw = KwAwAT,
(3)
where the temperature difference (AT) was the difference between the inside and neighbor volume temperatures.
The total heat transfer coefficient (Kw) is calculated using Eq. (5) as follows:
1 1 ^ L 1
— = — + > — + —. Kw h j-; t ho
(5)
The heat transfer rate from the product is estimated as follows:
Qp = (mpCp AT ))
cooling '
AT = Tp - Tin.
(6) (7)
Some amount of heat transfer occurred when the door of the cold room was opened to the outside or a place with a different temperature than the inside of the cold room. This situation brought an additional cooling load to the system, which was calculated as [13]:
Qa = HdV(( - hn )pai
(8)
The heat transfer rate from the heat sources (lamps were taken into consideration as heat sources) is determined as follows [13]:
Qh = Q = 3.6 • np.	(9)
1.2 Refrigeration Cycle Calculations
The refrigerant charge rate is calculated as follows:
mm =	(10)
h - h4
The condenser capacity and isentropic compression work of the compressor, respectively, are expressed as follows:
Qco
= m ( h2 - h,), , = ih(h2 -\).
(11) (12)
The COP is used to quantify the performance of refrigeration cycles. It is the ratio of the heat removed from the space to be refrigerated to the work to be consumed to come through this mission. The COP value rises as the power expenditure decreases. It is obvious that this situation will be beneficial to reducing the environmental pollution caused by power production systems. The COP of the refrigeration system's cycle can be determined by the following equation:
AT = T t - T .
out m
(4)
COP = Qevap
(13)
comp
2 RESULTS AND DISCUSSION
Freight transportation has increased its importance during the last century. This tendency is anticipated to have nonstop progress in the coming years because of increasing population and demand. According to recent trends in industry, all transportation methods will experience heavy expansion in the world regarding volume, cost, and energy consumption. The parameters of this work influence the cost and energy consumption issues significantly. In addition, the main advantage of this work is to reduce the pollution from diesel-driven refrigeration equipment in investigated vehicles. Lately, ozone depletion potential (ODP) and global warming potential (GWP) have become the most significant criteria in the improvement of novel refrigerants apart from the refrigerant chlorofluorocarbons (CFCs) and hydro chlorofluorocarbons (HCFCs), both of which have high ODP and GWP, owing to the their role in ozone layer depletion and global warming. Despite their high GWP, substitutes to CFCs and HCFCs such as HFC refrigerants with their zero ODP have been favored for practice in many industrial and domestic requests commonly for a decade. HFC refrigerants also have appropriate qualifications, such as non-flammability, stability, and similar vapor pressure to the refrigerants' CFCs and HCFCs.
Table 2. The detailed specifications of the cold room construction materials of (a) frigoship, (b) railroad car and (c) truck
Polyurethane thickness [m]
Fig. 1. Variation of evaporator capacity with polyurethane thickness
The refrigeration process is the main process of cold-room systems, which provide proper preservation temperatures to preserve the quality of products. The product is usually a type of perishable food within a certain isolated volume. In addition, many types of medical and perishable products are stored in the same way.
Walls and Ceiling of Frigoship					
Material	Stainless Steel Polyurethane				Stainless Steel
L [m]	0.007		0.15		0.007
k [W/mK]	15		0.024		15
ht [W/m2K]			12		
h, [W/m2K]			12		
K [W/m2K]			0.1558		
					
Floor of Frigoship					
Material	Stainless steel	Lean Concrete	Polystyrene	Lean „ . Alum Concrete	
L [m]	0.007	0.1	0.06	0.1	0.007
k [W/mK]	0.78	1.1	0.028	1.1	1.4
hi [W/m2K]			7		
hi [W/m2K]			12		
K [W/m2K]			0.3957		
					
Walls, Ceiling and Floor of Railroad car					
Material	Iron plate	Wood	Polyurethane		Stainless steel
L [m]	0.01	0.02		0.1	0.015
k [W/mK]	45	0.15		0.024	15
h, [W/m2K]			7		
h, [W/m2K]			40		
K [W/m2K]			0.23		
					
Walls, Ceiling and Floor of Truck					
	Glass Fiber	Aluminium	Plastic Foamed	Aluminium	Glass Fiber
Material	Reinforced Polyester				Reinforced Polyester
L [m]	0.005	0.005	0.15	0.005	0.005
k [W/mK]	0.752	205	0.03	205	0.03
h, [W/m2K]			7		
h, [W/m2K]			40		
K [W/m2K]			0.19		
h [kJ/kg
Fig. 2. The ideal vapor-compression refrigeration cycle [20]
1ш
а)
Refrigerant type
b)
Refrigerant type
Fig. 3. Refrigerant flow rate for different refrigerant types of a) frigoship, b) railroad car and c) truck
a)
R12	R134a	R22	R410A
Refrigerant type
b)
Refrigerant type
R12	R134a	R22	R410A
Refrigerant type
Fig. 4. Variation of COP for different refrigerant type of a) frigoship, b) railroad car and c) truck
Table 3. Variation of evaporator capacity with polyurethane thickness
Polyurethane thickness [m]
Qevap [kW]
Frigoship
Railroad car
Truck
0.30
30
3.87
33.43
0.20
34
4.42
35.09
0.15
39
4.96
36.72
0.10
47
6.00
39.90
0.05
70
8.85
48.81
0.03
97
12.14
59.35
Heat transfer rates from the environment are determined by using Eq. (1). Cold rooms' total cooling loads are considered unknown, and heat gain is calculated as 38.73 kW, 7.6 kW and 36.72 kW for frigoship, railroad car and truck by using Eq. (2), respectively. The total heat loss through each wall is estimated by means of Eq. (3). Total heat transfer coefficients of walls were calculated and given with the details in Table 2. The heat transfer rate from the product, air alteration and lamps are estimated
c)
c)
a)
• Refrigerant mass flow rate О Condenser sea water mass flow rate
Table 4. The characteristic points enthalpies of different refrigerants used in the cycle
Characteristic points enthalpies [kJ kg-1]
	Refrigerant	hi	h 2	h3	h4
	R12	342.1	388.1	239.2	239.2
Frigoship	R22	394.8	459.9	249.6	249.6
	R134a	383.4	442.3	256.4	256.4
	R410A	412.1	484.1	266.1	266.1
	R12	344.4	377.6	239.2	239.2
Railroad	R22	397.0	443.7	249.6	249.6
Car	R134a	386.5	429.0	256.4	256.4
	R410A	414.2	465.83	266.2	266.2
	R12	344.43	381.28	247.51	247.5
Truck	R22	397.0	449.06	260.31	260.3
	R134a	386.51	433.48	268.45	268.4
	R410A	414.26	471.59	281.68	281.6
S 0 6 -
40	60
RH [%]
:» бо
"ò 60 A
co &
о
'S 40 g
<D
S 30 -
о Wcomp • Qcond
0	20	40	60	80	100
С)	RH I'M
Fig. 5. Variation of mass flow rate of the sea water with sea water temperature a), the refrigerant and sea water mass flow rate with outside relative humidity b) and the evaporator capacity, the compressor work and the, c) condenser capacity with outside relative humidity
using Eq. (6), Eq. (8) and Eq. (9) for each vehicle, respectively.
It is well known that better insulation results in lower evaporator capacities in cooling applications. Variations of evaporator capacity with polyurethane thickness are given in Table 3 for investigated vehicles. It is found that the evaporator capacity increases as the polyurethane thickness decreases and this is plotted in Fig. 1.
Table 5. Variation refrigerant charge rate, condenser capacity compressor work and COP of different refrigerants for the same evaparator capacity
	Refrigerant	m [kg s-i]	Qevav [kW]	Qcond [kW]	W " comp [kW]
	R12	0.3764	38.73	17.31	56.04
Frigoship	R22	0.2667	38.73	17.35	56.07
	R134a	0.305	38.73	17.97	56.69
	R410A	0.265	38.73	19.06	57.79
	R12	0.0722	7.6	9.99	2.39
Railroad	R22	0.0516	7.6	10.00	2.41
Car	R134a	0.0584	7.6	10.08	2.48
	R410A	0.0514	7.6	10.25	2.65
	R12	0.3789	36.72	50.68	13.96
Truck	R22	0.2686	36.72	50.71	13.98
	R134a	0.3110	36.72	51.33	14.61
	R410A	0.2770	36.72	52.60	15.88
The ideal refrigeration cycle is shown in Fig. 2 and it is assumed that different refrigerants (R12, R22, R134a and R410A) are used in the cycle for conditioning cold rooms. In this investigation, Solkane [14] software is used in order to determine the required refrigerant properties. For the ideal cycle, characteristic point enthalpies, which are shown as 1, 2, 3 and 4, are given in Table 4. It should be noted that evaporation temperatures are assumed as -25 °C, -20 °C and -25 °C for frigoship, truck and railroad car, respectively. In addition, condensation temperatures are assumed as 10 °C, 35 °C and 48 °C, respectively. The refrigerant flow rate, condenser capacity, compressor work and COP are determined by using Eqs. (10) to. (13), respectively, and they are
b)
shown in Table 5. Refrigerant flow rate and refrigerant type comparison of different refrigerants is given in Fig. 3. It is found that the cycle working with R12 has the highest refrigerant flow rate among the other refrigerant types. COP comparison for different refrigerants is given in Fig. 4. It is found that the cycle working with R12 has the highest COP.
In addition, the variation of the mass flow rate of the sea water with sea water temperatures and the refrigerant and sea water mass flow rates with outside relative humidity were given in Figs. 5a and b, respectively. Also, the variation of the evaporator capacity, the compressor work and the condenser capacity with outside relative humidity was investigated and shown in Fig. 5c.
It should be noted that additional explanations, figures and tables on the subject exist in authors' previous publications [15] to [20] in the literature. Moreover, it should be known that time dependent local analyses need experimental works as shown in Ahmet et al. [6]. Otherwise, average values of related parameters can be taken from fundamental sources as in this investigation.
3 CONCLUSION
In this investigation, cold rooms that are located in the frigoship, railroad car and truck were designed. this is mainly aimed at keeping frozen fish (in frigoship) and meat (in railroad car and truck) at -18 °C and -20 °C, respectively. For this purpose, cold rooms' total heat gain is calculated and the ideal vapor-compression refrigeration cycle calculations were performed for R12, R22, R134a and R410A refrigerants.
It can be easily determined that polyurethane thickness is a significant parameter for evaporator capacity, so the effect of polyurethane thickness on evaporator capacity was investigated. It was observed that cooling loads of frigoship, railroad car and truck cold rooms decreased from 97 kW to 30 kW, from 12.14 kW to 3.87 kW and from 59.35kW to 33.43 kW with increase of insulation thickness of the wall from 0.03 m to 0.3 m, respectively.
The refrigerant flow rate, evaporator capacity, condenser capacity, compressor capacity and COP of refrigeration cycle having R12, R22, R134a and R410A as working fluid was determined. R12 showed higher COP when they were compared with other refrigerant types. It should be noted that the COP of refrigeration cycles working with R12 was estimated as 2.24, 2.63 and 3.17 for frigoship, truck and railroad cold rooms, respectively. Unfortunately, R12 cannot be used as an alternative refrigerant because of
its higher ODP and GWP values. Instead of these refrigerants, their alternatives R134a can be used in the applications with its COP values of 2.16, 2.51 and 3.15 for frigoship, truck and railroad cold rooms, respectively.
4 ACKNOWLEDGEMENT
The authors would like to thank the "Research Chair Grant" National Science and Technology Development Agency (NSTDA), the Thailand Research Fund and the National Research University Project (NRU) for the support.
5 NOMENCLATURES
A	area [m2]
Cp	specific heat for the temperature after
	freezing process [J kg-1K-1]
Hđ	air alteration coefficient
hn	enthalpy [kJ kg-1]
hin	inside air enthalpy [kJ kg-1]
h	inner convective heat transfer coefficient
	[W m-2K-1]
ho	outer convective heat transfer coefficient
	[W m-2K-1]
k	thermal conductivity [W m-2K-1]
Kw	total heat transfer coefficient [W m-2K-1]
L	length [m]
l	freezing heat of product [W kg-1]
m	mass flow rate [kg s-1]
mmeat	mass of product [kg]
P	heat emitted by lamp [kW]
Qh	heat transfer from condenser [kW]
Ql	heat transfer from evaporator [kW]
T	temperature [°C or K]
AT	temperature difference [K]
t cooling	cooling time [s]
tl	working time of lamp [hours/day]
V	volume [m3]
Greek symbol:	
P	density [kg m-3]
Subscripts:	
a	air alteration
oeil	ceiling
CL	cooling load
oomp	compressor
oond	condensing/condenser
evap	evaporating/evaporator
f	floor
gen	generation
h	heat sources
in	inside
l	lamp
n	neighbor
out	outside
p	product
w	wall
T	total
6 REFERENCES
[1]	Moureh, J., Menia, N., Flick, D. (2002). Numerical and experimental study of airflow in a typical refrigerated truck configuration loaded with pallets. Computers and Electronics in Agriculture, vol. 34, no. 1-3, p. 25-42, D0l:10.1016/S0168-1699(01)00178-8.
[2]	Tso, C.P., Yu, S.C.M., Poh, H.J., Jolly, P.G. (2002). Experimental study on the heat and mass transfer characteristics in a refrigerated truck. International Journal of Refrigeration, vol. 25, no. 3, p. 340-350, D0I:10.1016/S0140-7007(01)00015-9.
[3]	Hoang, M.H., Laguerre, O., Moureh, J., Flick, D. (2012). Heat transfer modelling in a ventilated cavity loaded with food product: Application to a refrigerated vehicle. Journal of Food Engineering, vol. 113, no. 3, p. 389-398, DOI:10.1016/j. jfoodeng.2012.06.020.
[4]	Moureh, J., Tapsoba, M.S., Flick, D. (2009). Airflow in a slot-ventilated enclosure partially filled with porous boxes. Part 1 - Measurements and simulations in the clear region. Computers & Fluids, vol. 38, no. 2, p. 194-205, DOI:10.1016/j. compfluid.2008.02.006.
[5]	Glouannec, P., Michel, B., Delamarre, G., Grohens, Y. (2014). Experimental and numerical study of heat transfer across insulation wall of a refrigerated integral panel van. Applied Thermal Engineering, vol. 73, no. 1, p. 194-202, DOI:10.1016/j.applthermaleng.2014.07.044.
[6]	Ahmed, M., Meade, O., Medina, M.A. (2010). Reducing heat transfer across the insulated walls of refrigerated truck trailers by the application of phase change materials. Energy Conversation and Management, vol. 51, no. 3, p. 383-392, DOI:10.1016/j.enconman.2009.09.003.
[7]	Tan, H., Li, Y., Tuo, H., Zhou, M., Tian, B., (2010). Experimental study on liquid/solid phase change for cold energy storage of liquefied natural gas (LNG) refrigerated vehicle. Energy, vol. 35, no. 5, p. 1927-1935, DOI:10.1016/j.energy.2010.01.006.
[8]	Cleland, A.C. (1996). Package design for refrigerated food: The need for multidisciplinary project teams. Trends in Food Science & Technology, vol. 7, no. 8, p. 269-271, DOI:10.1016/0924-2244(96)20006-8.
[9]	Morganti, E., Gonzalez-Feliu, J. (2015). City logistics for perishable products. The case of the Parma's Food Hub. Case Studies on Transport Policy, vol. 3, no. 2, p. 120-129, DOI:10.1016/j.cstp.2014.08.003.
[10]	Smale, N. J., Moureh, J., Cortella, G. (2006). A review of numerical models of airflow in refrigerated food applications. International Journal of Refrigeration, vol. 29, no. 6, p. 911930, DOI:10.1016/j.ijrefrig.2006.03.019.
[11]	Duret, S., Hoang, H.M., Flick, D., Laguerre, O. (2014). Experimental characterization of airflow, heat and mass transfer in a cold room filled with food products. International Journal of Refrigeration, vol. 46, p. 17-25, DOI:10.1016/j. ijrefrig.2014.07.008.
[12]	Kou, L., Luo, Y., Ingram, D.T., Yan, S., Jurick, W.M. (2015). Open-refrigerated retail display case temperature profile and its impact on product quality and microbiota of stored baby spinach. Food Control, vol. 47, p. 686-692, DOI:10.1016/j. foodcont.2014.07.054.
[13]	Özkol, N. (1999). Applied Cooling Technology, Chamber of Mechanical Engineers, Istanbul.
[14]	SOLKANE (2015). Refrigerant software, version 8. Solvay Fluor GmbH, Hannover.
[15]	Dalkilic, A.S., Kürekci, N.A., Kincay, O., Wongwises, S. (2013). Fundamental basis and application for cold-room project design: A Turkish case study. Arabian Journal for Science and Engineering, vol. 38, no. 5, p. 1115-1130, DOI:10.1007/ s13369-012-0534-5.
[16]	Dalkilic, A.S., Wongwises, S. (2011). Comparison of various alternative refrigerants for vapour compression refrigeration systems. ASME/JSME 8th Thermal Engineering Joint Conference, p. T10191-T10191-11, DOI:10.1115/ AJTEC2011-44267.
[17]	Kayki, M., Kayaci, N., Celen, A., Dalkilic, A.S., Wongwises, S. (2012). A yacht-type chiller design. ASME 2012 International Mechanical Engineering Congress & Exposition, p. 15331541, DOI:10.1115/IMECE2012-85736.
[18]	Ozturk, M., Goktepe, E., Dalkilic, A.S., Celen, A., Wongwises, S. (2013). Fundamental basis and application of cold-room project design: A case study of fisherman. ASME Summer Heat Transfer Conference.
[19]	Alkan, B., Celen, A., Cebi, A., Dalkilic, A.S., Wongwises, S. (2014). Refrigerated railroad car design for shipping frozen meat using alternative refrigerants. ASME Fluids Engineering Summer Meeting, p. V01CT16A002, DOI:10.1115/ FEDSM2014-21068.
[20]	DalkiliQ, A.S. (2012). Theoretical analysis on the prediction of coefficient of performance (COP) of vapour compression refrigeration systems using various alternative refrigerants. Journal of Thermal Science and Technology, vol. 32, p. 67-79.
Vsebina
Strojniški vestnik - Journal of Mechanical Engineering
letnik 62, (2016), številka 4 Ljubljana, april 2016 ISSN 0039-2480
Izhaja mesečno
Razširjeni povzetki (extended abstracts)
Hubert Kosler, Urban Pavlovčič, Matija Jezeršek, Janez Možina: Adaptivno robotsko razigljevanje
na osnovi izmere 3D položaja in orientacije ulitka z laserskim triangulacijskim senzorjem	SI 35
Shijun Ji, Huijuan Yu, Ji Zhao, Xiaolong Liu, Mingxu Zhao: Ultranatančna obdelava sinusoidne
površine obročev z veliko amplitudo po tehnologiji Slow Tool Servo	SI 36
Andrej Škrlec, Jernej Klemenc: Ocenjevanje parametrov Cowper-Symonds-ovega in Johnson-Cook-
ovega materialnega modela z uporabo Taguchi-jevih ortogonalnih polj	SI 37
Dan Ni, Minguan Yang, Bo Gao, Ning Zhang, Zhong Li: Nestacionarni tok in tlačni impulzi v
črpalki za hladilno vodo jedrskega reaktorja	SI 38
Deniz Duran, Celalettin Karadogan: Neposredna določitev koeficienta Coulombovega trenja s
preizkusi stiskanja valja	SI 39
Ahmet Selim Dalkilig, Ali Celen, Alican Cebi, Somchai Wongwises: Vpliv vrste hladiva in debeline
izolacije na kopenske in ladijske hladilne sisteme	SI 40
Osebne vesti
SI 41
Adaptivno robotsko razigljevanje na osnovi izmere 3D položaja in orientacije ulitka z laserskim triangulacijskim senzorjem
Hubert Kosler1 - Urban Pavlovčič2* - Matija Jezeršek2 - Janez Možina2 1 Yaskawa Motoman, Slovenija 2 Univerza v Ljubljani, Fakulteta za strojništvo, Slovenija
Nastanek igle pri ulitkih predstavlja velik problem, saj lahko poškoduje delavca in zmanjšuje kakovost samega izdelka. Ročno razigljevanje je dolgotrajen, zahteven in drag proces odstranjevanja materiala igle, zato je potreba po avtomatizaciji procesa očitna.
Večina sistemov za robotsko razigljevanje temelji na premikanju obdelovanca z robotsko roko pod orodjem za razigljevanje. Pot robotske roke je običajno določena s postopkom robotskega učenja, kjer operater na referenčnem obdelovancu določi množico točk, prek katerih nato robot premika obdelovanec. Pri tem se predpostavlja, da so vsi nadaljnji obdelovanci enake oblike in se nahajajo na isti lokaciji. Medtem ko je prvi prepostavki pri razigljevanju ulitkov lahko zadostiti, je točen položaj obdelovanca odvisen predvsem od njegovega vpetja v robotsko roko. Običajen pristop je zagotovitev ponovljivega vpenjanja obdelovanca v robotsko roko, kar je v primeru ulitkov zelo težko zagotoviti. Drug pristop, ki je predstavljen v tem prispevku, pa temelji na adaptaciji poti robotske roke glede na izmerjen položaj trenutno vpetega obdelovanca.
Korekcija poti robotske roke temelji na 3D merjenju površine obdelovanca in poravnavi te površine na površino, ki je bila izmerjena na referenčnem obdelovancu. Za 3D merjenje se uporablja laserski triangulacijski profilomer, s katerim najprej izmerimo površino referenčnega obdelovanca. Na tem obdelovancu izvedemo tudi učenje robota. Kasneje, po pritrditvi vsakega naslednjega obdelovanca, izmerimo tudi njegovo površino. To površino poravnamo na površino referenčnega obdelovanca. Rotacije in translacije, ki so bile potrebne za poravnavo, uporabimo za adaptacijo premikov robota. Tako ti ustrezajo položaju obdelovanca, ki je trenutno vpet v robotsko roko.
Za karakterizacijo predstavljene metode smo izvedli dva eksperimenta. Pri prvem smo površino istega ulitka v istem vpetju izmerili desetkrat in tako ugotovili natančnost metode. Ta je bila 0,03 mm v X smeri, 0,01 mm v Y smeri in -0,07 mm v Z smeri. Pri drugem eksperimentu smo ulitek v istem vpetju z robotsko roko premikali od 0,0 mm do 4,0 mm s korakom po 0,5 mm iz začetnega položaja v vsaki izmed osi, na koncu pa še v vseh treh sočasno. Rezultate translacij, izmerjenih s predstavljeno metodo, smo nato primerjali z dejanskimi pomiki robotske roke. Tako smo izmerili točnost korekcije. Rezultati kažejo, da je bila povprečna razlika med izmerjenim in dejanskim položajem robotske roke 0,23 mm s standardnim odklonom 0,12 mm. Povprečne razlike v posameznih oseh so bile 0,05 mm (0,10 mm) za X, -0,02 mm (0,03 mm) za Y in -0,03 mm (0,15 mm) za Z os. Čeprav obdelovanci niso bili namenoma rotirani, so bile izmerjene manjše rotacije, ki pa niso presegale 0,02°.
Izmerjene površine so v povprečju vsebovale 77582 (493) točk. Za poravnavo je sistem potreboval povprečno 1,52 s. Ob redčenju na 6 % izvornega števila točk se je čas poravnave skrajšal na 0,22 s, pri čemer razlike med rezultati 3D poravnave niso bile statistično značilne. Rezultati tako kažejo, da je metoda korekcije položaja dovolj točna za uporabo v sistemih za adaptivno robotsko razigljevanje. Primerna pa je tudi za uporabo na vseh podobnih sistemih, ker se srečujemo z odstopki pri pozicioniranju obdelovancev.
Ključne besede: adaptivne robotske obdelave, razigljevanje, laserska triangulacija, korekcije napake položaja, lokalizacija
Prejeto v recenzijo: 2015-03-01 Prejeto popravljeno: 2015-09-23 Odobreno za objavo: 2015-12-29
Ultranatančna obdelava sinusoidne površine obročev z veliko amplitudo po tehnologiji Slow Tool Servo
Shijun Ji - Huijuan Yu - Ji Zhao - Xiaolong Liu - Mingxu Zhao
Univerza v Jilinu, Fakulteta za tehniške vede in strojništvo, Kitajska
Obdelava prostih površin zahteva napredno tehnologijo. Značilen primer površin prostih oblik so sinusoidni obroči. Sinusoidni obroči velikih dimenzij z amplitudo sinusa makro velikostnega reda, podmikrometrsko natančnostjo oblike in nanometrsko površinsko hrapavostjo imajo pomembno vlogo pri sodobnih optičnih sistemih in natančnih kalibratorjih, njihova izdelava pa bi bila s tradicionalnimi postopki odrezavanja le zelo težavna. Članek predstavlja izdelavo sinusoidne površine obroča s podmilimetrsko amplitudo po postopku enotočkovnega diamantnega struženja (SPDT) s podporo tehnologije STS (Slow Tool Servo).
Predhodne študije so bile usmerjene v amplitudo sinusoidnih površin na mikro-/nanoravni. V naši študiji je bila izdelana sinusoidna površina obroča z amplitudo 0,4 mm. Za izdelavo želene površine so bile analizirane kinematične lastnosti osi stroja za oceno izvedljivosti procesa struženja, pri katerem največja hitrost in največji pospešek linearnih osi ne smeta prekoračiti mejnih vrednosti obdelovalnega stroja. Pot diamantnega orodja je treba za neprekinjeno in gladko površino pripraviti vnaprej, pri čemer se uporabi ravninska spirala. Kontaktne točke orodja se nato pridobijo s projiciranjem ravninske spirale na želeno površino v smeri osi Z. Za ustvarjanje ukazov za premike CNC-stroja je treba pretvoriti rezalne stične točke v točke na rezalnem orodju, oziroma opraviti kompenzacijo zaokrožitev na orodju.
Ker ima vsaka rezalna stična točka edinstveno ukrivljenost vzdolž rezalne poti na obdelovani površini, se lahko kompenzacija orodja namesto z odmikom stičnih točk orodja za vrednost polmera orodja R v smeri Z doseže z odmikom stičnih točk orodja za vrednost polmera orodja R v smeri normalnega vektorja. Na kakovost končnega izdelka močno vpliva geometrija diamantnega orodja in pravilna izbira rezalnega orodja je zato pomembna za načrtovanje obdelave.
Ukrivljenost orodja se spreminja v različnih točkah in da ne bi prišlo do motenj, morata biti zaokrožitev rezalnega orodja in kot konice v vseh rezalnih točkah manjša od minimalne vrednosti polmera zaokrožitve in kota loka. Naslednji pomemben parameter pri izbiri geometrije orodja je prosti kot, ki mora biti v vseh rezalnih točkah večji od največje vrednosti kota nagiba.
Uporabljeno je bilo diamantno orodje s polmerom zaokrožitve 0,506 mm, kotom konice 120° in prostim kotom 10°, ki izpolnjuje opredeljene kriterije. V eksperimentih je bil uporabljen stroj Nanoform 250. Za meritev kakovosti površine obdelane komponente sta bili izmerjeni in analizirani površinska hrapavost in natančnost oblike obdelanega sinusoidnega obroča z mikroskopom za velike povečave (Contour GT) in sistemom za precizno merjenje oblik (Taylor Hobson).
Za vrednotenje kakovosti površine je bila opravljena linearna meritev hrapavosti. Površinska hrapavost (Ra) v smeri meritve je približno 7,5 nm. Opravljena je bila tudi ploskovna meritev hrapavosti površine (Sa), ki znaša 4,6 nm. Natančnost oblike PV sinusoidne površine obroča, obdelanega s tehnologijo STS, pa znaša 0,274 ^m. Rezultati eksperimentov razkrivajo, da je bila s teoretično analizo učinkovito izboljšana natančnost obdelave. Ključne besede: sinusoidna površina obroča, STS, kinematične lastnosti, ultranatančna obdelava
Ocenjevanje parametrov Cowper-Symonds-ovega in Johnson-Cook-ovega materialnega modela z uporabo Taguchi-jevih ortogonalnih polj
Andrej Škrlec* - Jernej Klemene
1 Univerza v Ljubljani, Fakulteta za strojništvo, Slovenija
Za zmanjšanje stroškov raziskav in razvoja izdelkov (R&R), se v začetnih fazah razvojnega procesa uporabljajo kompleksna numerična vrednotenja in simulacije. Z njimi lahko napovemo, kako se bo izdelek obnašal v določenih obratovalnih pogojih. Za čim bolj točno napovedovanje obnašanja izdelka je potrebno izdelati 3D numerični model izdelka, poznati obratovalne pogoje in materialne lastnosti izdelka. Te lastnosti upoštevamo v konstitutivnih materialnih modelih, ki jih uporabljamo v numeričnih simulacijah. Zelo pogosto se v numeričnih simulacijah uporabljajo konstitutivni materialni modeli, pri katerih vrednosti njihovih parametrov ne moremo enostavno določiti. Zato je potrebno z izvajanjem eksperimenta in optimizacijskim postopkom poiskati ustrezne vrednosti parametrov konstitutivnih materialnih modelov. Optimalne vrednosti teh parametrov so tiste, pri katerih dobimo enak odziv izdelka kot v eksperimentu.
V članku je predstavljena metodologija za ocenjevanje parametrov materialnih modelov, ki opisujejo obnašanje materiala pri povišanih hitrostih specifičnih deformacij. Postopek ocene parametrov materialnih modelov, ki opisujejo vpliv hitrosti specifičnih deformacij, je osnovan na preizkusu vtiska kroglice v pločevinast preizkušanec. V ta namen smo izdelali napravo, ki jo lahko namestimo na obstoječe preizkuševališče. Na osnovi eksperimentalnih rezultatov smo ocenjevali parametre Cowper-Symonds-ovega in Johnson-Cook-ovega materialnega modela, ki se ju največkrat uporablja za numerične simulacije izdelkov, ki so podvrženi obremenitvam pri povišanih hitrostih specifičnih deformacij. Za oceno optimalnih vrednosti parametrov materialnih modelov smo med seboj povezali numerične simulacije in Taguchi-jeva ortogonalna polja.
Taguchi-jeva ortogonalna polja smo uporabili za določitev posameznih kombinacij vrednosti parametrov omenjenih materialnih modelov. S kombinacijami vrednosti parametrov smo izvedli numerične simulacije in njihove rezultate primerjali z rezultati eksperimenta. Za primerjavo rezultatov in posledično izbor najboljše kombinacije parametrov smo določili cenilno funkcijo. Ker so območja možnih vrednosti parametrov materialnega modela lahko velika in ker smo želeli raziskati večji prostor možnih rešitev, je bilo potrebno izbrati več vrednosti posameznega parametra. Tako bi dobili zelo veliko število kombinacij, kar je časovno neustrezno, saj je potrebno izvesti simulacijo za vsako kombinacijo rešitev.
S Taguchi-jevim ortogonalnih poljem lahko število možnih kombinacij zmanjšamo tako, da ne zmanjšamo števila možnih vrednosti posameznega parametra. Taguchi-jevo ortogonalno polje smo najprej uporabili na širšem območju raziskovalnega prostora. Ker je bilo na širšem območju večje število optimalnih kombinacij, smo nato Taguchi-jevo ortogonalno polje uporabili še na zmanjšanem območju. Tako dobljene vrednosti parametrov Cowper-Symonds-ovega in Johnson-Cook-ovega materialnega modela smo primerjali z vrednostmi iz literature. Izkaže se, da so vrednosti iz predstavljene metode primerljive vrednostim iz literature, kljub temu da so bile dobljene z drugačnimi preizkusi in ob dejstvu, da so bile hitrosti specifičnih deformacij pri našem eksperimentu višje za faktor 10 ali več od referenčnih podatkov. Predstavljena metoda ocenjevanja parametrov s Taguchi-jevim planom preizkusov je časovno efektivna še posebej, če imamo večje število parametrov materialnih modelov z relativno širokimi območji možnih vrednosti posameznega parametra.
Ključne besede: obnašanje materiala, hitrost specifične deformacije, metoda končnih elementov, eksplicitna dinamika, Taguchi-jeve ortogonalne matrike
Nestacionarni tok in tlačni impulzi v črpalki za hladilno vodo jedrskega reaktorja
Dan Ni* - Minguan Yang - Bo Gao -Ning Zhang -Zhong Li
Univerza Jiangsu, Fakulteta za energetiko in elektroenergetiko, Kitajska
Črpalke za hladilno vodo jedrskega reaktorja (RCP), ki skupaj s sredico reaktorja in generatorji pare za prenos toplotne energije predstavljajo srce jedrskega reaktorja, poganjajo kroženje hladilne vode v primarnem krogu. Te črpalke so edini rotacijski del v jedrskem otoku, zato spadajo med komponente prvega razreda jedrske varnosti. Poleg tega so tudi glavni porabnik energije v jedrski elektrarni ter morajo zagotavljati dolgotrajno delovanje brez prekinitev in motenj. Notranji tokovi v črpalki za hladilno vodo jedrskega reaktorja z mešanim tokom so eden glavnih problemov pri razvoju teh črpalk za velike tlačne vodne reaktorje. Nestacionarni tokovi v črpalki, še posebej v določenih predelih, so zelo pomembni za varnostno analizo jedrskega reaktorja, saj lahko povzročijo vibracije, ki ogrožajo integriteto črpalke. Pri konstruiranju tovrstnih črpalk je zato nujna celovita analiza in napovedovanje tlačnih impulzov zaradi interakcij med rotorjem in statorjem. Za izboljšanje učinkovitosti in stabilno obratovanje črpalk za hladilno vodo jedrskih reaktorjev je nujen podroben opis zapletenih nestacionarnih struktur notranjih tokov.
Pri teh črpalkah gre za posebno konstrukcijo s sferičnim ohišjem, ki določa značilne in kompleksne tokovne strukture v črpalki. Članek predstavlja študijo notranjih nestacionarnih tokov v modelni črpalki RCP z mešanim tokom pri nazivnih pogojih na podlagi metode simulacije velikih vrtincev (LES) in eksperimentov. Merilna mesta za merjenje tlačnih impulzov in nestacionarnega toka so bila postavljena na vhodu kolesa, v kanalih difuzorja in na sferičnem ohišju. Podatki o tlačnih impulzih in zapletenih tokovnih strukturah so bili nato združeni za opis nestacionarnega toka v črpalki.
Notranje tokovne strukture so težko določljive in zapletene zaradi posebnosti sferičnega ohišja teh črpalk. Raziskava je pokazala, da se tokovne strukture pri nazivnem pretoku v različnih kanalih občutno razlikujejo in da so tesno povezane s položajem kanala difuzorja glede na šobo. V kanalih difuzorja v bližini šobe prihaja do ločevanja toka in do povratnih tokov. Na predel levo od praznilne šobe ohišja pomembno vplivajo intenzivne interakcije med rotorjem in statorjem. Desno in na sredini šobe ohišja pa so nestacionarne tokovne strukture bolj zapletene kot drugje zaradi obsežnega ločenega toka (prehodni in krožni tok) in močnega vrtinčenja. Ugotovljeno je bilo, da imajo vrtinčni in tlačni spektri praktično enake vzbujalne frekvence. Tlačni impulzi v določenih predelih so tako odvisni od odlepljanja vrtincev na zadnjem robu lopatice difuzorja modelne črpalke za hladilno vodo jedrskega reaktorja.
Nadaljnje raziskave bodo usmerjene v odvisnosti med nestacionarnimi tokovnimi strukturami in tlačnimi impulzi v zadevnih črpalkah v vseh delovnih pogojih, uporabljene pa bodo tehnike merjenja tlačnih impulzov in PIV.
Ključne besede: črpalka za hladilno vodo jedrskega reaktorja, nestacionaren tok, tlačni impulzi, sferično ohišje, vrtinčenje, obsežen ločeni tok, nestacionarne tokovne strukture
SI 38
*Naslov avtorja za dopisovanje: Univerza Jiangsu, Fakulteta za energetiko in elektroenergetiko, Zhenjiang 212013, Kitajska, nxm0424@163.com
Neposredna določitev koeficienta Coulombovega trenja s preizkusi stiskanja valja
Deniz Duran* - Celalettin Karadogan
Center odličnosti za preoblikovanje kovin pri univerzi Atilim, Turčija
Razvita je bila nova metoda za neposredno določitev koeficienta Coulombovega trenja pri aplikacijah kovanja v hladnem s preizkusi stiskanja valja. Pri klasičnem pristopu se meri deformirana neenakomerna oblika stisnjenega preizkušanca in meritve se nato primerjajo z rezultati kalibracijskih simulacij, ki pokrivajo različne koeficiente trenja. Klasični pristop je zelo zamuden in zahteva vnaprejšnje poznavanje krivulj tečenja materiala. Težave nastopijo zaradi nehomogene porazdelitve spremenljivk polja, ki so posledica trenja na stiku med ploščo in preizkušancem. Nova metoda odpravlja potrebo po podrobnem poznavanju krivulje tečenja ter lahko zamenja kalibracijske simulacije za vse materiale in koeficiente trenja s preprostim polinomskim metamodelom.
Metamodel je sestavljen na podlagi osnosimetričnih simulacij s končnimi elementi za preizkus stiskanja valja premera 10 mm in višine 15 mm. Realistično območje vedenja materiala in koeficientov trenja je bilo izbrano po podatkih iz literature. Izbrana območja za K in n v potenčni formuli ter p kot koeficient Coulombovega trenja so bila 250 < K < 1250 (MPa), 0,05 < n < 0,5 in 0,05 < p < 0,3. Po velikem številu simulacij v tem območju parametrov se je izkazalo, da je premaknjena stična površina ICA pri stiskanju na polovično višino, ki predstavlja lateralni del površine v stiku s ploščami, le funkcija koeficienta trenja in eksponenta deformacijskega utijanja. Ta omejena odvisnost omogoča opis metamodela z enostavnimi polinomskimi relacijami.
Enačba (2) podaja koeficient Coulombovega trenja p v odvisnosti od širine ICA u pri znanih koeficientih an, bn, cn in dn. Ti koeficienti so funkcije eksponenta deformacijskega uüjanja n, skladno z enačbo (3).
Zgornji metamodel zahteva le eksperimentalno določitev eksponenta deformacijskega uüjanja materiala in širine ICA na stisnjenem valju u. Eksponent deformacijskega utrjanja je treba določiti z mehanskim preizkusom, meritev u pa je preprosta in natančna s tehnikami digitalne korelacije posnetkov (DIC). ICA po kompresijskem testu se lahko navadno določi brez težav na podlagi razlike kontrastov. Vrednost u se izračuna z notranjimi in zunanjimi krogi na končni ploskvi stisnjenih valjev, ki pa niso idealni zaradi nehomogenosti materiala in stične površine. Za te površine je zato treba izračunati ekvivalentne polmere. Območja se v ta namen od zunaj obrišejo z zlepki in določijo se polmeri ustreznih ekvivalentnih krogov. Ekvivalentna vrednost u je razlika med obema ekvivalentnima polmeroma. MATLAB-ova skripta za kvazisamodejni izračun vrednosti u je podana v prilogi. Skripta uvozi sliko stisnjenega preizkušanca in dela z milimetrsko razdaljo umerjanja. Najprej se poiščejo slikovne točke, ki ustrezajo končnim točkam razdalje umerjanja. V naslednjem koraku se po zunanjosti notranjega kroga (t. j. po površini, ki ustreza končni ploskvi začetnega valja) in po zunanjosti zunanjega kroga (t. j. po končni skupni kontaktni površini) razporedi ustrezno število znamenj. Iz koordinat postavljenih znamenj se periodično konstruirajo zaprti kubični zlepki, iz katerih se nato izračuna površina znotraj teh sklenjenih krivulj.
Predlagana metoda torej izkorišča uveljavljeni inverzni pristop in vse njegove možnosti. Praktična uporaba je neposredna in preprosta ter odpravlja potrebo po inverzni numerični analizi za širok razpon materialov. Ko je območje koeficienta trenja omejeno na razpon 0,1 < p < 0,25, je napaka zaradi polinomske interpolacije omejena največ na 5 %, povprečna napaka pri reprodukciji osnovnih podatkov pa ne presega 1,5 %. Točkovna napaka reprodukcije v celotnem območju 0,05 < p < 0,3 je največ 10 %, pri čemer je celotna povprečna napaka največ 2 %. Glavna prednost predstavljene metode je v njeni praktičnosti, saj za karakterizacijo trenja ne zahteva nobene posebne opreme. Pri karakterizaciji trenja se lahko uporabijo tudi kompresijski preizkusi za vrednotenje krivulj tečenja. Koeficient trenja, pridobljen po predlagani metodi, je uporaben tudi za inverzno vrednotenje kompresijskega testa pri določevanju korektur. Metoda uporablja običajne fotografske tehnike in je primerna za različne materiale. Rezultati eksperimentalne študije dokazujejo, da je predlagana metoda ustrezna in dovolj točna zamenjava za zamudno inverzno analizo.
Ključne besede: koeficient Coulombovega trenja, preizkus s stiskanjem valja, spremenjena kontaktna površina, digitalna analiza posnetkov
*Naslov avtorja za dopisovanje: Center odličnosti za preoblikovanje kovin pri univerzi Atilim, Turčija, dnz.duran@gmail.com
39
Vpliv vrste hladiva in debeline izolacije na kopenske in ladijske hladilne sisteme
Ahmet Selim Dalkilig1* - Ali Celen1 - Alican Qebi1 - Somchai Wongwises2 1 Tehniška univerza Yildiz, Fakulteta za strojništvo, Turčija 2 Tehniška univerza Thonburi kralja Mongkuta, Tehniška fakulteta, Tajska
Transportna vozila imajo pomembno vlogo pri ohranjanju kakovosti ter podaljševanju roka trajanja svežih, pokvarljivih in zamrznjenih izdelkov. Pri vsakem procesu hlajenja je bistveno zagotavljanje dovolj nizke temperature za minimalno metabolno in mikrobiološko razgradnjo. Z drugimi besedami: projektanti snujejo za temperature pod lediščem, ki preprečujejo rast bakterij. Projektanti morajo pri snovanju hladilnih procesov upoštevati več pomembnih podrobnosti. Vedeti morajo npr. kateri kompresor in kondenzator bosta najprimernejša za sistem hladilnice. Upoštevati morajo tudi učinkovitost hladilnice, čas transporta od hladilnice do prodajnih mest, vrste izdelkov, ki se lahko skladiščijo skupaj, gradbene materiale in sezonske spremembe temperature na lokaciji.
V tej raziskavi so bili projektirani hladilna ladja s 3391 m3 skladiščnega prostora za 57.000 kg zamrznjene ribe, železniški vagon s 93 m3 skladiščnega prostora za 4080 kg zamrznjenega mesa in tovornjak z 113 m3 skladiščnega prostora za 23.000 kg zamrznjenega mesa. V članku so podane glavne dimenzije vozil. Zamrznjena riba in meso se hrani pri temperaturi -18 oz. -20 °C.
Hladilna obremenitev hladilnice je določena kot vsota toplote, ki prehaja skozi stene, iz izdelkov, krožečega zraka in toplotnih virov. Članek navaja vse pripadajoče enačbe, upoštevani pa so tudi toplotni prejemki neznanega izvora. Celotna hladilna obremenitev sistema je bila povečana za 10 %, kar ustreza stopnji odvajanja toplote iz hladilnice na uparjalniku. Izračunan je masni tok hladiva, moč kondenzatorja in izentropno kompresijsko delo kompresorja, za kvantifikacijo učinkovitosti hladilnih ciklov pa je uporabljeno hladilno število (COP). Hladilno število je razmerje med toploto, odvedeno iz hladilnice, in vloženim delom. Vrednost COP raste z zmanjševanjem porabe električne energije. Takšen razvoj pripomore k zmanjšanju onesnaževanja okolja, ki ga povzročajo sistemi za proizvodnjo električne energije.
Končno so bile projektirane hladilnice za shranjevanje hrane v hladilni ladji, železniškem vagonu in tovornjaku. Upoštevan je bil tudi vpliv hladiv in debeline izolacije hladilnih sistemov v hladilnicah. V analizah sta bili uporabljeni hladivi R12 in R22 ter njuni alternativi R134a in R410A. Določena je bila moč uparjalnika, moč kondenzatorja, delo kompresorja, pretok hladiva in hladilno število (COP) projektiranih hladilnih sistemov posameznih vozil. Ugotovljeno je bilo, da imata hladivi R134a in R410A nekoliko manjše hladilno število in zahtevata več dela na kompresorju za temperaturo kondenzacije kot hladivi R12 in R22. Hladilnice v hladilni ladji, tovornjaku in vagonu delujejo z vrednostmi COP 2,24, 2,63 in 3,17. Predlagano hladivo R134a se lahko uporablja v hladilnicah hladilne ladje, tovornjaka in vagona z vrednostmi COP 2,16, 2,51 in 3,15. Preučen je bil tudi vpliv debeline izolacije stene na moč uparjalnika. Ob povečanju debeline izolacije sten iz 0,03 m na 0,3 m se je hladilna obremenitev pri hladilni ladji, železniškem vagonu in tovornjaku zmanjšala za 69, 68 oz. 43 %. Ključne besede: hlajenje, alternativna hladiva, hladilno število, hladilnica, izolacija, kondenzator, uparjalnik
DOKTORSKE DISERTACIJE, ZNANSTVENO MAGISTRSKO DELO
Na Fakulteti za strojništvo Univerze v Ljubljani so obranili svojo doktorsko disertacijo:
•	dne 7. marca 2016 Marko JERMAN z naslovom: »Fenomenološko modeliranje lednega abrazivnega vodnega curka v industrijskih procesih rezanja« (mentor: prof. dr. Mihael Junkar, somentor: doc. dr. Henri Orbanić);
Kadar pri postopku obdelave z abrazivnim vodnim curkom (AVC) mineralni abraziv zamenjamo z ledenimi zrni, ustvarimo ledni abrazivni vodni curek (LAVC). Glavna prednost tehnologije je čistost postopka in nizke temperature obdelave. V okviru disertacije smo razvijali prototip stroja za obdelavo z LAVC in na različnih stopnjah razvoja prototipa na njem izvajali eksperimente, s katerimi smo popisovali fenomene procesa. Na ta način smo določevali mesta segrevanja vode v sistemu in merili, kako ohlajanje vode vpliva na učinkovitost in kakovost obdelave z vodnim curkom. Na podlagi rezultatov smo z metodo dimenzijske analize postavili semi-empirični model za napovedovanje učinkovitosti in kakovosti obdelave z upoštevanjem temperature curka. Z uporabo numerične analize smo ugotavljali temperaturne razmere znotraj vodne šobe, na podlagi katerih smo sklepali o temperaturnih pogojih, katerim so izpostavljena ledena zrna. Model smo validirali s pomočjo eksperimentov in literature. V sklopu raziskav smo razvili optično metodo za zasledovanje razpada vodnega curka, na podlagi katere smo ocenjevali vpliv spreminjanja parametrov na kakovost curka. Izdelali smo napravo za izdelavo ledenih zrn, ki smo jih uporabili kot abraziv. Na podlagi raziskav smo izdelali prototip stroja za rezanje z LAVC in izvedli prve poskuse rezanja;
•	dne 8. marca 2016 Eva OBLAK z naslovom: »Nano-tribološke lastnosti mejnih površinskih filmov« (mentor: prof. dr. Mitjan Kalin);
Razvoj novih aplikacij in sistemov je povezan s funkcionalnimi, ekonomskimi in ekološkimi izboljšavami, ki jih lahko dosežemo z novimi nekonvencionalnimi materiali (npr. DLC prevleke) in okolju prijaznejšimi mazivi. Za učinkovito uporabo le teh je potrebno novo razumevanje, do kakšnih interakcij in tvorjenja tribofilma prihaja med njimi. Prav lastnosti nastalega tribofilma imajo predvsem v sistemih, ki obratujejo v režimu mejnega mazanja, ključen vpliv na trenje. V raziskavi smo se osredotočili na kontakte jeklo/jeklo, jeklo/a-C:H in jeklo/SiDLC, mazane s tremi olji z različno reaktivnostjo aditivov in okoljsko sprejemljivostjo (olja razreda E6 in
E7), z mineralnim baznim oljem z aditivi ZDDP in s čistim mineralnim baznim oljem. Preučevali smo nanomehanske in nanotribološke lastnosti tribofilmov, ki so še vedno zelo slabo ovrednotene. Z meritvami lastnosti tribofilma v različnih časovnih intervalih smo ovrednotili rast in razvoj tribofilma. Z mikroskopom na atomsko silo (AFM) smo merili lastnosti tribofilma, kot so: morfologija, nanohrapavost, nanotrenje, nanotogost, debelina filma in adhezija. Z nanovtiskovalcem pa smo določili njihovo nanotrdoto. Nastanek tribofilma smo potrdili na vseh izbranih površinah. Poleg tega smo pokazali, da je rast tribofilma močno odvisna od izbranih aditivov in površine. Ovrednotili smo odvisnosti med lastnostmi tribofilmov na nanoskali in makrotrenjem. Odvisnosti so po krajšem času drsenja veliko bolj očitne, kot po daljšem času drsenja, ko se razvije bolj kompleksen tribofilm in vplivajo na makrotrenje še drugi neizmerjeni parametri in prihaja do medsebojnih vplivov med parametri. DLC prevleka dopirana s silicijem se je izkazala kot primernejša izmed obeh DLC prevlek, saj je zaradi vsebnosti silicija sposobna tvoriti učinkovitejše tribofilme, s katerimi dosegamo dobre tribološke rezultate, tudi z olji z zmanjšano vsebnostjo aditivov;
• dne 11. marca 2016 Martin PETKOVŠEK
z naslovom: »Kavitacija v termično občutljivih kapljevinah« (mentor: izr. prof. dr. Matevž Dular, somentor: prof. dr. Branko Širok);
Doktorsko delo obravnava fizikalni pojav kavitacije v termično občutljivih kapljevinah. Kavitacija je pojav, okarakteriziran s prehodom iz kapljevitega v plinasto agregatno stanje, ter nazaj v homogeno kapljevino pri približno konstantni temperaturi. Kavitacija navadno nastane na mestu lokalnega padca tlaka v bližini kavitacijskega jedra. Ko mehurček raste, latentna toplota iz okolice prehaja v mehurček, kar ustvari termično mejno plast. Posledica je lokalno znižanje temperature kapljevine, ki rezultira v padcu parnega tlaka. Pojav je poznan kot termični zamik ali termodinamski učinek in ima pri razvoju kavitacije zaviralno vlogo - nižji parni tlak odloži nadaljnji razvoj mehurčka, saj je zdaj za njegovo rast potreben večji tlačni padec. Termodinamski učinki se navadno lahko zanemarijo v kapljevinah, katerih kritična temperatura je precej višja od delovne temperature oz. termodinamski učinki postanejo izraziti, ko je delovna temperatura kapljevine v bližini kritične temperature kapljevine -na primer pri kriogenih kapljevinah.
V	okviru doktorskega dela je bil fenomen, termodinamski učinek, eksperimentalno opazovan v vodi temperature približno 100 °C, kjer so se v kavitirajočen toku merila temperaturna polja s pomočjo infrardeče termografije. Eksperimentalno pridobljeni rezultati bodo pripomogli k boljšemu razumevanju samega pojava in k razvoju ter validaciji naprednih numeričnih modelov, ki upoštevajo termodinamski učinek v kavitirajočem toku. Pravilna numerična napoved termodinamskih učinkov kavitacije je torej velik temeljni in aplikativni problem;
•	dne 11. marca 2016 Eniko PETER z naslovom: »Izboljšanje kakovosti procesa izdelave z metodami večkriterijskega optimiranja« (mentor: doc. dr. Davorin Kramar);
Doktorska disertacija obravnava razvoj in aplikativno potrditev modela za nadzor kakovosti in optimiranje tehnologije izdelave. Na osnovi hipoteze, da je z razvojem omenjenega modela možno doseči stabilno tehnologijo izdelave, je bila obravnavana tehnologija izdelave batne izvrtine. Model za nadzor kakovosti in optimiranje tehnologije izdelave je primeren tako za serijsko kot maloserijsko proizvodnjo. Model vsebuje metode statističnega nadzora procesa, prepoznavanja vzrokov in posledic, načrtovanja eksperimentov, optimiranja in ekonomske analize. Z aplikativno uporabo modela je bil ugotovljen in minimaliziran negativen vpliv neproduktivnih operacij na kakovost produkta. V okviru izhodiščnih raziskav je bilo določeno tehnološko okno procesa povrtavanja. V nadaljevanju je bila tehnologija izdelave batne izvrtine empirično modelirana, večkriterijsko optimirana, optimiran model pa implementiran v redno proizvodnjo. Ekonomska analiza je v zaključku potrdila učinkovitost ter ekonomičnost modela za nadzor kakovosti in optimiranje tehnologije izdelave;
•	dne 30. marca 2016 Marko BOJINOVIĆ z naslovom: »Računalniško simuliranje procesov toplotne obdelave podevtektoidnega jekla z diodnim laserjem« (mentor: prof. dr. Boris Štok, somentor: doc. dr. Nikolaj Mole);
V	disertaciji je predstavljen numerični model in validirana računalniška simulacija laserske toplotne obdelave podevtektoidnega jekla 50CrV4. Za numerično reševanje matematičnega modela je uporabljena metoda končnih elementov v komercialnem programskem paketu ABAQUS nadgrajenim z lastnimi podprogrami. Računalniška
simulacija je validirana v treh soodvisnih sklopih, to je validacija temperaturnega polja na površini obdelovanca, validacija trdote po globini obdelovanca in validacija spremembe geometrije obdelovanca. S pomočjo validirane računalniške simulacije je prikazan vpliv tehnoloških parametrov pri laserskih toplotnih obdelavah na temperaturno polje, trdoto ter deformacijsko-napetostno polje na obdelovancu v obliki kvadra velikosti 100 mm x 100 mm x 15 mm. V delu je še posebej predstavljen vpliv hitrosti segrevanja in vpliv začetne mikrostrukture na kinetiko faznih transformacij in trdote pri laserski toplotni obdelavi podevtektoidnega jekla; • dne 31. marca 2016 Gorazd KRESE z naslovom: »Določitev ravnovesne temperature kot energijske karakteristike stavbe na podlagi analize rabe energije« (mentor: prof. dr. Vincenc Butala, somentor: doc. dr. Matjaž Prek);
Delo obravnava razvoj semiempiričnega modela za določitev ravnovesne temperature stavbe z upoštevanjem vpliva sončnega obsevanja in vlažnosti zraka. Model temelji na predpostavljeni funkcijski odvisnosti med toplotno obremenitvijo stavbe in močjo sistema za vzdrževanje parametrov notranjega okolja. Sestavlja ga 11 prostih parametrov, ki so določeni s pomočjo optimizacijskega postopka na osnovi podatkov o odjemni električni moči za hlajenje stavbe in izbranih meteoroloških spremenljivk za lokacijo stavbe. Pravilnost delovanja razvitega modela je bila preverjena na vzorcu realnih stavb kot tudi na umetno generiranih setih podatkov. Dokazano je bilo, da razvita metodologija za določitev in analizo ravnovesne temperature omogoča obravnavo le-te kot energijske karakteristike stavbe.
ZNANSTVENO MAGISTRSKO DELO
Na Fakulteti za strojništvo Univerze v Ljubljani je z uspehom zagovarjal svoje magistrsko delo:
• dne 22. marca 2016 Andrej POD NAR z naslovom: »Korekcija profila rotorske lopatice Kaplanove turbine za zmanjševanje kavitacije na lopaticah turbine« (mentor: prof. dr. Branko Širok)
V spomin akademiku prof. dr. Janezu Pekleniku In Memoriam Academician prof. Janez Peklenik, Ph.D.
V 90. letu starosti nas je zapustil akademik profesor dr. Janez Peklenik, zaslužni profesor Univerze v Ljubljani.
Rodil se je leta 1926 v Tržiču. Domače obrtniško okolje mu je zbudilo zanimanje za strojništvo, ki ga je vodilo skozi vse življenje. V šoli je dobro napredoval, vendar mu je pričetek 2. svetovne vojne leta 1941 onemogočil nadaljevanje šolanja na gimnaziji. Zato se je izučil za orodjarja v tovarni letalskih delov LGW v Kranju, iz katere se je kasneje razvila Iskra. To je bila prva strokovna stopnica, ki jo je uspešno prestopil. Leta 1944 se je kot mlad borec Prešernove brigade in kasneje Jeseniško-Bohinjskega odreda pridružil NOV Na to obdobje svojega življenja je bil vedno zelo ponosen.
V povojnih letih je najprej dopolnil zamujeno gimnazijsko izobrazbo in opravil veliko maturo. Šolanje je nadaljeval s študijem strojništva na Univerzi v Ljubljani, kjer je z odliko diplomiral leta 1954. Že v času študija je izkazal nagnjenost k raziskovalnemu delu, za kar je dvakrat prejel Prešernovo nagrado. Med študijem si je tudi nabiral izkušnje kot konstruktor v industriji obdelovalnih strojev.
Po končanem študiju je odšel v Laboratorij za obdelovalne stroje na Tehniški visoki šoli (RWTH) v Aachnu, ki je vodilna raziskovalna ustanova s področja proizvodnega strojništva. Po dveh letih raziskovanja je leta 1957 z odliko doktoriral pod mentorstvom prof. Opitza. Z delom je nadaljeval na isti inštituciji kot znanstveni sodelavec. Leta 1961 je habilitiral z delom »Natančnostna vprašanja v kibernetizaciji proizvodnje« in bil izvoljen v docenta. Naslednje leto je bil kot gostujoči izredni profesor povabljen na Carnegie-Melon Univerzo v Pittsburghu, ZDA.
Leta 1964 se je preselil v Birmingham, kjer je bil na tamkajšnji univerzi izvoljen v naziv rednega profesorja in kjer je osnoval prvo katedro za obdelovalne sisteme na svetu. Isto leto je bil izvoljen za rednega profesorja
Academician Janez Peklenik, Ph.D., Professor Emeritus of the University of Ljubljana, has passed away at 89 years of age.
Born in the town of Tržič in 1926, Peklenik was inspired by the trade and craftsmen in his home environment to pursue mechanical engineering, which drove him throughout his life. He was a keen student, but the beginning of WWII in 1941 prevented him from attending secondary school at a gymnasium. Undeterred, he became a qualified tool maker at the LGW (Luftfahrgeratewerke Krainburg) airplane parts factory in Kranj, a company that later evolved into Iskra. This was the first successful step in his professional career. In 1944 he joined the national liberation struggle as a young fighter in the Prešeren Brigade and later in the Jesenice-Bohinj Detachment, and this was a part of his life he was very proud of.
In the years after the war he finally received his secondary education and took his matriculation exams at the highest level. He then went on to continue his studies in mechanical engineering at the University of Ljubljana, where he graduated cum laude in 1954. Already during his years at the university he showed interest in research work, which twice earned him the Prešeren Award. While studying he also gained experience as a technician in the machine tooling industry.
After graduating he went to the Laboratory for Machine Tools at the RWTH technical college (Rheinisch-Westfälische Technische Hochschule) in Aachen, the leading research institute for production engineering. He earned his doctorate in 1957 after two years of research, under the mentorship of Professor Opitz. Employed as a researcher, he continued working at the same institution. Peklenik earned his habilitation in 1961 with the title "Questions of Precision in the Cybernetisation of Production", and was promoted to assistant professor. The next year he received an invitation to lecture as a guest associate professor at Carnegie-Melon in Pittsburgh, USA.
In 1964 he moved to Birmingham, where he was appointed as a full professor and where he established the first department for manufacturing systems in the world. The same year he was elected full professor at the University of Illinois in Urbana-Champaign.
While he was abroad, Peklenik still worked with the University of Ljubljana. In 1972 he returned home
tudi na Univerzi v Illinoisu v kraju Urbana-Champain, ZDA.
Že med službovanjem v tujini je sodeloval z Univerzo v Ljubljani. Leta 1972 pa se je za stalno vrnil v domovino, kjer je na Fakulteti za strojništvo Univerze v Ljubljani kot redni profesor osnoval Katedro za tehnično kibernetiko, obdelovalne sisteme in računalniško tehnologije in istoimenski laboratorij. Vpeljal je vrsto novih, sodobnih predmetov ter modularni način študija. Osnoval je študijsko smer Mehatronika in Podiplomsko šolo za avtomatizacijo in proizvodno kibernetiko. Pod njegovim mentorstvom je študij zaključila preko sto diplomantov, preko petdeset magistrov in preko trideset doktorjev znanosti.
Med službovanjem na Fakulteti za strojništvo je opravljal vrsto vodilnih funkcij. Dolga leta je bil predstojnik katedre. V letih 1974 do 1976 je bil dekan Fakultete za strojništvo, v letih 1987 do 1990 pa rektor Univerze v Ljubljani.
Po upokojitvi leta 1996 je še vrsto let sodeloval z laboratorijem in bil angažiran na številnih drugih področjih. Z leti pa je ta zavzetost počasi ugašala in zadnja leta je preživel v krogu svojih najbližjih.
Prof. Peklenik je predstavil rezultate svojega dela v obsežnem opusu, ki ga sestavljajo številne vrhunske in mnogokrat citirane objave in vrsta domačih in mednarodnih patentov. Njegovo znanstveno delo obsega zelo širok spekter raziskovalnih dognanj na področju strojništva.
Njegovo zgodnje raziskovalno delo je potekalo na področju fizikalnih osnov brušenja, iz česar je tudi doktoriral. Razvil je novo metodo merjenja temperature na konici zelo hitro gibajočih se brusilnih zrn in pri analizi eksperimentalnih rezultatov prvi upošteval naključni značaj brusilnega procesa. Kot prvi je uvedel statistično vrednotenje merskih rezultatov. To so temeljna dela, ki jih še danes prebirajo in citirajo številni raziskovalci s področja brušenja. Delo na področju brušenja je zaokrožil z izvirno metodo za označevanje trdote brusov, ki jo je tudi patentiral.
Raziskovalno delo je nadaljeval na področju avtomatizacije obdelovalnih sistemov. Ko je služboval v Birminghamu, je sodeloval na številnih projektih z vrhunsko angleško industrijo. Med drugim tudi na projektu MOLINS 24, v okviru katerega je bil razvit prvi računalniško krmiljen fleksibilni obdelovalni sistem (FMS) na svetu. Že pri raziskavah brušenja je spoznal, da deterministični modeli obdelovalnih procesov niso ustrezni. Zato je pri svojih raziskavah izhajal iz naključnega značaja procesov in ustrezne statistične obdelave eksperimentalnih podatkov, s čimer je postavil nova izhodišča za raziskovalno delo na področjih modernih proizvodnih tehnologij. Iz tega obdobja izhajajo najvidnejši uspehi njegovega dela
for good to found, as a full professor, the Department of Technical Cybernetics, Manufacturing Systems, and Computer Technology, along with the Laboratory of the same name. Peklenik introduced new, modern subjects, and a modular method of studying. He designed the Mechatronics study track and a postgraduate school for automation and production cybernetics. Under his mentorship more than 100 undergraduates received their degrees, over 50 their master's, and Peklenik oversaw 30 doctoral dissertations.
During his time at the Faculty, Peklenik worked in a number of leadership positions, serving as Department Head for many years. From 1974 to 1976 he was dean of the Faculty of Mechanical Engineering, and University of Ljubljana rector from 1987 to 1990.
After his retirement in 1996, Peklenik still worked with the laboratory and was busy on several other fronts as well. Over time he reduced his commitments and spent his final years in the company of his friends and family. Professor Peklenik presented the results of his work in a staggering magnum opus, comprised of extraordinary, oft-cited publications and a whole heap of domestic and international patents. His scientific work spans a broad range of discoveries and findings in mechanical engineering.
The first studies he made were on the physical foundations of grinding, which is also what his dissertation explored. He developed a new method of measuring temperature on the tip of rapidly-moving abrasive grains, and in his analysis of the experimental results was the first to consider the random nature of the grinding process. Peklenik was also the first to introduce statistical evaluation of measurement results. These were the ground-breaking studies that are still cited today by researchers and experts on grinding. He then turned his research on grinding into an original method for marking the hardness of grinding wheels, which he also patented.
Peklenik's research continued on the automation of manufacturing systems. While working in Birmingham, he collaborated on several projects at the highest levels of English industry. This included the MOLINS 24 project, which developed the world's first computer-controlled flexible manufacturing system (FMS). Already in his work on grinding, Peklenik figured out that deterministic models of manufacturing processes were inappropriate. This compelled him to devise an approach based on the random nature of processes and suitable statistical analysis of experimental data, which he used to set new fundamental starting points for research in modern production technologies. The clearest successes of his work in this period were in real-time identification and adaptive control of machining processes and systems, and the description
na področjih sprotne identifikacije ter adaptivnega krmiljenja obdelovalnih procesov in sistemov, opisa in karakterizacije naključnih lastnosti tehničnih površin ter površinskih vmesnikov.
Ko se je računalniška tehnologija začela intenzivno razvijati, se je zopet med prvimi zavedal industrijskega pomena uvajanja računalniško podprtih proizvodnih in krmilnih tehnologij in njihove integracije. Pri tem je spoznal, da je potrebno spremeniti Taylorjevo paradigmo proizvajanja in začel je iskati ustrezne rešitve. Bil je med prvimi, ki so se začeli zavedati problema naraščajoče kompleksnosti v proizvodnji in pomena vloge človeka, oziroma kot ga je sam imenoval, subjekta, pri tem. Kot odziv na ta spoznanja je razvil teorijo proizvodne kibernetike, ki omogoča sistematično načrtovanje, realizacijo in operiranje kompleksnih sistemov v katerih so prisotni subjekti, kot je to značilno za proizvodne sisteme.
Za svoje znanstvene dosežke je prejel mnogo visokih mednarodnih odlikovanj, med drugimi Taylorjevo nagrado mednarodne akademije CIRP (1960), Okoshi medaljo na Japonskem (1974), Taylorjevo nagrado SME v ZDA (1980) in Georg Schlesingerjevo nagrado v Nemčiji (1988). Njegova mednarodna uveljavljenost je razpoznavna tudi iz številnih predsedovanj na znanstvenih konferencah in seminarjih ter vabljenih predavanj ki jih je imel na vrsti univerz po celem svetu. Za njegov prispevek sta mu Univerza v Birminghamu v Veliki Britaniji (1972) in Nankinški aeronavtični inštitut na Kitajskem (1982) podelila naslov častnega profesorja.
Od slovenskih priznanj je potrebno omeniti Kidričevo nagrado (1974), naziv Ambasador znanosti Republike Slovenije (1992), Nagrado Republike Slovenije za znanstveno-raziskovalno delo (1996), naziv Zaslužni profesor Univerze v Ljubljani (1996), ter naziv Častni občan občine Kranj (2012).
Na osnovi rezultatov svojega dela je bil izvoljen v več domačih in mednarodnih akademij. Bil je redni član Slovenske akademije znanosti in umetnosti, ustanovitelj in častni predsednik Slovenske inženirske akademije, redni član Mednarodne akademije za proizvodno inženirstvo CIRP, Ruske inženirske akademije ter akademije Academia Europaea.
Profesor Peklenik je bil vrhunski znanstvenik svetovnega merila, ki je v svoji karieri pomembno prispeval k razvoju znanosti na področju strojništva. Z njimi smo izgubili velikega učitelja in vzornika. Bil je nesporno véliki človek. Za njim pa ostaja njegovo delo in ponosni smo, da lahko to delo tudi uspešno nadaljujemo. Ohranili ga bomo v globokem spominu.
Prof. dr. Peter Butala
and characterization of random properties of technical surfaces and surface interfaces.
When computer technology began developing intensively, he was again among the first to realize the importance of industrial computer-aided production and control technologies and their integration. He discovered in doing so that it was necessary to change Taylor's production paradigm, and he began searching for good solutions. Once again he was among the first to realize the problems associated with the ever more complex nature of production and the importance of the role of humans, or "subjects", as he called them. As a reaction to these findings, he developed the theory of production cybernetics, which makes possible systematic planning, implementation and operation of complex systems in which human "subjects" are present, as is characteristic of production systems.
His scientific achievements earned him the F.W. Taylor Medal of CIRP (1960), the Japanese Okoshi Medal (1974), the Taylor Research Medal in the USA (1980), and the Georg Schlesinger award in Germany (1988). His international recognition was evident from the several positions he held at professional conferences and seminars, as well as the countless lectures he was invited to give throughout the world. For his contribution he was appointed honorary professor of University Birmingham, United Kingdom (1972), and honorary professor of Nanjing University of Aeronautics and Astronautics, China, (1982).
Among his Slovenian awards and recognitions, the Kidrič Prize in 1974 stands out, as well as his title as Ambassador of Science of the Republic of Slovenia (1992), the Republic of Slovenia award for scientific research (1996), the title of Professor Emeritus of the University of Ljubljana (1996), and the title Honorary Citizen of Kranj (2012).
The results of his work earned him elections to numerous domestic and international academies. He was a full member of the Slovenian Academy of Sciences and Arts, founder and honorary member of the Slovenian Academy of Engineering, a full member of the International Academy for Production Engineering, of the Russian Academy of Engineering, and the Academia Europaea.
Professor Peklenik was a world-class scientist, whose career was a vital contribution to the development of knowledge in mechanical engineering. With his passing we have lost a great teacher and role model. He was without question an extraordinary man, who accomplished a great deal of work, and we are proud to be continuing in his footsteps. We will keep his memory deep in our hearts.
Professor Peter Butala, Ph.D.
Information 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 10 pages. Longer contributions will only be accepted if authors provide juastification 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.
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 vestnik - Journal of Mechanical Engineering, vol. 57, no. 7-8, p. 531-538, D0I: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 D e sorption 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). 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 240.00 EUR (for articles with maximum of 6 pages), 300.00 EUR (for articles with maximum of 10 pages), plus 30.00 EUR for each additional page. The additional cost for a color page is 90.00 EUR. 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
Hubert KosLer, Urban Pavlovčič, Matija Jezeršek, Janez Možina:
Adaptive Robotic Deburring of Die-Cast Parts with Position
and Orientation Measurements Using a 3D Laser-Triangulation Sensor
Shijun Ji, Huijuan Yu, Ji Zhao, Xiaolong Liu, Mingxu Zhao: Ultra-Precision Machining of a Large Amplitude Sinusoidal Ring Surface Based on a Slow Tool Servo
Andrej Škrlec, Jernej Klemenc:
Estimating the Strain-Rate-Dependent Parameters of the Cowper-Symonds and Johnson-Cook Material Models using Taguchi Arrays
Dan Ni, Minguan Yang, Bo Gao, Ning Zhang, Zhong Li:
Flow Unsteadiness and Pressure Pulsations in a Nuclear Reactor
Coolant Pump
Deniz Duran, Celalettin Karadogan:
Determination of Coulomb's Friction Coefficient Directly
from Cylinder Compression Tests
Ahmet Selim Dalkilic, Ali Celen, Alican Cebi, Somchai Wongwises: Effect of Refrigerant Type and Insulation Thickness on Refrigeration Systems of Land and Sea Vehicles
252
9770039248001