JOURNAL OF
MECHANICAL ENGINEERING
STROJNISKIVESTNIK
year 2008 volume 54
Contents
Strojniški vestnik - Journal of Mechanical Engineering volume 54, (2008), number 2 Ljubljana, February 2008 ISSN 0039-2480
Published monthly
Editorial
Alujevič, A.	80
Papers
Rihtaršič, J., Šubelj, M., Hočevar, M., Duhovnik, J.: Flow Analysis Through the Centrifugal Impeller
of a Vacuum Cleaner Unit	81
Čas, J., Klobučar, R., Hercog, D., Šafarič, R.: Teleoperation of SCARA with Neural Network Based
Controller	94
Imrak, C.E.: Artificial Neural Networks Application in Duplex/Triplex Elevator Group Control System 103 Lemeš, S., Zaimović-Uzunović, N.: Using Buckling Analysis to Predict Wrinkling in Incremental
Sheet Metal Forming	115
Dašić, P., Natsis, A., Petropoulos, G.: Models of Reliability for Cutting Tools: Examples in
Manufacturing and Agricultural Engineering	122
Ljubenkov, B., Đukić, G., Kuzmanić, M.: Simulation Methods in Shipbuilding Process Design	131
Bakhshandeh, K., Rajabi, I., Rahimi, F.: Investigation of Stress Concentration Factor for Finite-Width Orthotopic Rectangular Plates with a Circular Opening Using Three-Dimensional Finite Element Model	140
Bezgovšek, J., Grabec, I., Mužič, P., Govekar, E.: Development of a Snake-like Robot	148
Corrigendum	154
Instructions for Authors	155
Editorial
In this issue of our Journal of Mechanical Engineering you will find eight scientific papers, one of them being accepted recently by the Publishing council, and the remainder by the old Editorial board in 2007. At present we have a large quantity of papers already accepted and/or being refereed abroad. If we continue to publish at the same rate of eight papers per issue, within next 10 issues, this means that papers submitted during 2008 can only be accepted for 55th volume. Bearing in mind that next year is exceptional, "once on the blue moon" (December 31, 2009) when I do wish to retire from my active job as University Professor, it can simply mean that my new deputy editor Professor Vincenc Butala is going to take over as soon as reasonably achievable.
This year 2008 we also commemorate the 100 years of the Eminent Professor Bojan Kraut, who was the founding Editor of our Journal of Mechanical Engineering, back in 1955, being also
well known by his pocket Mechanical Engineering Handbook (fist issued 1954 by TZ Litostroj). Professor B. Kraut was born in Kamnik, on April 12, 1908, and graduated in 1932 at TF in Zagreb. From 1936 to 1946 he has been employed at the Railway wagons factory in Slavonski Brod. In 1948 to 1950 he has been Technical director of Litostroj, and also part time Professor at the Department of Mechanical Engineering in Ljubljana, and at the newly formed Faculty of Mechanical Engineering in 1962 he became full Professor for Technology of Metals in Railways till his retirement in l972. In 1978, he has been awarded the title of an Eminent Professor by the University of Ljubljana, and Dr.h.c. by the University of Maribor in 1984. He died on August 22, 1991 in the age of 83 years, and Journal of Mechanical Engineering is one of his legacies.
Editor
Prof. Dr. Andro Alujevič
Prof. Bojan Kraut
Paper received: 3.1.2008 Paper accepted: 12.3.2008
Flow Analysis Through the Centrifugal Impeller of a Vacuum
Cleaner Unit
Janez Rihtaršič* - Matjaž Subelj - Marko Hočevar - Jože Duhovnik University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
Flow through the centrifugal impeller of a vacuum cleaner unit was investigated. This was done using experimental measurements and numerical simulations. Two different operational points were selected for the experiment.
Experiments were performed using flow visualization with the use of the passive seeding particles. The rotational speed was reduced by using the theory of dynamical similarity where water was used as a flow medium instead ofair. Numerical simulations were performed by using compressible and incompressible fluids. Flow through impeller channels is presented as particle trajectories and relative velocities for experimental approach and as flow streamlines and relative velocities for numerical approach. © 2008 Journal of Mechanical Engineering. All rights reserved.
Keywords: flow visualization, centrifugal impeller, dynamical similarity, numerical simulations
0 INTRODUCTION
Aerodynamic geometry of vacuum cleaner units is a result of compromise between diverse working conditions; from high volume flow to minimum volume flow caused by full dust bag or clogged filters. Wrong geometry results in unsatisfactory performance, increased noise and surge or stall effects. To achieve a wide area of optimum operational conditions of the vacuum cleaner, it is necessary to know and understand the properties of the fluid flow at all operation conditions. This is possible to achieve through experimental measurements and/or with numerical analysis.
Numerical simulations of flow are often and successfully used for the flow visualization and flow analysis, meanwhile experimental visualizations are less frequently conducted. However, both approaches supplement each other and they are indispensably needed to understand some of the phenomena. In these cases, experimentally gathered data are used as boundary conditions for numerical simulations and sometimes, especially outside of the optimum working conditions, experimental data are sometimes the only source of information.
There are many different methods of experimental flow visualization. Some use smoke or coloured fluid other use solid additives such as small particles [1] to [5].
Two different methods of experimental flow visualization with the use of solid particles are presented in references [1] and [2]. Experimental visualization of flow through centrifugal impeller was described by Pedersen and Jacobsen [1], where Particle Image Velocimetry (PIV) method was used. They performed measurements at nominal flow and at one fourth of the nominal flow. They noticed that at operation in nominal conditions, flow follows the blade curvature without separations, whereas at operation at one fourth of the nominal flow, separations of the flow from the suction side of the blades were observed. Another approach to experimentally visualize the flow through the impeller is by capturing trajectories of the seeding particles. Such method was applied in Irabu et al. [2], where polystyrene particles were used to visualize the flow through the volute. The lather method was also used for visualization of flow through impeller at two different working conditions.
1 DYNAMICAL SIMILARITY OF THE FLOW
Application of the visualization technique to capture trajectories of the individual particles dictated the reduction of flow velocities in the centrifugal rotor. On the basis of the theory of flow similarity, this was achieved using water as a flow medium. In the follow-up, we will discuss the
*Corr. Author's Address: University of Ljubljana, Faculty of Mechanical Engineering, Aškerčeva 6, SI-1000 Ljubljana, Slovenia, janez.rihtarsic@lecad.uni-lj.si	81
geometrical and dynamical similarity and the selection of the measurement point.
Similarity of flows is achieved by assuring the geometrical similarity where we assume that the losses through the impeller are not dependent upon the magnitude of the flow (one impeller is a true geometrical scale of the other) and with assurance of the dynamical similarity [6] and [7]. To fulfill the requirements of the dynamical similarity, test conditions should be chosen to satisfy the Reynolds, Euler, Thoma, Froude and Weber similarity numbers, however this is usually impossible. Thoma, Froude and Weber similarity numbers deal with cavitation. Here, such conditions were selected that no cavitation occurred in all operational points and thus these numbers were not considered. The Reynolds and Euler similarity numbers are explained later.
Flow through two geometrically similar impellers is equal if the velocity fields of both flows are identical. We used the identical impeller for measurements as it is used in the vacuum cleaner unit (Fig. 2), thus the geometrical similarity was assured. By fulfilling dynamical similarity and with similarity of the Reynolds number, we can expect development of similar boundary layers and separations of the flow, as they appear with the use of air as a flow medium. Compressibility conditions were not satisfied but their influence on stream lines and relative velocities was estimated to be low. This was found by comparing CFD simulations where compressible and incompressible flows were compared (chapters 4.1.2 and 4.2.2).
When operating in the same flow medium at two different rotational speeds, the ratio of the volume flows is proportional to the ratio of the rotational speeds and the ratio of the static pressures is proportional to the square ratio of the rotational speeds. Ratio of the static pressures is also proportional to the ratio of the densities of the flow medium if the rotational speed is constant [8]:
Qa2
Pa2 = Pal
Eal n
^ Qal =-
' Qa2
2
Pa2- . n 2
Pal
Pa2 _Pa . p _Pw p
- _--> Pw _--pa2
Pw Pw	Pa
(1) (2) (3)
rotational speed, static pressure, flow medium density, air flow medium index, water flow medium index.
From Equations (1) to (3), we get dependence between rotational speed, when the flow medium is air, and rotational speed, when the flow medium is water (Eq. 4):
n	[min-1]
p	[Pa]
p	[kg/m3
a	
w	
PW .Е
Pal Pw
(4).
Because we used identical impeller for both flow medium, the change in Reynolds number, depends only upon the change of rotational speed and upon the change of the kinematic viscosity. Relation between rotational speeds and viscosities is given in:
nw
n
(5)
V [mm2/s] kinematic viscosity.
Using water as a working medium, we were able to adjust Reynolds number by changing the water temperature. However, due to the selection of the test rig materials and due to capturing frequency of the CCD camera, we were not able to achieve the same Reynolds number as when using air as a flow medium. We heated the water to 60° C. This way, we were able to obtain similarity of Reynolds number of around 31% for both volume flows.
Euler number was fulfilled because the measurement conditions were selected in accordance to the Equation (4).
2 CFD MODEL
Numerical analysis was conducted with commercial code CFX 5.7 by Ansys Inc. It is a 3-dimensional unstructured mesh code that solves a set of Reynolds-averaged Navier-Stokes equations (RANS) [9]. In the RANS approach, the mass and the momentum conservation equations together with the equations of the turbulence model form a closed set of equations. The mass conservation equation reads:
dp d
dt + dxj^0
(6).
2
nw _ na2 _
n
V w2 _ V al
V al V
n
n
n
^ na2 _
n
Here p is the density of the fluid and u is velocity of the fluid. The momentum conservation equation is:
d , 4 d / \ dP d
— (put )+-(pu.u, ) =--+-
dt	dx,	dx, dx.
d
dP д
Mf
dut
эХТ
Эх,
(7),
where P denotes modified pressure and /л denotes viscosity.
Лeff =Л+Л,
The closure of the set of the equations is simplified to the determination of the turbulent viscosity л, which is usually determined by introducing the equations of the turbulence model. For this case, a two equation k-e turbulence model was applied [9] and [10]. Additional equations read:
d(pkuj )
d(pk )
dt
дх,
_d_ dx,,
/ \ Лt
л + —
dk
dx,
+ Pk-Pe (8)
and d(pe) dt
d(e
dx,
_d_ dx ,
Mi
Л +
Oe
de
dx,
+ -(Ce!Pk - Ce2pe) k
(9),
where k is the turbulent kinetic energy, e is the turbulent kinetic energy dissipation rate, Pk represents the generation of turbulent kinetic energy due to the mean velocity gradients, ok, O-,
C, and C, are model coefficients. The numerical
-1 -2
model uses a finite volumes, coupled implicit, pressure based solution technique. A highresolution differentiating scheme was used, which locally adjusts the discretisation to be as close to Second-Order as possible, while ensuring boundedness of the solution.
2.1. Simulation
In the numerical study of the vacuum cleaner unit, a steady state approach was used with a single rotating frame of reference. A stage interface was used between rotating impeller and both steady inlet region and diffuser. Stage interface is modeled using a mixing plane - flow quantities are averaged in the pitch wise direction through the mixing plane, whereas their actual distribution is maintained in axial direction. An unstructured tetrahedral mesh was used with hexahedral layer mesh at the wall boundaries. The analysis was carried out through the whole suction unit, including diffuser and motor part (Fig. 1). The overall mesh consisted of477,149 nodes; one impeller blade passage was modeled with 30,070 nodes. A periodic boundary condition was used to model full geometry. A mesh sensitivity study was carried out to ensure that appropriate mesh density was chosen for the simulations. While finer mesh (640,178 nodes) did not show any considerable improvement in accuracy of flow quantities, coarser mesh (375,983 nodes) exceeded the predicted inlet pressure and air efficiency of vacuum cleaner unit by as much as 20%.
Fig. 1. Mesh model of the impeller, diffuser, return channels and the electrical motor.
Outlet region is not shown here.
The walls were modeled as adiabatic smooth walls and the scalable wall functions were used as a boundary condition. Measured values of mass flow rate and static temperature were prescribed at the inlet, static pressure and temperature at the outlet. Outlet boundary condition was prescribed as an opening at the exhaust openings of motor housing.
Two sets of simulations were carried out using air as a compressible and incompressible working fluids. In the compressible case, we used air as an ideal gas and in the incompressible case we used air with constant properties at 25 °C. Of the whole volume flow / pressure characteristics curve of the suction unit, we analyzed two operating points that match the measurement points of the experiment. Numerical results for both selected volume flows show very good agreement with the pressure/volume flow characteristics of the suction unit with discrepancy of inlet pressure and aerodynamic efficiency values being within 4%. Residual target of 10-4 was achieved in converging calculations, usually within 300 iterations.
3 EXPERIMENTAL SETUP 3.1 Test Rig
Measurements were conducted on a shrouded impeller with nine 2D blades (Fig. 2). The cover of the vacuum cleaner unit and the shroud of the impeller were manufactured from perspex, which enabled us to follow the particles through the entire impeller channel.
Impeller was submerged in the test rig as shown in Figure 3. A glass tank was used with dimensions 400x400x350 mm. The flow medium
was circulating inside the test rig so that working conditions were constant. Testing impeller was operating inside the original vacuum cleaner unit, driven by external asynchronous motor. Asynchronous motor was controlled by a frequency regulator.
Seeding particles with radius of 0.5mm were made out of polypropylene to have the same density as working medium. Particles were introduced into the impeller at radius R0=0 mm. For visualization we used non-interlaced CCD camera Sony XC HR-50. Resolution of the camera was 640x480 pixels, which equals 140x105 mm. Sampling frequency is f =60 Hz and exposure time is t = 510-4 s. Camera was positioned under a glass tank in the axis of the impeller. For illumination we used ring illuminator Vega Velum DC. With appropriate positioning of the illumination source and the camera, the reflections of light were reduced to a minimum. The images were digitized and recorded using a frame grabber NI 1409 and a NI Vision software module running on a Labview software platform.
Temperature was measured using Pt-100 resistance thermometer and maintained in the liquid vessel at the temperature region ±1°C.
3.2 Selection of Measurement Point
Figure 4 shows pressure/volume flow characteristics curve of the vacuum cleaner suction unit under operation with air. Measurements with water were conducted in dynamically equivalent points at high volume flow Q1=1.46Qn, n =336.8 rpm (measurement point 1) and nominal flow ß2=ßn, n2 =359.5 rpm (measurement point 2).
Measurement points 1 and 2 were obtained by setting the rotational speeds which were calculated
R2
127°
a)	b)
Fig. 2. Geometry of the centrifugal impeller; a) impeller 3D cross section and b) impeller blade angles
in accordance with Equation 4 and by setting appropriate water column height that corresponded to selected static pressure. Both measurement points are on the same pressure/volume flow characteristics.
The sampling frequency of the camera and the dimensions of the test rig were considered when selecting the pressure/volume flow characteristics. Thus, we were able to record sufficient number of
CCD camera (SONY XC HR-50)
Fig. 3. Closed loop test rig
consecutive images of the same particle at high volume flow, and the height of the water column was kept moderate at high pressures (nominal volume flow).
In each of the two measurement points, we conducted 15 measurements and for each of the measurement we recorded a series of 200 consecutive images.
35 30 25
77 20
Q.
a. 15
10
5
				
				
			2	
				
				
				\ 1
				
10
20
30
40
50
Q [l/s]
Fig. 4. Pressure/volume flow characteristics curve of the vacuum cleaner suction unit under operation with air. Dots show measurement point at a high volume flow Q=1.46Q (measurement point 1) and at nominal flow Q2=Qn (measurement point 2).
3.3 Image Analysis
The purpose of image analysis was to establish occurrences of particles in the centrifugal rotor.
For the purposes of the analysis, a computer program called TRACER was written in C++.
Tracer enabled to rotate pictures in correlation with impeller revolutions. This enabled to follow the particle through the seemingly static impeller blade channels. Furthermore, this feature also enabled to make a computing corrections of the impeller frequency, which was varying during the measurements within the range of 1 Hz.
Fig. 5. Trajectories of a single measurement at nominal volume flow
Fig. 6. Trajectories of a single measurement at nominal flow after the rotation into a single blade
channel
Measurement error due to changing frequency was already reported in reference [1], where its influence was much greater because of the longer measuring time.
On average, we captured 16 trajectories per measurement and the measurement itself lasted 3.3 secs. Trajectories were captured through different impeller blade channels (Fig. 5). Occurrences of particles were defined with accuracy of 1 pixel which equals 0.22 mm. Trajectories were further rotated into a single impeller blade channel (Fig. 6).
4 RESULTS OF MEASUREMENTS AND COMPARISON WITH CFD SIMULATIONS
Results of experimental measurements are based on image analysis, performed with the TRACER program. From the established occurrences of the particles we were able to define trajectories of a single particle and its relative velocity. Altogether, 246 trajectories were captured at measurement point 1 and 199 trajectories were captured at measurement point 2.
4.1 Trajectories
4.1.1 Experimental Results
At each flow rate, 15 measurements were conducted and for each measurement, an average trajectory was calculated. In this way, we got 15 average trajectories in each of the two measurement points. To calculate average trajectory, we divided the length of the channel into several sections. The
measurements at Q1=1.46Qn in n1 =336.8 rpm. Dashed line indicates the middle line of the blade channel.
number of the sections equals the lowest number of detected occurrences of particles of individual trajectory in the measurement. If one trajectory had several occurrences in one section an average location of this particle in the observed section was calculated. In this way, particles with lower relative velocity do not have greater influence on average trajectory than particles with higher relative velocity, which have consecutively lower number of occurrences of particles.
Envelopes of average trajectories are shown in Figures 7 and 8. Hatched areas indicate region of occurrence of average trajectories.
It can be seen that average trajectories move from the pressure side towards the suction side of the impeller blade channel when moving from high flow at measurement point 1 (Fig. 7) towards the nominal volume flow at measurement point 2 (Fig. 8).
Comparison between the average trajectory of particles in each of the measuring points and the blade curvature is shown in Figure 9. Tangential angles of average trajectories and of the impeller blades are given in table 1. Angles were measured at radii 22, 28, 34 and 40 mm. The relative flow angles were not determined at inlet and at outlet of the channel. The reason for this is the method of measurement, which enables determining trajectories from occurrences of particles on visualization images. Because the acquisition frequency of the camera was limited to 60 Hz, images with occurrences of particles at the inlet and outlet were seldom available. From Figure 9 we can notice congruence between the average trajectory at nominal flow and the impeller blade. Average trajectory at large flow is more back swept
measurements at Q2=Qa in n2 =359.5 rpm. Dashed line indicates the middle line of the blade channel.
Table 1. Tangential angles of the average trajectories and the impeller blade at diferent radii
	ß [°] R=22mm	ß [°] R=28mm	ß [°] R=34mm	ß [°] R=40mm
Impeller blade	33.4	34.7	29.9	28.4
Average trajectory in measurement point 1	25.4	34.3	27.7	28.6
Average trajectory in measurement point 2	35.5	31.1	32.4	27.8
Fig. 9. Average trajectory of particles at high flow Q1=1.46Qn (dashed line) and at nominal flow Q=
(continuous line)

in comparison to the impeller blades. Average trajectories also indicate negative incidence angle at nominal flow and positive incidence angle at high flow.
At both measurement points, we observed trajectories that indicate existence of turbulent flow. At nominal flow in measurement point 2, they represented 7% and at high flow at measurement point 1 they represented 3 % of all recorded trajectories. The reason for lower number of turbulent trajectories at high flow in comparison with nominal flow can also lie in the fact that we
have lower number of captured occurrences of particles at high flow and thus less detailed trajectory can be drawn.
Figures 10 and 11 show four characteristic turbulent trajectories that occurred in each of the two measurement points. As we observe from Figures 10 and 11, trajectories at large and nominal flow detach from the pressure side of the blade, which indicates vortices that appear on the last two thirds of the blade. Observations also indicate that these vortices are not steady but they are appearing and disappearing with time.
Fig. 10. Trajectories that indicate turbulence at high flow Ql=1.46Qn
Fig. 11. Trajectories that indicate turbulence at nominal flow Q2=Qn
4.1.2 CFD Analysis
Numerically calculated streamlines for high and nominal flows are given in Figures 12 and 13. Numerical results when using compressible fluid (Figs. 12a and Fig. 13a) and when using incompressible fluid (Figs. 12b and 13b) show negligible differences when comparing streamline plots. However, there is an observable difference in the overall air efficiency of a vacuum cleaner
unit. Difference in efficiency between measured and numerical results for compressible fluid was 3.72 % at nominal flow and 4.45% at high flow, while the difference in efficiency between measured and numerical results for incompressible fluid was 11.12% at nominal flow and 23.00% at high flow.
We can observe at experimental evaluation as well as at numerical simulation that there is detached fluid on the pressure side of the blade
a)	b)
Fig. 12. Numerically calculated streamlines at high flow Q1=1.46Q using a) air as a compressible working fluid and b) air as an incompressible working fluid
a)	b)
Fig. 13. Numerically calculated streamlines at nominal flow Q2=Qa using a) air as a compressible working fluid and b) air as an incompressible working fluid
channel. At both operational points experimental results agree with numerical results well.
4.2 Relative Velocity
By following the particles through impeller blade channels, we were able to calculate relative velocities. The path was determined with the distance between two successive occurrences of the particle and the time was determined by sampling frequency. When calculating relative velocities in this manner, one has to assume that velocity between two successive occurrences of particles is constant and that the flow through the impeller is two dimensional.
4.2.1 Experimental Results
From the experimental results, the velocities were estimated from two consecutive occurrences of the particles. The velocity was assigned to the middle point between both consecutive occurrences of particles. Due to this, the density of points with assigned velocity near the inlet and the outlet of the blade passage was low. This resulted in a lower accuracy of the inlet and outlet relative velocities (Figs. 14 and 15). The accuracy depends on the number of velocity points and this was not the same in all parts of the blade channel. The mesh in Figures 14 and 15 was calculated using the 9th order two dimensional polynomial regression. The results were calculated using Equations (1) to (5) from water to air as a medium to enable comparison with CFD results.
Figure 14 shows relative velocities at high volume flow. The relative velocity at the inlet channel is higher on the suction side of the blade. This region extends to approximately one third of the blade length. Further to the outlet region near the suction side of the blade, the relative velocities are low and the peak of relative velocities shifts towards the center of the channel and to the pressure side. There are areas of small relative velocities near the pressure side of the blade channel and at the tip of the pressure side of the blade at the outlet. However, this could be due to the low number of data used for meshing in this region.
Figure 15 shows relative velocities at nominal volume flow. The most important feature is the presence of a region with low relative velocities at the suction side of the blade which extends from approximately one third of the blade to the outlet of the blade. At the inlet, the measured relative velocities are higher, although this could also be due to the low number of velocity points in this region. At the outlet of the channel, there are higher velocities on the pressure side of the blade. There is a region of low relative velocities in the middle of the pressure side, similarly to what happened to a smaller extent in the case of high volume flow (Fig. 14).
4.2.2 CFD Analysis
The relative velocities were calculated. Results are shown in Figures 16 and 17. Figure 16 shows relative velocities at a high volume flow. On the
[m s*-1]
Velocity		
		95
		85
		75
		65
		54
	_	44
J	■	34
J	■	24
	■	14
	1	4
[m		5л-1]
Fig. 14. Relative velocities at high flow Ql=1.46Qn
Fig. 15. Relative velocities at nominal flow
Q=Q„
pressure side of the blade there is a region of secondary flow with very low relative velocities, starting from the tip of the blade and extending to approximately two thirds of the blade length. The peak of relative velocities amplitude is located at approximately one fourth of the channel close to the suction side of the blade. The relative flow
velocity at the outlet of the channel is low on its suction side and high on its pressure side. The maximum velocities exceed 100 m/s.
Figure 17 shows relative velocities at nominal volume flow. Here, the point of the highest velocity is close to the tip of the blade on its suction side. The position of the highest relative velocities
[m s*-1]
[m s*-1]
a)	b)
Fig. 16. Numerically calculated relative velocities at high How Q1=1.46Qn using a) air as a compressible working fluid and b) air as an incompressible working fluid. Plots show relative velocities
at middle surface between the hub and the shroud.
Fig. 17. Numerically calculated relative velocities at nominal flow Q2=Qn using a) air as a compressible working fluid and b) air as an incompressible working fluid. Plots show relative velocities at middle
surface between the hub and the shroud.
has moved from the suction side to the center of the channel. The same phenomenon of the secondary flow at the pressure side of the blade and distribution of velocities at the outlet of the channel was established as in Figure 16. Maximum velocities exceed 85 m/s.
No significant differences in relative velocities amplitudes and distribution were established between compressible and incompressible flows. This supports the use of experimental technique based on incompressible water flow visualization.
Comparison between experimental results and numerical simulation shows good agreement between both methods. Using both methods, we observe a region of high relative velocities along the suction side of the blade located at approximately the same position. Both methods similarly predict a region of low relative velocities at the end of the blade at the suction side. We also notice a region of low relative velocities along the center of the pressure side of the blade and a region of high relative velocities at the outlet of the channel at the pressure side of the blade.
5 DISCUSSION
In the follow-up, we will discuss the above results and compare results of experimental work with results of CFD analysis.
Particle size and material selection is influenced by inertia, gravity and buoyancy forces of the particles. To ensure following of the particles to the flow, the selected density of particles matched the density of flow medium. The diameter of particles was dictated by the illumination technique. We used continuous illumination and thus, a relatively large size (radius r = 0.5mm) of particles was necessary, in contrast to the size of particles used in PIV methods [1]. Using the theory of similarity, influence of the inertia of particles was reduced by lowering the number of rotations.
The compressibility of air should be taken into account generally as Mach number increases to 0.3 or higher. In the current research, we achieved maximum velocities of up to Mach 0.6. By using water as a medium, we were not able to consider compressibility and we relied on the comparison between CFD results with compressible and non compressible flows. As shown in Figures 12, 13, 16 and 17, differences between both cases exist, but
they are generally negligible. Although there was no significant difference between streamlines and relative velocities, there is a noticeable difference between the suction efficiencies of the two flows.
The flow is also influenced by the Reynolds number. In our case, the Reynolds number of the measurements in water was approximately 31% of the corresponding Reynolds number in air.
The relative velocities as a result of CFD analysis are shown in Figures 16 and 17. At both operational flows, we noticed that average relative velocities calculated using CFD method are on the average a little higher than the corresponding measured velocities shown in Figures 14 and 15. In our opinion, the main reason for this is the selection of the plane where relative velocities are shown. This was at the center of the channel for the CFD, while in case of experimental techniques; velocity of particles was averaged over the entire width of the impeller flow channel.
A comparison with CFD showed that the presented experimental method enables measurements of flow properties inside the impeller with a reasonable accuracy.
These results indicate further possibilities for the changes of the geometry of the impeller to achieve better flow conditions.
6 CONCLUSION
Numerical simulation and flow visualization through the vacuum cleaner centrifugal impeller at nominal (Qn) and at high flow (1.46Qn) were compared.
Flow visualization was conducted with the use of the theory of dynamical similarity where experimental flow medium was water. The experiment confirmed good matching between the calculated and measured values for pressure difference, rotational speed and volume flow, when using air as a flow medium and when using water as a flow medium.
Both experimental visualization and numerical simulation indicate detachment of the flow after one third of the blade length from the pressure side of the blade. Detachment occurs at high flow as well as at nominal flow. Average trajectory of the particle shows congruence with blade curvature at nominal flow, meanwhile the average trajectory at high flow is more back swept at trailing edge and has a positive incidence angle
at the leading edge. There is a good agreement in distribution of relative velocities between the pressure and the suction side of the channel obtained by experimental measurement and numerical simulation at a high flow. The similarity decreases when comparing distribution of relative velocities at nominal flow.
An advantage of such approach to experimental flow visualization is that it enables observing the flow in different operation conditions. Because it is possible to get numerical solutions for almost all operation conditions, it is necessarily to get also a confirmation of these results experimentally and thus gain the confidence in numerical computation. Although we achieved relatively good comparison of the results between experimental measurement and numerical simulations, one has to be aware that it is most often impossible to fulfill all conditions for dynamic similarity. The influence of conditions that are not fulfilled should be evaluated and considered when evaluating the final results.
There is also possibility of automation of the presented experimental method by capturing of the positions of the particles in the flow, which would give the results in real time. The practical value of such method would be recognized in cases where there are numerous changing flow conditions.
7 REFERENCES
[ 1 ] Pedersen, N., Jacobsen, C.B. PIV Investigation of the Internal Flow Structure in a Centrifugal Pump Impeller. 10th International Symposium on Applications of Laser Techniques to Fluid Mechanics, Lisbon, 2000.
[2]	Irabu, K., Tamazato, E., Teruya, I. Velocity measurements of flow around a volute tongue of vaneless radial diffuser. Proc. 3rd Asian Symposium on Visualization, Japan, 1994.
[3]	Širok, B., Grabec, I. Chaotis properties of a pulsating two-phase solid particles/gas flow. Strojniški vestnik - Journal of Mechanical Engineering, 34 (1988), p. 162-165. (In Slovenian).
[4]	Novak, M., Širok, B., Hočevar, M., Oberdank, K. Simultaneous measurement of velocity fluctations and smoke tracer fluctations of turbulent flow over a back-facing step. Journal of Flow Visualization and Image Processing, vol 6., 1999, p. 221-229.
[5]	Dular, M., Bachert, R., Širok, B., Stoffel, B. Transient simulation, visualisation and PIV-LIF measurements of the cavitation on different hydrofoil configurations. Strojniški vestnik - Journal of Mechanical Engineering, 51 (2005), p. 132-145.
[6]	Lazarkiewicz, S., Troskolanski, T., A. Impeller pumps. Pergamon Press, Warszawa: Wydawnictwa naukowo-techniczne, 1965.
[7]	ISO 60193, Hydraulic turbines, storage pumps and pump turbines - Model acceptance tests, IEC 60193 (1999-11).
[8]	Bleier, P.F. Fan handbook, selection, application and design. New York: McGraw - Hill, 1998.
[9]	Ferziger J. H., Perić M. Computational Methods for Fluid Dynamics, 2nd ed. Springer Verlag, 1999.
[10]	Davidson, P.A. Turbulence, an introduction for scientists and engineers. New York: Oxford University Press, 2004.
Paper received: 20.3.2007 Paper accepted: 19.12.2007
Teleoperation of SCARA with Neural Network Based
Controller
Jure Cas* - Rok Klobučar - Darko Hercog - Riko Safarič University of Maribor, Faculty of Electrical Engineering and Computer Science, Slovenia
This paper describes the development of neural network based controller for the teleoperation of Selective Compliance Assembly Robot Arm (SCARA). The SCARA is controlled and teleoperated via the internet. Presented experiment is used by students at the University of Maribor as a remote educational tool. Application is based on MATLAB/Simulink and LabVIEW software packages. MATLAB/Simulink and developed library DSP-2 Library for Simulink are used for neural network control algorithm development, simulation and code generation. The executable code is downloaded to the Digital Signal Processor (DSP). The DSP controls through the analog and digital I/O the real process and maintain the data connection with the laboratory server The LabVIEW virtual instrument (VI) is used as a client server application for the teleoperation. LabVIEW VI provides the ability for parameter tuning, signal monitoring, on-line analysis and via Remote Panels technology also teleoperation by using the internet browser (Internet Explorer). The main advantage of a neural network controller is the exploitation of its self-learning capability. For example: when friction or an unexpected disturbance occurs, the user of a remote application does not need any information about the changed robot dynamics, because it is estimated independently of the remote user.
© 2008 Journal of Mechanical Engineering. All rights reserved.
Keywords: robotics, neural network controller, remote experiment system, LabVIEW, MATLAB, Simulink
0 INTRODUCTION
The theoretical development of a continuous neural-network sliding-mode controller (CNNSMC) based on the theory of continuous sliding-mode controllers is presented in the paper. Derived equations of the CNNSMC were verified on the real laboratory teleoperated SCARA mechanism, which is used as remote educational tool. An inverse dynamic model of a robot arm mechanism is needed to design a robot controller but such a modelling of dynamic system is not an easy task. Besides that, the users of the remote robot experiment are unable to physically remove the disturbing causes of a robot's dynamics and, what is even more obstructive, the disturbance dynamics is completely unknown to "the user of the experiment". The consequence is that the conventional procedures for estimating of robot dynamics, which are used in non-teleoperation robot applications, can not be used for the purposes of teleoperating applications.
When using the CNNSMC as a control algorithm however, unknown robot dynamics is not a control problem anymore, because the
approximation of changed robot dynamics is computed independently of the remote user. The additional advantage of application is the possibility for the on-line tuning parameters of CNNSMC, which improves the performance of a robot's controller and makes application more appropriate as an educational tool.
Students are able to learn how the parameters of CNNSMC have an affect on its performance because of the possibility of on-line tuning parameters and observing the responses of regulation using graphs and numerical indicators. They can improve their knowledge of the neural network control algorithm and, what is more important, with the use of internet connections experiments can be executed on a real SCARA from somewhere outside a laboratory .
The developed software, which consists of CNNSMC and graphic user interface (GUI), is being executed on a DSP. DSP is via analog outputs and digital inputs connected with servo-electronics, which is designed for the generation of the required current of DC driving motors and for measuring the position of robot axes with incremental encoders. On the other hand, the DSP-2 robotic
*Corr. Author's Address: University of Maribor, Faculty of Electrical Engineering and Computer Science, Smetanova 94	17, SI-2000 Maribor, Slovenia, jure.cas@uni-mb.si
controller via serial bus and by using library ComVIEW connected with the laboratory server, where the developed VI is running. By using the technology Remote Panels the laboratory server is also responsible for publishing the developed VI on web and for serving the remote client of experiment. Hardware implementation of experiment is shown in Figure 1.
1 SYSTEM COMPONENTS
The DSP-2 Library for Simulink [1] is a MATLAB/Simulink add-on toolbox, which provides rapid control prototyping support for DSP-2 robotic controller. The library contains a set of device driver's blocks for all available I/O ports of the DSP-2 robotic controller, including blocks for analog I/ O, digital I/O, incremental encoder and blocks for communication between the PC and DSP-2 robotic controller (TCP/IP, USB in RS232). Simulink Real Time Workshop and the DSP-2 Library for Simulink enable developers to model applications in the Simulink block-diagram environment and, after successful simulation, quickly verify the designed algorithm with the DSP controller or DSP-2 robotic controller on a real device.
The LabVIEW virtual instrument named ComVIEW has been used for serving fly data visualization and parameter tuning tasks for the DSP-2 robotic controller. When the DSP target is selected in the Simulink model, a LabVTEWvirtual instrument is automatically created from the ComVIEW template VI during binary code generation. A ComVIEW template contains an empty front panel
and a fully functional block diagram. The block diagram implements functions for VI initialization, executable code download to the DSP-2 robotic controller, functions for transmitting and receiving messages between the PC and the DSP-2 robotic controller. During VI creation, numerical controls and indicators are automatically added to the VI front panel template, where the number of controls and indicators depends on the number of DSP communication blocks used in the Simulink model. Links between DSP signals and VI front panel objects are established programmatically using the DSP Connection Manager Window. This window appears on the PC immediately after the downloadable binary code starts executing on the DSP target system. Using mouse clicks, the user can create links between the VI front panel indicators and the DSP output signals, and links between the VI front panel controls and the DSP input signals or DSP parameters. When these links are set, a communication link is established between the VI running on the PC and the code being executed on a DSP controller. Whenever the controls on the VI front panel are changed, LabVIEW automatically downloads them via serial port, to the DSP controller. At the same time, all arrived DSP output signals are read from the PC serial port and displayed on VI control panel as graphs numerical indicators. In the scope mode, a small portion of code running on the DSP controller handles data acquisition and storage management. The selected DSP signals are, firstly, captured and then stored in the temporary controller's memory. After that, the captured data is transferred to the PC.
Fig. 1. Schematic implementation of the developed experiment for teleoperation of the SCARA via internet
Fig. 2. On the DSP-2 robotic controller the developed software is being executed
Remote Panels is a LabVIEWadd-on toolkit developed by National Instruments that enables the viewing and controlling of LabVIEW VI's over the Internet. Using this toolkit, the LabVIEW VI can be published on the internet with no additional programming. Afterwards, the virtual instrument can be remotely observed or controlled by using standard web browser. The remote user can fully access the user interface that appears on the web browser and, consecutively, has complete control of the remote application. Other users can point their web browser to the same URL to view a remote experiment. To avoid confusion, only one client can control the application at a time.
ComVIEW VI and the LabVIEW server are run on the same lab server. ComVIEW VI performs communication between the lab PC and the DSP-2 robotic controller, while the LabVIEW server enables remote operation of the ComVIEW VI. User of telerobotics application must have a LabVIEW Run-Time Engine installed on home computer in order to perform remote experiments.
The DSP-2 robotic controller [2] (Fig. 2) is composed of a DSP-2 controller and a DSP-2 addon robotic board. The key components of the DSP-2 controller are the floating point digital signal processor, used for control algorithm execution, and the Xilinx FPGA, which implements peripheral interfaces. The DSP-2 robotic controller contains all the necessary peripheral for 4-axes robot control i.e. this system has 16 digital inputs, 8 digital outputs, 4 analog inputs/outputs, and 4 incremental encoder interfaces.
3 DESCRIPTION OF THE SCARA
The SCARA (Fig. 3) has two links. Position of the first link in the joint space is presented as
the first degree of freedom ( вг ), and position of the second link in the joint space is presented as the second degree of freedom ( Q2 ) and Q is the vector of the real positions.
Each link is driven by a DC motor, and a gear-box with a transmission ratio 173/19. By using gear-boxes, any nonlinear inertia influences of the presented robot mechanism are decreased, but not completely. Additional dynamic nonlinearities are brought to the system as friction, which is proportionally large (about 20% of maximum torque) in the used SCARA.
Robot links are driven by a DC motor with nominal torque 28.4 [mNm]. The sinus signal from the incremental encoder is by servo electronics transmitted to 400 pulses per joint rotation.
A dynamic model of the S CARA mechanism, with two degrees of freedom is described by the Lagrange equation of motion, as follows:
t = M(&) в+h(q, в)+F(è)+Tn (i),
where T is a vector of the drive torque on the robot's joints, M is an inertial matrix, h is a torque vector due to the centrifugal forces, centripetal forces, and Coriollis forces, F is a torque vector due to friction forces and Tn stands for a torque vector due to unknown disturbances, Q and Q are vectors of real velocities and accelerations of the robot joints. This well known mathematical note of a robot mechanism dynamics (1) can be expressed with an n-dimensional state-space system of equations with regard to the control value u:
i = f ( x, t) + B( x, t) u + d ( x, t)
(2).
New terms are defined as:
x£»n, u e , B(x, t) = B(x, t) + DB(x, t) (3),
Fig. 3. The SCARA is shown
where d is an unknown disturbance, B is an actual input matrix, B is an estimated input matrix, u is a control vector, x is a state space vector of mechanism, and t stands for time.
4 DEVELOPMENT OF CNNSMC
The main advantages of the sliding mode control are the robustness to parameter uncertainty, loading disturbance, and fast dynamics response. However, these properties are valid on the sliding surface under the conditions
of modelling imprecision absence, external disturbances, and switching time delays. The control law of the sliding mode technique contains a discontinuous component and may, on the sliding surface, excite those high frequency dynamics neglected in modelling, due to high control activity. Practical sliding mode control implementations exhibit high frequency oscillations in the plant output, called chattering. This cause harmful effects such as torque pulsation in robot electrical drive and, consequently, imprecise positional control of a robot tip. Various methods have been proposed to eliminate chattering. The most commonly cited approach is smooth approximation of the switching element by saturation, so that a narrow boundary layer is introduced near the sliding surface. The control law is constituted by two components: a continuous law derived from the plant model using the Lyapunov theory, and a saturation law [3]. The first component is required to slide-down on the sliding surface, while the saturation law controls the unknown parts of the robot dynamics.
A number of researchers have suggested methods for alleviating the chattering effects and improve precision in the sliding mode position control of robot mechanisms. In [4] are developed perturbation estimation schemes; other authors suggested supplementing the control input with a predictive correction term [5]. Recently, soft computing techniques fuzzy logic [6], neural network [7] and genetic algorithm [8] have been used to estimate the robot non-linear dynamics, thus needing only the continuous law derived from the robot dynamics using the Lyapunov theory.
Our approach uses a neural network as an estimator for a part or, even complete, robot dynamics. We decided to use a neural network because of its high convergence speed in robot dynamics estimation, robust sliding-mode control scheme, and as little preliminary knowledge on the estimated mechanism dynamics model as possible. In our case, only nominal (average) values were used for inertia matrix parameters. Differences between actual inertia matrix parameters and nominal inertia matrix parameters represent structured uncertainties. Torque terms due to Coriollis forces and friction forces were neglected and they represent non-structured uncertainties. If we did not have a robust control scheme, the robot
behaviour would be unpredictable during the first few moments of learning.
A well-known mathematical note of robot mechanism dynamics (Eq. 1) is transformed into an я-dimensional state-space system of equations with regard to the control value u (Eq. 2), because the theory of Lyapunov for searching the control law can only be used in the following way.
Our goal is to prove the stability of function for the robot system. This means that after transient time, defined with parameters of the matrix G, the difference between the actual and the desired vector of state, space variables x and xr will equal zero and will be stable for all disturbances. Function a( x, t) = 0 will be stable if the Lyapunov function V>0 and the first Lyapunov time derivative of function is V < 0. According to the following definition:
a( x, t) = G( x(t) - xr (t)) = G( x - xr ) (4), where x and x are the vectors of the desired and
г
actual state space variables and G is the matrix defining the control of system dynamics, we cannot prove the robot system's stability (Eq. 2). Nevertheless, we can look for suitable conditions for control law u, where the robot system will be stable. This is done in the following way.
For the simplest Lyapunov function V, to determine the control law u, the Equation (5) has been selected:
V	= aTa /2	(5), The following equation is derived from(5):
V	= aT S	(6).
Owing to the fact that V is not always less than zero for all x, x, and G, the first desired Lyapunov negative time function derivative is defined as:
V--
-aTDa
(7),
where D is a diagonal matrix with positive diagonal elements. If the Equation (7) and the derivative of Lyapunov's function (Eq. 6) are made equal, the result is:
aT (Da-
ti ) = 0
(8).
The Equation (8) is valid if both, or at least one, of the multiplicators equals zero. Since the first multiplicator, the term a T does not equal zero,
the control law can be calculated on the basis of the second multiplicator:
Ds + S = 0
(9).
If the Equation (4) is differentiated and the Equation (2) is inserted into the recently calculated derivative, we obtain the following result:
S = G( f + Bu + DBu + d - xr ) (10).
After the Equation (10) has been inserted into the implementation condition of control law (Eq.(9)), the result is as follows:
-1 [
(GB) I G( f + DBu + d -
Ds (11).
x )-
r'
Since the term ( f + AB • u + d ) is unknown and not measurable it is, therefore, approximated with the neural network N = [ol ... o,]rby changing Equation (11) into:
u = -(GB)-1 [G( N -&r ) + Ds] (12).
Since the term ( f + DBu +d) is unknown and non-measurable, a classic supervised weight learning of the neural network can not be used. Therefore, a so-called on-line estimator has been developed for estimating a learning signal (the difference between the target and the output of a neural network). The result after Equation (4) has been differentiated, is as follows:
x = CT1š-
(13).
After the Equations (12) and (13) have been inserted into the basic equation of mechanism dynamics, the result is as follows:
S + Ds = G( f+DBu + d) - GN = G(Z- N)(14),
where is substituted Z = f + DBu + d .
By using the derivative of Lyapunov's function and the Equation (14), the condition, where the system controlled by the developed control law (Eq. 12) remains stable, is assured:
V = sTS = sTG( f+DBu + d - N)-sTDs < 0 (15).
The next condition (16) has been developed from Equations (14) and (15). To make V < 0 and consequently а ^ 0 possible, the condition expressed in Equation (16) has to be fulfilled during the time during in which the neural network is approximating the unknown part of robot dynamics ( f + ДBu + d) :
|G( f + DBu +d)-GN = |Ds + S| <|Ds\ (16).
To learn about the output layer of neural network with two layers, a modified BPG rule has been used:
net . = E wjoj + b, O- = g(net})
j
g(net) = 1- 2/(1-
E = (Da
s)T(Ds -
e~netЬ Dwij = s9E/ dwj (17),
s )/2 = (G(Z
N)T )(G (Z -
N) )/2
Dw, = edE/ dnet.dnet./ dw.. = edE/ dnet.o . = i j	i i i)	i J
= edE/ do. I do. I dnet.o . -
where: 11 1 J
- Dw j = edE/ do.g'(net.)o .
dE/dot = d/dot [(GZ- GN)T(GZ- GN)]/2 =
= d/do. [(GN)T(GZ- GN)]^ (18).
^ dE/doì = d/doì [(GN)]T (Ds + ss)
To learn the hidden layer of a neural network the traditionally back propagation rule is used.
From now, the detailed equations of control law for a SCARA mechanism with two degrees of freedom will be derived:
T = M в + h + Gf + Tn
(19),
where T, h, Gf and Tn are column vectors of the 2 by 1 dimension, M is the matrix of the 2 by 2 dimension, and в = [öj 02 ]T is the column vector of the 2 by 1 dimension of two axes of the SCARA. The previous equation can also be rewritten in the following form:
i = f( x, t) + B( x, t).u + d ( x, t) (20),
where:
= 1в ^ в, в,
i =|в, в0 в& в,
f =
(21)
A4 -
| h + IG f + F + Tn
1 0 0 1
M -1
and where M, -ft, Gt and Tn are estimated values of M, h, Gf and Tn.
Only nominal or average parameters are selected for the matrix Ml. This means that all 4 parameters of the matrix are constant while the robot hand is moving. This is, of course, only a rough simplification of how things really look; for
it is a common fact that the parameters of the matrix vary according to individual axis movements in a robot's working space. Because of this, the unknown variable part exists and is estimated by the neural network (Eq. 12). The dimension of vector /is 4x1 and the dimension of the matrix is 4x2. The control law u of the 2x1 dimension is illustrated in the following equation:
u = -(GB)- G N - xr ) + Ds ] where the variables are defined as:
0
(22),
G =
x = r
x = |q
r I 1r
Kp1		0	Kv
0	Kp 2		0
D		D1	0 "
		0	D2.
1r	0, 2r		01 1r
K
V 2
(23).
u2r
0 I
2r 1r 2r ] s = G( x - xr )
and
Coefficients of the matrices G in D should be selected in such a way that they enable the fastest convergence of neural network algorithm possible.
The column vector N is of the 4x1 dimension and represents the outputs of the neural network with 2=1-4. The learning procedure for all the weights of an output layer is:
Dw1 j = e j [ Kp1 0]( Ds-Dw2 j = e j [0 KP2 ](Ds-Dw3, = e j [ Kv1 0]( Ds-Kv 2 ](Ds-
'"3 j
V4 J
Dw4 , = e j [0
-s ) g '(net1)oj s)gl(net2 )oJ s ) g '(net3)oJ -s ) g '(net4)oj
where:
net. -i
= £ wfj + bj J
(24),
(25)
and where J = 1...20, I = 1...6, l = 1...4, and g' = 1 is the first derivative of the output function.
The neural network inputs are: two actual positions, two actual velocities and two differences between the reference and the actual position in joint coordinates.
Figures 6 and 7 show, that the developed CNNSMC is appropriate as a control scheme for the robot control. In the Figure 6, the positional error for the positional reference values in joint space is shown. The reference value is described as follows:
"1 ref
2 ref
3sin(0.5-1)
(26).
At the beginning of test experiment the positional error is proportionally huge. The reason is a random value of the neural network weights at the beginning of experiment. After first few moments neural network becomes stable and the positional error during movement of robot is approximately 5 mm. After 30 seconds the reference values in joint space ( 0 f, 02 mf ) are constant. Measured static error after 30 seconds is equal to zero.
For measurement of sum square positional error (SSEp) of the robot in Cartesian space, the Equation (27) is used. The result (Fig. 6) shows the convergence ability of a CNNSMC for the above described experiment.
2000
SSEp = ^ ((Xref i - Xi)2 + (Yf i - Y)2)(27X
where X
and Yre
are reference
vref_ i	ref _
coordinates in Cartesian space, X. and Y. are actual coordinates in Cartesian space and i is the number of measurements. It is determined with the equation:
Fig. 4. Control scheme of CNNSMC for SCARA teleoperation is shown
T,
2s

0.001s
■ = 2000
(28),
where Ts is the sample time and T is the time period of measurement.
5 IMPLEMENTATION OF EXPERIMENT
Telerobotic application is designed for educational purposes of students. Only predefined conditions for using it are reliable internet connection and installed LabVIEWs free library LVRun Time Engine. When a user is intending to take a control over the telerobotic application, the internet address of application is written to the internet browser program. The falling menu appears using the right-mouse click. There the option Request control of VI must be selected.
After control is established, the dialog window appears. The user name and user e-mail address must be written there, because the results of the experiment are sent via an e-mail. When the button SEND E-MAIL on the control panel is chosen, the results are sent to a defined e-mail address. The results are attached in folder rezultati.zip, which consists of:
rezultati.rtf the front panel and block diagram of the used virtual instrument is attached in this folder.
rezultati.datihe text-data of all measured signals is contained in this folder (in column form). The first column is time data, and the other columns contain the measured data. This format can easily be imported to MATLAB, where it can be observed in graph form.
14
12
E
Ш. 10
CL
Ш
5
6 4
2
Fig. 5. Measured dynamic error for reference position movement (Eq. 26) and measured static error (after the 30s period of experiment is elapsed) of the SCARA s tip
After the initial procedure of a web-based educational tool is completed, the front panel of the virtual instrument appears (Fig. 7). By using the switch on the top left rectangle of the control panel, it is possible to choose from three different running modes of the SCARA. The holding robot in the fixed-point is the first running mode, rotating robot links with sine reference in a joint-space is the second running mode, and defining the top coordinates of the robot arm in the Cartesian space by mouse clicking on the working area, which is presented as the grid on control panel, is the third running mode.
When the second running mode of the robot mechanism is chosen, the user is able to change the rotation frequency of the SCARA's axes from zero to half-radian per second in joint space and the amplitude of rotation from zero to pi radians and two graphs appear on the control panel of the virtual instrument. The referenced and actual position values of the first robot axis are shown on the upper graph, and the referenced and actual position values of the second axis are shown on the bottom graph. In addition, a positional error of robot, the SSEp in the joint space and the Cartesian space, values of some weights of neural network and the first output from neural network, can also be shown on both graphs.
In the middle of the control panel, within the small rectangles, the numerical values of the robot's position in the joint space, the values of some neural network's weights, the value of the first output from the neural network, the values of
0.45 0.4 0.35 0.3
ČT
0.25
Ш
и no
in 02
0.15 0.1 0.05 0
1 [s]
Fig. 6. The SSEp of SCARA s tip for the reference position movement (Eq. 26)
sum square errors in the joint space and Cartesian space, and the positional error value of the robot's top in Cartesian space, can be observed in the same control panel.
By using the numerical controls on the left middle rectangle on the control panel, the user can change the values of CNNSMC parameters at predefined intervals. The performance of CNNSMC changes when changing these parameters. The work is to find the best parameters, which would ensures that the performance of CNNSMC is optimal. Values of parameters are shown in Table 1.
6 HYPOTHETICAL USE OF THE DESCRIBED EXPERIMENT IN REALISCT MANUFACTURING ENVIRONMENT
The research field of production technology is often focused in the development of intelligent distributed production plants. The intelligent production plant is a machine, a part or whole
production plant, which is able to self-adapt of modifications of the production process and environment, which has also influence on the production process. By involving the optimization algorithms (e.g. genetic algorithm), the intelligent production plant is able to optimize the production process as well.
The developed experiment for teleoperating of the SCARA by CNNSMC is possible to implement in realistic manufacturing environment as a part of intelligent distributed production plant. The self-adapting ability of CNNSMR and possibility of teleoperation the experiment via the internet connection are the qualities, which classify the developed experiment in the class of the intelligent production plant.
If the control system dynamics of hypothetic production process is changed because of the load modifications and the variation of the friction forces value in bearings (usage and environment temperature), the developed CNNSMC is able to self-adapt on the dynamics changes. But, when the
Fig. 7. By using the standard Internet browser, the remote user is able to access the control panel of
experiment
Table 1. The remote user of experiment is able to set the values of CNNSMC's parameters within the predefined intervals
Symbol of parameter	Value of parameter interval	Short description of parameter
Kn	[2, 20]	velocity coefficient (1st degree of freedom)
KV 2	[2, 20]	velocity coefficient (2nd degree of freedom)
Kp1	[20, 200]	positional coefficient (1st degree of freedom)
Kp 2	[20, 200]	positional coefficient (2nd degree of freedom)
D	[20, 400]	dynamic coefficient (1st degree of freedom)
A	[20, 400]	dynamic coefficient (2nd degree of freedom)
nr	[2, 20]	number of neurons in hidden layer
eta	110-8	learning rate
system dynamics is changed even more, the self-adapting ability of CNNSMC can not optimize the regulation process. In that case, the remote user is able to fix the parameters of proposed controller to optimize the performance of CNNSMC again. The value of the reference position and speed of the production plant degrees of freedom is possible to set by teleoperation process as well.
7 CONCLUSION
The paper shows the intelligent interface, based on the neural network control approach for adaptation to environment disturbances (e.g. friction etc.) for telerobot application. The intelligent interface is used as a web-based educational tool for teaching neural network robot control techniques in telerobotics. When the robot dynamic is changed, this information is unknown to the user and because of this he is unable to dispatch it manually. This is not a problem for control however, because when CNNSMC is used, the approximation of changed robot dynamics is computed independently of the remote's user. The presented application's additional advantage is the possibility of on-line tuning parameters using CNNSMC, with the numerical controls of GUI, thus improving the performance of the robot's control and also making application more appropriate as an educational tool for students. With the possibility of on-line tuning parameters and observing responses of control using graphs and numerical indicators, students are able to learn how the parameters of CNNSMC affect performance.
8 REFERENCES
[1]	Hercog, D., Čurkovič, M., Edelbaher, G., Urlep, E. Programming of the DSP2 board with the MATLAB/Simulink. Proceedings IEEE ICIT
2003,	December 2003, p. 709-713.
[2]	DSP-2 web page: www.ro.feri.uni-mb.si/ projekti/dsp2.
[3]	Slotine, J., Li, W. Applied nonlinear control. Englewood Cliffs, NJ: Prentice-Hall, 1991.
[4]	Curk B., Jezernik K. Sliding mode control with perturbation estimation: application on DD robot mechanism. Robotica, vol. 19, 2001, p. 641-648.
[5]	Kaynak, O., Denker, A. Discrete time sliding mode control in the presence of system uncertainty. Int. J. Control, vol. 11, 1993, no. 7, p. 665- 678.
[6]	Rojko, A., Jezernik, K. Sliding-mode motion controller with adaptive fuzzy disturbance estimation. IEEE trans. ind. electron., vol. 51,
2004,	no. 5, p. 963-971.
[7]	Šafarič, R., Jezernik, K., Pec, M. Neural network control for direct-drive robot mechanisms. Engineering application of artificial intelligence, vol. 11, 1998, no. 6, p. 735-745.
[8]	Fujisawa, S., Obika, M., Yamamoto, T., Kawada, K., Sueda, O. Speed control of 3-mass system with sliding mode control and CMAC. IEEE International Conference on Systems, Man and Cybernetics, vol. 5, 2004, p. 44004407.
Paper received: 12.1.2007 Paper accepted: 19.12.2007
Artificial Neural Networks Application in Duplex/Triplex Elevator Group Control System
C. Erdem Imrak Istanbul Technical University, Turkey
Artificial neural networks can offer the better solution to the passenger call distribution problem when compared to the conventional elevator control systems. Therefore, the application of neural networks in elevator group control system is discussed. The significance of introducing artificial neural networks is presented. Elevator group control systems with neural networks can predict the next stopping floors to stop by considering what has been learnt by processing the changes in passenger service demand pattern. This paper deals with the use of artificial neural networks for the distribution of the most suitable cars to the floors by considering the passenger service demand. Artificial neural networks are applied in Duplex/ Triplex group control systems for improving passenger waiting time. The backpropagation algorithm is used for training neural networks. The elevator traffic analysis and simulation results are presented and compared to conventional elevator control systems. © 2008 Journal of Mechanical Engineering. All rights reserved. Keywords: elevators group control systems, neural networks, backpropagation
0 INTRODUCTION
The service provided by the elevator system of any modern building is necessary for the efficient functioning of the building. This service has to be not only reliable, but also satisfactory to the passengers. Elevator group control systems can be utilized to control the cars as a group and convey the passengers to their destinations comfortably and promptly. It is viewed as a combination of on-line scheduling, resource allocation and stochastic optimal control problems similar to robotics or automated manufacturing systems. However, elevator control systems are harder to deal with than other systems because the states of the system are dynamically changed in large state space and coming events are unpredictable in many cases [1]. In offering a solution to elevator control problem, the traditional control systems often yield unsatisfactory results because they lack in considering number of technical characteristics and possibilities to be taken in to account. They also possess limitations and their flexibility is still restricted even if they are adapted to utilize computers. The selection and distribution of the most suitable cars in the building is a function of the call assignment. Landing and car calls are often allocated to suitable cars by taking into account of
the minimum cost concept that operates by performing a trial allocation to all available cars and allocating the call to the car giving the lowest cost. The criteria for determining a suitable cost function depends on either quantity of service and/ or quality of service. The quality of service is a measure of the elevator capacity consumed to serve a specific set of calls, indicated by total journey times of all the cars [2].
Elevator control problem has been studied for a long time: in the last two decades artificial intelligence techniques such as fuzzy logic, neural networks, genetic algorithms and evolutionary algorithms were introduced. Although considerable research studied has been carried out to explore artificial intelligence applications in elevator control area [3] to [9], a little work has been done in neural network application in elevator group control.
Neural networks have been applied to tackle the selection and distribution of the most suitable cars which is a function of assignment of calls. Artificial neural networks which simulate the way that human neural networks operate and learn, and mimic the hardware structure of the brain in a simplified way consist of artificial neurons that perform a simple task and are fully interconnected [1]. A neural network is used to predict the number
*Corr. Author's Address: Istanbul Technical University, Faculty of Mechanical Engineering, Gümü§suyu, TR-34437, Istanbul, Turkey, imrak@itu.edu.tr	103
of landing calls likely to be made and the cars' likely destination for the time of service in order to optimize the car allocation efficiency. The prediction is made by combining predictions from historic data learnt during corresponding periods of previous days and real time data learnt during a short interval. They can be envisaged as a black box acting in response to a set of inputs of position, direction and calls produces an output representing the next stopping floor. Artificial neural networks can be placed in the conventional control system to improve the elevator performance [11].
In this study, the artificial neural networks has been inserted into an elevator group control systems for minimizing the average passenger waiting time and developed to model the behavior of the building population and also to adapt the algorithm to changes in passenger demand. For training the neural networks, the backpropagation method was used in the present work. The aim of this application is to optimize some of the parameters that a controller uses in order to minimize the average passenger waiting time. The data used here are from a simulator program developed by the author [9]. The elevator traffic analysis has been carried out using the software written in Turbo Pascal version 7.0 by examining the simulation results obtained. Simulations are run for several levels of interfloor traffic demand for the improved Duplex/Triplex control system (ICS). In the simulation an existing building having 3 elevators with 10 floors was selected as a model case. The normalized performance figures, with respect to the interfloor demand, are compared to other traffic control algorithms which are priority timed fixed sectoring algorithm (FSPTS), the fixed sectoring common sector system (FSCSS).
1 STRUCTURE OF THE ELEVATOR GROUP CONTROL
Building traffic becomes more and more diversified and complicated as building has more various functions and intelligence. Elevator group control algorithm, therefore, become more complicated as they become more intelligent. The goal of elevator control is to provide operational management of a group of elevators, by selecting cars to meet landing calls and achieve passengers' destinations pleasantly and promptly. An efficient elevator control system has four properties [2]: to
provide even service to every floor in a building, to minimize the passenger journey time in the car, to minimize the passenger waiting time, and to serve as many passengers as possible given time. The basic function of elevator control system is to assign an appropriate car to occurred landing calls to quickly respond them. The car to be assigned will be selected with consideration of many factors such as passenger waiting time, passenger journey time and energy consumption. Among them passenger waiting time is the dominant factor in elevator group control [8].
In this study, the Duplex/Triplex group control system is selected as an example of the conventional elevator group control system. The chosen system is suitable for groups of two or three cars in low-rise buildings and serves car calls and landing calls according to the directional distributive control principles by considering the current position, distance between landing call floor and car position, and the direction of car travel [2] and [9]. In order to perform directional distributive control, the nearest car continuously is to be searched with the call is allocated. Therefore it provides a single control algorithm designed to serve interfloor traffic in low-rise buildings.
Initial design of the elevator control system in building the main concern is to fulfill its requirements in conditions of predominantly up-peak traffic. Under up-peak traffic conditions, the Duplex/Triplex group control algorithm presents poor performance. To investigate improvement of this control algorithm, artificial neural networks were introduced into the algorithm and optimization was carried out to many parameters such as passenger waiting time and performance figure [10] and [11].
Landing calls, car calls and car positions are the inputs of the Duplex/Triplex control system. The elevator configuration such as number of floors, arrival rates, car capacities etc. is the parameters of the control algorithm. The outputs of the group control are the car direction and the next stopping floors. In this study, the next stopping floor problem in elevator traffic control is the explicit problems which has been tackled. It is stated that a given car system, with a known position, direction commitment and registered up and down landing calls as inputs find the next stopping floors as an output of the network. To solve this problem artificial neural network is developed
and applied to the Duplex/Triplex elevator group control algorithm.
2 NEURAL NETWORK APPLICATIONS ON
THE NEXT STOPPING FLOOR PROBLEM
The quality of service in multi-stores building is indicated by the average passenger waiting time. Thus, the aim of this study is reduce the average waiting time for the different traffic pattern. For this reason, it is proposed that neural networks learn the likely destinations from one floor to another and then use this information to allocate cares to calls and have been employed to predict elevator demand on a day-to-day basis. This study is an exploration in possibilities of the various fields in which neural networks can be applied in elevator control. Feedforward neural networks embedded in the car control and call distribution module is shown in Figure 1.
Neural network embedded elevator control module is situated in the simulation program and allocates the both landing and car calls to suitable car according to the control algorithm. Inverse-Stop Passenger method was proposed by Al-Sharif [12]. This module can continuously learn passenger arrival rate patterns throughout the day of the elevator system, and predict the passenger arrival rates for each floor and destination in the building. In an elevator system, when a landing call is registered at a given floor, the information the group control system receives is the floor identification and desired direction of travel. The next stopping floor as the destination is not known. The information could be used by the elevator control algorithm to send the best car to answer the landing
call. Conservative approaches can not provide such flexibility and autonomous behavior. As a preferred embodiment of this technique, population behavior is modeled using neural network approach as described by Rumelhart and McClelland [13]. Artificial neural network application in identification of vertical traffic pattern uses a feedforward neural network to perform this task [14]. Another neural network application in elevator control intends to be a new method using artificial neural networks to predict the response time of an elevator to a landing call.
2.1 The Structure of Artificial Neural Networks
Artificial neural networks consists several units that carry out a simple task and are fully interconnected. This layout allows reproducing some features of the brain that are not usual on computers [15]. A typical simulation of the action of a single neuron is shown in Figure 2. The strength of each synaptic junction is represented by a weight, with a positive sign for excitatory connections, and negative otherwise. The output of a neuron is easily computed by using vector multiplication of the input and weights, and adding the trainable bias to each neuron for speeding up convergence. A bias is included in the neurons to allow the activation functions to be offset from zero. This feature can be easily incorporated by adding an additional weight connected to +1, to each neuron. The additional bias term also determines the spontaneous activity of the neuron. This can also be achieved as setting the threshold values for the sudden onset of a high firing rate, hence the term non-linear threshold element.
Car allocation
Fig. 1. Overall diagram of an elevator control system
Input
Desired oulpul j
Fig. 2. A typical model of a neuron
After an input has been applied as a stimulus, it is propagated through each unit on output is generated. This output is then compared the desired output, and an error signal is computed for each output. The net input vector a. to neuron j is a non-linear function of the inputs x. of the neurons that are connected to j and the weights wjj on these conditions is:

w
i * j
(1).
Using sigmoid non-linearity at each node, the neuron has a real-valued output yj that is a nonlinear function of the net input is: 1
y> =
1 + e
(2),
where 0 it the threshold value for the artificial
j
neurons. Improving the performance sigmoid temperature a is used and slightly higher than one which makes the middle portion of the curve comparatively flatter. The node characteristics of an artificial neuron are thus determined by Equation (2).
The presented system can be adapted to do this purpose because of all the input values such as number of floors are defined as parametrically. A typical multilayer neural network used consists of input and output layers are configured by number of floors and two hidden layers as depicted in Figure 3. An increase the number of hidden layers
may improve generalization capacity. In the present case, after a few trials with the network with one hidden layer, good convergence could not be achieved. The network used, therefore, has two hidden layers as this guarantees it can learn any pattern. There are also no formulae to determine the number of neurons to arrange in hidden layers. The most used procedure is to try several configurations. To improve the generalization, the number of nodes in hidden layers are double that of the number of input layer nodes [16] and [17]. The development of successful neural networks requires considerable time and effort.
In fully connected neural networks, the inputs of the network gets in are: number of car calls, number of up landing calls, number of down landing calls and car directions. The hidden layers are computational layers and fully connected to input layer. There are no direct connections between the input and output layers. Case information flows forward from the input layer to the output layer through the hidden layers. Each layer consists of processing element called neurons [18].
The program uses the multilayer backpropagation method for determining the next stopping floors and calculates the outputs by multiplying the input matrix by the weight matrix to produce the output vectors. Then it sets the highest output to "1", while setting all other outputs to "0". It compares this output with the desired output, and then adjusts the weights.
a
i=i
Backward pass
Input layer
Output layer
First hidden	Second hidden
layer	layer
Fomard pass
Fig. 3. The architecture of multilayer neural network
2.2 The Backpropagation Algorithm	one. The learning algorithm is an iterative
algorithm, intending to reduce the error at the output. The sequence of steps involved is as follows.
Step 1: Initialize all of the weights to small random values. This breaks the network symmetry and ensures that the network does not become saturated by large values of weights.
Step 2: Apply an input vector to the input layer and specify the desired network outputs. The inputs can be new on each trial or samples from a training set can be presented cyclically until weights stabilize.
Step 3: Calculate the network outputs by performing a forward through the layers, using Eqs. (1) and (2) to calculate the outputs of individual neurons. The outputs of the input layer neurons are made equal to the corresponding component of the input vector.
Step 4: Adjust the weights of the network in a way that minimizes the error using recursive algorithm at the output nodes by:
Feedforward neural networks in this study use the backpropagation algorithm, as it is widely used and very popular. It provides predictive information and a systematic method for training multilayer neural networks that has been applied to the control of elevator call allocation [1] and [19]. It does not have feedback connections, but errors are backpropagated during training. Thus, data are processed from the input to output layer, whereas knowledge is performed by the algorithm of minimizing the square error by backward movement [20] and [21].
Backpropagation neural network is an exampled of supervised learning and applicable for the multilayer neural network for it considers weight parameter in all layers. There are two layers of neural cells between input layer and output layer as shown in Figure 3. The Sigmoid non-linear nodal function is used when the backpropagation learning algorithm is used to train the feedforward network which is intended for the next stopping floor prediction.
The backpropagation algorithm is also based on the generalized delta rule which calculates the gradient to find the deepest descent direction, searching for a minimum in the error surface. It calculates first the error in the output layer, then evaluating the contribution of each hidden neuron in that error. Thus, the corresponding weights are modified in order to reduce the mean squared error between the actual network output and the desired
Wj(t +1) = Wj(t) + gyj Ej + eAwy(t +1) (3),
where w.. is the weight from i-th node to j-th node, g is a learning rate which is normally sets less than 1 to prevent to oscillation of weight values, e is a momentum coefficient and is normally set to around 0.9 and where:
AWj(t +1) = Wj(t +1) - Wj(t)
(4),
and E is the error at j-th node:
Ej = - yj(1-yj )( yj - dj
for an output layer, and:
ej=- yj(l - yj )E E

(5),
(6),
for an intermediate layer neuron.
Step 5: Repeat by going to Step 2.
The major disadvantages of this algorithm are its slow convergence, the possibility of being stuck at a local minimum instead of a global and the memorization instead of generalization. Various methods have been suggested to speed up convergence time. Two of these involve the inclusion of a momentum term and a learning rate term in the weight adjustment equations. Using this method, the network tends to follow the bottom of narrow gullies in the error surface rather than oscillating from side to side.
There are also problems if the number of training pattern is small because this patterns using a big network in order to avoid the over fitting of the training data. Some modifications to the original algorithm tend to overcome these difficulties. Increasing the number of hidden layers can also improve the generalized capacity.
performance and scale the adjustments to weights. They involve change stop the weight update procedure, a variation to improve computational power and an explicit statement of the backpropagation algorithm for two hidden layers. A simple change to training law that sometimes results in much faster training is the addition of a momentum term. It is an attempt to try to keep the weight change in process and maintains the same direction the weight variation, ensuring stability in the learning. For constant learning rate, the number of epochs with respect to momentum can be seen in Figure 4. A momentum term is added to speed up the learning stage and takes value of 0.9 in this study.
Learning rate as another modification parameter can be used to control the amount which weights are changed during the training and be used as one of convergence. This convergence is drawn according to epoch with constant momentum parameter as shown in Figure 4. Therefore, the learning rate is set at 0.5 for training the network. As can be seen from Figure 4, care has to be taken in choosing the learning rate and the momentum, in order to ensure convergence. Higher learning rate values make convergence faster, but learning rate values higher than 0.5 must be avoided, as well as momentum values higher than 0.9.
2.3 Modification Parameters
2.4 Training the Network
Modification parameters can be made to the backpropagation algorithm which may improve its
The process of learning is known as training. There are three phases of training. In first phase,
0.4	0.6
Learning rate
(a)	(b)
Fig. 4. Modification parameters: (a) Momentum with constant learning rate (b) Learning rate with
constant momentum
Table 1. Parameters of neural networks training
	n 4	Neural model	backpropagation
Input layer	3n + 3	Topology	Feedforward
1st Hidden	2(3n + 3)	Supervision	Supervised
2nd Hidden	2(3n + 3)	Learning rule	Generalized delta rule
Output layer	n	I/O information	Binary
Input layer	Linear	Momentum	0.9
1st Hidden	Sigmoid	Learning rate	0.5
2nd Hidden	Sigmoid	Sigmoid temperature	0.5
Output layer	Linear		
Number of floors Number of layers
Number of neurons
Transfer function
supervised learning is used to train the network to predict the next stopping floor. The training set is first presented to the network, until it learns the internal representation of training pairs and then the error at the output node is reduces along the steepest descent direction. The initial weights and the thresholds are randomly generated at the beginning. At the end of the first phase of training, weights from the input layer to the hidden layer of the network are fixed. In phase two, the output layer of the network is retrained to emulate the existing controller. In phase three, single weights in the output layer of the network are perturbed, and the resulting performance is measured. The weights are then modified in the direction of improved performance. Neural networks to be trained wherein the parameters specified in Table 1 were employed.
Because of the size of the training pattern, it is kept a separate file produced by the training pattern program. A training pattern consists of training pairs of inputs and outputs. The inputs are
car calls, up and landing calls and car direction. All the inputs were introduced for several epochs. The outputs are the next stopping floors. A schematic diagram of an elevator system of a group of cars is illustrated in Figure 5. A group of cars has to handle a group of passengers in a building with a given number of floors. In this figure, the elevator cars are represented as boxes in the diagram. In addition, a black circle indicates a car call, a black triangle indicates a landing call assigned to an elevator, and a white triangle indicates an unassigned landing calls.
When considering even small number of stops, the size of the training pattern becomes extremely large to be calculated manually. Thus it decodes the direction and the position linearly, and the calls in binary [12]. The module is used to classify the input vectors and to convert them in binary information as shown in Figure 5, by a processing module active code/binary code. Thus, the neural network receives only binary data that represent a favorable situation for applications in
KJoor	Carl	Car 2	Car 3
© T			
			
©V			
	»		
	\		
		\	\
G л			>
Car calls	Car direction	Up landing calls	Down landing calls	Next stopping floor
<£K5><D<D©0)G	*	<£>®0©<2>Q)G	ФФФФФОС	
				
				
0000100	001	0010000	1000010	0000010
Fig. 5. Elevator system schematic diagram
large systems. This represents a reduction on the computational effort necessary to realize the training and an improvement on the quality of the analysis. This form permits to represent the contingences by a binary code and reduce the quantity that is an adequate way of providing a faster training and the analysis more reliable [9] and [12].
In order to train the neural network to be able to calculate the next stopping floors, the set of possible patterns has to be evaluated in advance, and stored in a file to be used later by the training program. While calculating the situations, it finds the next stopping floor corresponding to each simulation, and then writes all these variables to a file. In the sample file record given is Fig. 5, "0" represents no call and "1" represents a call to be registered. The length of record lines depends on the number of floors. The first 7 digits represent the car calls, the next 3 represents the direction (001 for down and 100 for up), the next 7 represents the up landing calls (each digit represents one call), and the next 7 represents the down landing calls and the last 7 digits represents the calculated next stopping floors. This format is used to give flexibility in order to changing the number of floors during input stage.
Each training pattern was randomly selected to be fed into the network, continuing until all the patterns were done with. This constitutes one epoch of learning. The network was trained for 1000 epochs, until the mean squared error dipped below 0.01 then the program stops [9]. After 1000 iterations the minimum reached which corresponds to the model that minimizes the fit for the validation data. Further iterations give an increasing criterion and the model becomes worse. Once the training stage is successfully finished and the trained network is tested, the network is employed in the prediction of the next stopping floors by considering what has been learnt in normal operation mode.
2.5 Overtraining the Network
In neural networks, one of the major pitfalls is overtraining, analogous to curve fitting for rule based trading systems. Overtraining refers to the reduction in generalization ability that can occur as networks are trained, and may be avoided by early stopping the training stage in time, before an
adaptation to the noise starts. It is difficult to detect an optional stopping moment without the use of a test set, pruning, reducing the size of a network after it is trained [12], [15] and [17].
Overtraining occurs when a network has learned not only the basic mapping associated with input and output data, but also the subtile nuances and even the error specific to the training set. Overtraining can be detected during the training stage by the use of a test set. Once the data set has been selected, it is separated into training and testing subsets. In the testing mode, the network is fed new inputs and utilizes the representation it had previously learned to generate associated outputs without changing its weights.
Overtraining is reduced by limiting the number of neuron to minimum necessary, increasing the training set size, using cross-validation training to identify and stop training before over fitting occurs. The cross-validation set was used to monitor generalization during training.
3 THE ELEVATOR CONTROL SIMULATION
In this section the elevator control simulation example is considered. A widely used method for elevator system performance prediction is the use of simulation. This can usually produce a set of curves, which describe the behavior of a certain variable against the change in system loading. One of the most useful curves which are used to study elevator system performance is the average waiting time curve versus interfloor demand. Interfloor demand as percentage can be found by dividing the passenger arrival rate by the handling capacity of the system [2]. Simulations are run for several levels of interfloor demand for the improved control system. At the end of each simulation, the average waiting time is recorded and the performance figure calculated by dividing the average waiting time by the interval to normalize it.
3.1 The Structure of Software
The simulation program is executed on a pseudo building with various numbers of elevator cars. The simulator was written by the author. A modular construction has been used as shown in Fig.6 as this allows ease of programming, extension, modification and defined interfaces
between modules and databases. Conversational programming techniques have been used to take advantage of an interactive digital computer. All inputs are checked for validity and where an error is detected the user is prompted as to the correct reply. The principal features of the simulator are that it can simulate up to 3 cars and 10 floors.
The lift configuration module is an interactive conversational program which allows the setting up of all the data defining the system to be simulated. Three cars with various contract capacities are used to service the whole building. In this program interfloor distance, contract speed, time for opening and closing the door, passenger transfer time are chosen by users. After that the training module is executed to be able to calculate the next stopping floors, the set of possible patterns has to be evaluated in advance, and stored in a file to be used in neural networks module. In this program, feedforward neural network is used to allocate the calls to the best cars during day period by means of weight matrix. The software uses the backpropagation learning algorithm to training the network.
The traffic analysis and control module which is embedded in the simulation program assigns the landing calls to the convenient cars. This module decides which car to be dispatched to which floor. In the traffic analysis and control module, the executive is the simplest part in the module and is used to control which sub-module is to run. After each sub-module has been used, the executive is re-entered to determine the next phase in a simulation. The executive also controls the access to the traffic analysis and control module. It generates landing and car calls. The car capacity is selected by users. The goal of the module is to obtain the next stopping floor, which is important in the allocation of the cars. In the first step, the landing calls and car calls, next arrival floors and next destination floors are determined. The number of passengers in the car, which one of the factors in selecting the next arrival floor is determined by means of Inverse-Stop Passenger method [12]. After two steps the next stopping floors can be determined using feedforward neural networks.
The output of simulation is the distribution of calls to the most suitable cars and car direction. The results are kept a separate file and saved in a file to be used later by the training program. Performance calculations have been arranged to present graphs to the user. They include round trip time, average waiting time and the normalized system performance figure which are used for a performance comparison of elevator control systems.
The flow chart of simulation program which employs both conventional Duplex/Triplex control and the neural network algorithm is depicted in Figure 6 [9].
3.2 Performance Criteria and Calculations
The performance objectives of an elevator system can be defined in many ways. One possible objective is to minimize the average waiting time, which is the time between a passenger's arrival and entry into a car. The passenger waiting time would be the best indicator of the quality of service that an installed elevator system could provide. To calculate the average waiting time according to car loads, one can find:
' 1.8P
2Ht + (S + 1)t + 2Pt AWT=-—---p
Fig. 6. The flow chart of simulation program
L
0.4 +
CC
-- 0.772
(7),
where P is the number of passengers enter the car at the main terminal, CC is the contract or rated capacity of car, H is the highest reversal floor, S is the expected number of stops, tv is the single floor transit time, ts is stopping time and t passenger transfer time and L is the number of cars within a elevator group. The highest reversal floor and the expected number of stops are determined by means of simulation program. The number of passengers in the car is calculated by means of Inverse-Stop Passenger method.
The normalized system performance figure is defined as the percentage of the simulated average passenger waiting time divided by calculated up-peak interval. The interfloor demand is defined as the interfloor arrival rate divided by calculated up peak handling capacity.
4 SIMULATION RESULTS
To investigate the improvement obtained by using the group control where neural network are embedded into the Duplex/Triplex control algorithm, simulations are run for several levels of interfloor traffic demand level for the improved Duplex/Triplex control system (ICS). The simulation program is executed for a building which population 1000 with 10 floors, three cars with contract capacity equal to each is used to service the whole building. Interfloor distance is fixed at 3.3 m, contract speed 1.0 m/s. Door operating time is 20 seconds while closing time 2.6 secs and passenger transfer time is 1.2 secs. Simulations output file keeps the record of the round trip time, the interval and the average passenger waiting time. The program runs on a 486PC with 8 MB RAM, 256MB hard disk configuration. The normalized performance
Table 2. Comparison of the system performance
figures, with respect to the interfloor demand are compared with other traffic control algorithms which are the fixed sectoring priority timed system (FSPTS) and the fixed sectoring common sector system (FSCSS) control algorithm. For a performance comparison, the results of each simulation performed under the same conditions have been tabled. Table 2 gives the comparison of the performance figure of the traffic control algorithms in a normalized form. Existing control algorithms as standards are used for comparison [9] and [10].
The simulation results showed the proposed method could provide a considerable improvement of the call assignment, if the method is implemented to an elevator control system for low-rise buildings in up peak traffic conditions. The ICS shows an improvement over the fixed sectoring priority timed system (FSPTS) and the fixed sectoring common sector system (FSCSS) control algorithm, when the interfloor demand exceeds 55%.
5 EXPANSIONS IN MATERIALS HANDLING SYSTEMS
Automated Guided Vehicle (AGV), one of the materials handling systems, is widely used in flexible manufacturing systems. Due to the increasing complexity of manufacturing environments coupled with the need for increased flexibility in AGV systems. The design and evaluation of AGV systems are highly complex because of the randomness and the large number of variables involved. Typically the design of material handling systems using AGVs includes determination of the vehicle guidepath layout, the traffic flow pattern, the number of vehicles required [22].
Interfloor demand	Performance figure [%]		
[%]	FSPTS	FSCSS	ICS
10	25.7	28.5	42.5
20	31.4	37.1	51.3
30	45.6	51.3	64.1
40	64.6	66.9	78.8
50	86.6	82.8	93.1
60	114.0	106.8	112.4
70	143.9	133.4	134.1
80	174.9	163.8	151.6
90	209.5	195.2	168.2
100	238.8	226.6	190.8
The problem of AGV consists of the operational control strategies of dispatching and routing of a set of AGVs and can be solved by using a lot of heuristic algorithms as Petri nets, fuzzy logic and neural networks. Hao and Lai [23] have described a neural network approach for the AGV problem and proposed neural network models to perform dispatching and routing tasks for the AGV under conditions of single or multiple vehicles which are based on Kohonen's self-organizing feature maps. Bostel et.al. have proposed a new navigation technique for AGV's based on neural networks [24].
The desired path along which the vehicle should move may be a fixed route, semi fixed route or an arbitrary route. Fixed route implies the installation of a fixed active guide like an energized cable or a passive one like a reflecting strip painted on the plant floor. It is essentially AGV programmed to convey materials or work in progress through a predefined route in horizontal environment. It is similar with elevator which carries the passengers and loads through a predefined route in vertical environment.
Vehicle travel time of AGV is a fundamental parameter for solving various flexible manufacturing system design problems. The feedforward neural networks embedded system presented in this paper possesses generality for transportation systems. It is mostly used in elevator engineering, but it could be used for single vehicle or multi vehicle AGV's control up to 3 vehicles in flexible manufacturing systems, if it is adapted and generalized for fixed route path problem using average waiting time and average travel time concepts as described in elevator traffic calculation. This control system could be applied the static production environment where the product mix or machine routings are assumed to be stable over time. It can also be used to determine and estimate the vehicle waiting times of AGV. It is believed that the models can be applied to a number of call-a-ride cases, such as AGV travel.
6 CONCLUSIONS
In this paper the outline of a program preparing feedforward neural networks has been introduced. It has been seen that the average waiting time which defines elevator performance goes down satisfactorily thanks to the neural network applied to elevator control systems and that it gives
better results compared to conventional methods. It could effectively support the conventional group control system and shorten the average waiting time by forecasting car position and using call distribution laws. The system results are shown in Table 2 compared to those obtained from conventional control algorithms. A 15 to 20% improvement in performance has been achieved by using this technique. It can be concluded that better results are reached when neural networks are combined with conventional techniques. On the other hand, the main drawback of the method is that it mainly applies to up-peak traffic. More research is also needed for down-peak, as it depends on building occupancy. This control system may be expanded to the other materials handling systems such as AGV's control problem with a few necessary arrangements.
Acknowledgement
The constructive comments of the referees are gratefully acknowledged.
7 REFERENCES
[1]	Miravete, A. New materials and new technologies applied to elevators. New York: Elevator World Inc, 2002.
[2]	Barney, G.C. Elevator traffic handbook. London: Spon Press, 2003.
[3]	Ho, M., Robertson, B. Elevator group supervisory control using fuzzy logic. Canadian Conference on Elevator and Computer Engineering 2, 11.4.4., 1994.
[4]	Prowse, R.W., Thomson, T., Howells, D. Design and control of lift systems using expert systems and traffic sensing. Elevator technology 4 (Ed. G.C.Barney). London: IAEE Publ, 1992, p. 219-226.
[5]	Qun, Z., Ding, S., Yu, C., Xiaofeng, L. Elevator group control system modeling based on object oriented Petri net. Elevator World, 49(2001)8, p. 99 - 105.
[6]	Chan, W.L., So, A.T.P. Dynamic zoning for intelligent supervisory control. Int. to Elevator Engineering, 1(1996), p. 1 - 10.
[7]	Siikonen, M-L. On traffic planning methodology. Elevator technology 10 (Ed. A. Lusting). London: IAEE Publ, 2000, p. 267-274.
[8]	Closs, G.D. The computer control of passenger traffic in large lift systems. PhD Thesis, University of Manchester, Institute of Science and Technology, Manchester, UK, 1970.
[9]	Imrak, C.E. Traffic analysis, design and simulation of elevator systems. PhD Thesis, Istanbul Technical University, Istanbul, Turkey, 1996.
[10]	Imrak, C.E., Barney, G.C. Application of neural networks on traffic control. Elevator technology 9 (Ed. G.C.Barney), IAEE Publ, London, 1998, p. 140-148.
[11]	Imrak, C.E., Ozkirim, M. The modeling and simulation of elevator group control systems for public service buildings. Preprints the 3rd IFAC Workshop DECOM-TT 2003, Istanbul, Turkey, June 2003, p. 159-164.
[12]	Al-Sharif, L.R. Predictive methods in lift traffic analysis. PhD Thesis, University of Manchester Institute of Science and Technology, Manchester, UK, 1992.
[13]	Rumelhart, D.E., McClelland, J.L. Parallel distributed processing. Vol. 1 and 2, Cambridge: MIT Pres, 1986.
[14]	So, A.T.P., Beebe, J.R., Chan, W.L., Liu, S.K. Elevator traffic pattern recognition by artificial neural network. Elevator technology 6 (Ed. G.C.Barney). London: IAEE Publ, 1995, p. 122-131.
[15]	Mukherjee, A., Deshpande, J.M. Application of artificial neural networks in structural design expert systems. Computer & Structures, 54(1995)3, p. 367-375.
[16]	Korn, A.G. Neural networks and fuzzy-logic control on personal computers and workstations. London: MIT Press, 1995.
[17]	Lisboa, R.G. Neural network current application. New York: Chapman and Hall Pub, 1992.
[18]	Alexandris, N., Chrissikopoulos, V., Vassilacopoulos, G. Lift - an expert system for lift system design. Elevator technology 2 (Ed. G.C.Barney). London: IAEE Publ, 1988, p. 1-9.
[19]	Rumelhart, D.E., Hinton, G.E., Williams, R.J. Learning representations by backpropagation errors. Nature, 323 (1986), p. 533-536.
[20]	Jovanovic, J., Krivokapic, Z. The application of an atypical neural networks when quantifying the modeling of environmental aspects. Strojniski vestnik - Journal of Mechanical Engineering, 52(2006)6, p. 392403.
[21]	Krivokapic, Z., Zogovic, V., Spaic, O. Using neural networks to follow the wear of a S390 twist drill. Strojniski vestnik - Journal of Mechanical Engineering, 52(2006)7-8, p. 437-442.
[22]	Koo, P-H., Jang, J. Vehicle travel time models for AGV systems under various dispatching rules. The International Journal of Flexible Manufacturing Systems, 14(2002), p. 249-261.
[23]	Hao, G., Lai, K.K. Solving the AGV problem via a self-organizing neural network. The Journal of the Operational Research Society, 47(1996), p. 1477-1493.
[24]	Bostel, A.J., Gan, W.W., Sagar, V.K., See, C.H. Generation of optimal routes in a neural network based AGV controller. Second International Conference on 'Intelligent Systems Engineering', Hamburg, Germany, September 1994, p. 165 -170.
Paper received: 19.3.2007 Paper accepted: 28.9.2007
Using Buckling Analysis to Predict Wrinkling in Incremental
Sheet Metal Forming
Samir Lemeš* - Nermina Zaimović-Uzunović University of Zenica, Mechanical Engineering Faculty, Bosnia and Herzegovina
The idea of using buckling analysis to predict wrinkling of sheet metal products manufactured by incremental forming is presented in this paper. We start from assumption that buckling analysis can be used to obtain geometry ofproduct after first forming operation. That shape and critical buckling force are calculated using finite element method in elastic domain. The analysis was performedfor various values of tool diameter. Results were tested through measurement ofproduct geometry after first forming operation. The analysis showed that upper tool diameter can be reduced without change in product quality expressed through number of wrinkles.
© 2008 Journal of Mechanical Engineering. All rights reserved.
Keywords: sheet metal forming, incremental forming, buckling, finite element methods
0 INTRODUCTION
Incremental sheet metal forming is usually being used in small-scale products manufacturing. These products are always manufactured with clamped edges and normally do not have geometrical imperfections related to wrinkling (Fig. 1).
On the contrary, when incremental forming is used to manufacture large products, such as spherical tank ends, wrinkling phenomena has significant influence onto quality of final products. The tool diameter is 2 to 5 times smaller than the forming part diameter, the edge of forming part is free, and the thickness of sheet metal plate is small. All these facts contribute to wrinkling occurrence. Figure 2 shows the machine for incremental
Fig.1. Incremental forming of small-scale products
forming of spherical tank ends up to 4 meters in diameter and sketch of 3 forming steps. The press uses two-part tool; upper tool is convex, and lower tool is concave, both having the same diameter. Step 3 in Figure 2 is repeated many times in order to incrementally obtain the desired form.
A number of researches were performed recently in order to minimize errors in incremental sheet metal forming. Mackerle [1] gave an exhaustive bibliography about application of finite element method in sheet metal forming simulation. The bibliography deals with material properties (texture, anisotropy, and formability), springback, fracture mechanics and calculation strategies, as well as with specific forming processes: bending,
Fig.2. Incremental forming press P2MF 200x4 -Sertom, Milan, Italy
*Corr. Author's Address: University of Zenica, Mechanical Engineering Faculty, Fakultetska 1, 72000 Zenica, Bosnia and Herzegovina, slemes@mf.unze.ba	115
extrusion, deep drawing, pressing, hydroforming etc.
Cao and Boyce [2] used forming limit diagrams as a criterion for disproportional tool load regimes. Their research showed that clamping force can be used to influence and to control wrinkling.
Wang and Cao [3] performed numerical analysis of wrinkling using modified energy approach, and investigated sensitivity of various input parameters and integration methods of finite element model onto buckling prediction. Their model was based on assumption that stress distribution in sheet metal is uniform. By defining appropriate boundary conditions, critical buckling stresses can be determined in sense of material properties, Poisson's ratio, sheet metal thickness and other geometry parameters. Wrinkling phenomena was than analysed by comparing thus obtained critical stress and real compressive stress calculated with finite element method. Their research showed that critical buckling stress of curved sheet metal depends on:
-	local curvature radius,
-	material properties,
-	sheet metal thickness,
-	dimensions and
-	load.
According to their research, the tool velocity has no influence on wrinkling.
Similar procedure was presented by Wang, Cao and Li [4], applied on bending of thin walled product edges. They offered an interesting conclusion that wrinkling is reduced when the length of bent edge is increased. That conclusion could be used to make an assumption: increasing starting sheet diameter could reduce wrinkling of spherical tank ends. They also investigated the influence of thickness and concluded that thickness has no influence onto number of wrinkles, but increased thickness leads to increase in critical length of bent edge.
Fig.3. Wrinkling of bent edges of thin walled products [4]
Pohlak et al. [5] observed this process in rapid prototyping. They performed numerical calculation in order to optimize manufacturing process.
Lee et al. [6] analysed the differences between static implicit and dynamic explicit integration method in time to compare numerical results in simulation of cold forming manufacturing process. They divided defects which occur in cold sheet metal forming into three categories: wrinkling, tearing (during deformation phase) and elastic springback (after tool removal).
Schafer and Schraft [7] compared various methods of flexible sheet metal forming with incremental deformation process. Their research was more focused on behaviour of manipulator robot rather than of the material itself. They proposed rapid tool movement, such that the tool acts onto sheet metal with strokes.
Ambrogio et al. [8] focused on material formability in incremental deformation process, particularly on estimation and compensation of elastic springback. They used integrated numerical/ experimental procedure to limit geometry deviations.
Kim and Yang [9] investigated buckling phenomena in deep drawing process using energy principle. They introduced "buckling factor" which is used to predict shape and location of wrinkles in sheet metal. Figure 4 illustrates wrinkling in deep drawing process.
Similar approach was used in this research. Buckling mode shapes were analysed in order to obtain tool parameters for incremental forming process resulting with minimum dimensional deviations of spherical tank end. The first forming operation, when buckling occurs, was analysed using finite element method.
Fig.4. Wrinkling in sheet metal deep drawing [9]
1 THEORETICAL BACKGROUND
Wrinkling occurs in areas which are not in contact with tool [3]. In incremental forming of spherical tank-ends, the only two contacts occur between the plate and the edge of lower concave tool, and between the plate and upper convex tool.
The following analysis assumes that tangential stresses before buckling are neglectable and plate thickness is uniform.
Well known Timoshenko's energy method with various combinations of boundary conditions was used in analysis of thin plates elastic buckling. The shape of deformed plate is presumed and critical buckling criterion can be obtained when internal energy of buckled plate equals the work performed by plane membrane forces. If internal energy for every possible deformation is larger than the work performed by membrane forces, the plate is in stable equilibrium.
The process usually starts with presumed function which describes plate deformation. It is usually a double sine function, with shape depending on plate shape (circular, rectangular, elliptical, etc.). For circular plate, the function can be:
d 2w 1 dw „
—2 +v--= 0,
dr2 r dr
w = w sin
m, n = 1,2,3,
(md)s
/ m(r- ra) Л
(rr- ra)
(1),
where wg is amplitude of wrinkles, m is number of wrinkles per perimeter, n is number of wrinkles in radial direction, ra is lower tool radius, rr is outer diameter of the plate.
The following boundary conditions apply:
w = 0,
(2)
Equations describing internal energy and work are very complex in this case and it is common to use numerical methods to solve this problem. Analytical calculations, as the one presented in [4], can show what influences critical buckling stress and wrinkling in incremental forming process. Proposed double-sine function, when visualised, corresponds to mode shapes obtained by experimental and numerical results presented here. Nevertheless, this function can be used only to confirm the assumption that buckling mode shapes play significant role in wrinkling occurrence during incremental forming.
2 EXPERIMENTAL SET-UP
The geometry of incrementally formed spherical tank (Fig. 5) was measured after first forming operation, in order to determine the buckling mode shape. Fig. 6 shows the layout of measurement points. Some parameters were varied to determine their influence onto buckling mode shape.
The measurement results are shown in Table 1. The number of lateral waves refers to modal diameters, while number of radial waves refers to modal circles. The press (Fig. 2) has only indication of pressure in hydraulic cylinder, which corresponds to press force. For deeper analysis, this pressure should be calculated or measured by means of force transducers.
Fig.5. Incremental forming of spherical tank end
Fig.6. Measurement points scheme
Table 1. Geometry measurement results
	Pressure	Lower tool	Upper tool curvature	Maximum	Number of	Number of
	[MPa]	radius [mm]	radius [mm]	amplitude [mm]	lateral waves	radial waves
a	10	350	700	163	4	0
b	10	350	650	156	6	2
c	10	350	250	136	6	0
d	13	350	250	185	4	0
Table 2. Material chemical composition
C	Mn	Si	P	S	Al	Cu	Cr	Ni	Mo	Ti	V	Nb	N
0.17	1.26	0.014	0.015	0.006	0.043	0.03	0.02	0.01	0.002	0.015	0.001	0.039	0.005
The desired buckling mode shape is the one with 0 lateral waves and 1 radial wave. Such a shape can be obtained by changing tool radii or pressure. Finite element analysis was then performed in order to obtain the combination of tool radius and press force which would lead to ideal buckling mode shape after first forming operation.
3 BUCKLING ANALYSIS
The tool radii, as the most controllable parameters, were varied in order to estimate their influence onto buckling and wrinkling. The
Fig.7. Simplified model used for numerical analysis
F [kN]
3000
r = 300 mm
50	100	150	200	250
f [mm]
Fig.8. Critical buckling force for f=50 to 250 mm
simplified model shown in Figure 7 was used. The plate support is circle with radius r (free rotations, fixed translations), the force F acts downwards on the surface with radius f The plate's outer radius is R and the thickness is d.
Material used is steel St 52-3 (S355J2G3) with chemical composition given in Table 2. The software used for analysis was UGS I-deas v.11. To relate the tool radius with wrinkling, it was important to sort the results according to buckling mode shapes.
The following values were used:
R = 1250 mm	d = 12 mm
f = 50 to 250 mm r = 200 to 600 mm F = 1000 kN
3.1 The influence of upper tool radius
The linear buckling analysis was performed for set of upper tool radii. Figure 8 shows analysis results for lower tool radius r=300 mm and variable upper tool radius f, sorted according to number of waves.
F [kN]
r = 400 mm
50	100	150	200	250
f [mm]
Fig.9. Critical buckling force for f=50 to 250 mm
0
200	250	300	350	400
r [mm]
Fig.10. Critical buckling force for
f =200 to 400 mm It is evident from Figure 8 that number of wrinkles does not change significantly with change in upper tool radius f
The analysis performed with lower tool radius r=400 mm and variable upper tool radius f showed even smaller influence of upper tool radius f onto critical buckling force (Fig. 9).
3.2 The Influence of Lower Tool Radius
Another analysis was then performed, varying the lower tool radius r values, keeping upper tool radius f=50 mm constant. Figure10 shows that increase in lower tool radius reduces critical buckling force.
Since results presented in Figure 10 showed increased influence of lower tool radius, similar analysis was performed with larger
F [kN] 5000
f = 50 mm
4500 4000 3500 3000 2500 2000 1500
f = 150 mm
F [kN]
3.3 Buckling Mode Shapes
The mode shapes which were used to sort the FEM analysis results are shown in Figure 12. The number of waves can be controlled with press force; the critical buckling force determines the buckling mode shape.
It is desirable to control manufacturing parameters in such a way that buckling mode shape shown in Figure 13 occurs. The analysis performed here showed that tool radius and press force can be used to change the mode shape, but it is rather hard to obtain desired shape, because different mode shapes have overlapping curves.
4 waves
0 waves
5000 4500 4000 3500 3000 2500 2000 1500 1000 500 0
6 waves
4 waves
200	250	300	350	400
r [mm]
Fig.11. Critical buckling force for f =200 to 400 mm upper tool radius and variable lower tool radius. Figure 11 shows analysis results for upper tool radius f=150 mm.
8 waves
Fig.13. Buckling mode shapes without waves
Fig.12. Mode shapes with 4, 6, 8 and 10 waves
Finite element method uses solution of eigenvalue problem to calculate both buckling and vibration mode shapes. In mathematical terms, every solver algorithm gives results sorted by size. It is important to visualise all numerical results, to be able to determine number of waves (wrinkles) from FEM results. For example, FEM solver gives 10 mode shapes accompanied with 10 critical force values, and visualisation must follow to determine the number of waves for each critical force value.
3.4 Experimental Results
To check the validity of numerical results, the geometry of incrementally formed spherical tank was measured after first forming operation. As numerical results led to conclusion that upper tool radius has low influence onto number of wrinkles, the measurements were performed for three different curvature radii of upper tool, between 250 and 700 mm. As Figure 14 shows, the larger curvature radius, the larger contact surface between upper tool and the sheet metal being formed.
According to results presented in Table 1, significant decrease in upper tool radius led to neglectable change in number of waves (wrinkles), as opposed to press force, which has strong influence onto buckling mode shape, in terms of both number of waves and buckling amplitude.
Therefore, both FEM and experimental analysis showed that upper tool radius has neglectable influence onto buckling mode shapes and number of wrinkles after first forming operation.
Figure 15 shows mode shapes modelled according to measurement results. The captions a, b, c, d correspond to mode shapes presented in Table 1.
ri < Г2
fi f2
Fig.14. Relation between curvature radius and contact surface
4 CONCLUSIONS
Linear elastic finite element analysis was used in this research in order to obtain buckling mode shapes for various radii of upper and lower tool. The analysis showed that number of wrinkles, which corresponds to number of waves in buckling mode shape, cannot be easily controlled by means of tool radius. The curve corresponding to zero-wave mode shape (the shape without wrinkles) propagates over entire domain and intersects other curves, as shown in Figures 8 to 11.
However, the upper tool radius is proven to have less influence onto critical buckling force and number of waves. That means the upper tool could be kept smaller, thus reducing the tool costs, keeping product quality the same.
It is hard to control process parameters to provide conditions when desirable mode shape occurs. The press force should be controllable if one wants to influence the buckling mode shape.
It is shown that conclusion given in [4] is not applicable in cases with free, unclamped edges.
Thin sheets wrinkle more easily and due to large lateral deformation. In contrary, thick plates wrinkle when deformation is limited; buckling occurs when major stress is compressive. In this manufacturing process upper surface is compressed and it buckles in a predictable mode shape. The higher mode shapes calculated by FEM correspond to wrinkles measured experimentally.
The major lack of this research is inability to control press force precisely, which is limitation introduced by press design. Further experiments should be performed on press with controllable force.
Fig.15. Measured buckling mode shapes
5 REFERENCES
[1]	Mackerle, J. Finite element analyses and simulations of sheet metal forming processes. Engineering Computations, Vol. 21 (2004) No. 8, p. 891-940.
[2]	Cao, J., Boyce, M.C. A Predictive Tool for Delaying Wrinkling and Tearing Failures in Sheet Metal Forming. Transactions of the ASME, Vol. 119 (1997), p. 354-365.
[3]	Wang, X., Cao, J. On the prediction of side-wall wrinkling in sheet metal forming processes. International Journal of Mechanical Sciences, 42 (2000), p. 2369-2394.
[4]	Wang, X., Cao, J., Li, M. Wrinkling analysis in shrink flanging. Transactions of the ASME, Vol. 123 (2001), p. 426-432.
[5]	Pohlak, M., Küttner, R., Majak, J. Modelling and optimal design of sheet metal RP&M processes. Rapid Prototyping Journal, 11/5 (2005), p. 304-311.
[6]	Lee, S.W., Yoon, J.W., Yang, D.Y. Comparative investigation into the dynamic explicit and the static implicit method for springback of sheet metal stamping. Engineering Computations, Vol. 16 (1999), No. 3, p. 347-373.
[7]	Schafer, T., Schraft, R.D. Incremental sheet metal forming by industrial robots, Rapid Prototyping Journal, 11/5 (2005), p. 278286.
[8]	Ambrogio, G., Costantino, I., De Napoli, L., Filice, L., Fratini, L., Muzzupappa, M. Influence of some relevant process parameters on the dimensional accuracy in incremental forming: a numerical and experimental investigation. Journal of Materials Processing Technology, 153-154 (2004), p. 501-507.
[9]	Kim, J.B., Yang, D.Y. Prediction of wrinkling initiation in sheet metal forming processes. Engineering Computations, Vol. 20 (2003), No. 1, p. 6-39.
Paper received: 13.3.2007 Paper accepted: 19.12.2007
Models of Reliability for Cutting Tools: Examples in Manufacturing and Agricultural Engineering
Predrag Dašić1* - Athanassios Natsis2 - Georgios Petropoulos3 'High Technological School in Kruševac, Serbia Agricultural University of Athens, Department of Agricultural Engineering, Greece 3University of Thessaly, Department of Mechanical and Industrial Engineering, Greece
Reliability of machining systems depends to a great extent on the reliability of cutting tool performance with the latter being the weak link due to wear Soil processing systems have many in common with machining systems and a major problem related to use of tillage equipment is ploughshare wear, as it markedly affects tillage quality and agricultural production economy. In the present study a machining and a soil processing example will be given for the determination and modelling of the reliability function R(T): for turning of 20CrMo5 steel alloy using Cubic Boron Nitride (CBN) and mixed ceramic cutting tools, and in a tillage operation using a straight toothed harrow. In the machining case, through comparative analysis of theoretical distribution models, which fit closer the experimental data, the Gaussian model is selected to represent the reliability function for both tool materials considered. In the tillage example, by carrying out an optimization procedure based upon the reliability function R(T) of the tillage system, a particular tool tooth geometry was found that establishes the maximum tooth working time possible for a reliable performance. © 2008 Journal of Mechanical Engineering. All rights reserved.
Keywords: machining systems, cutting tools, system reliability, agricultural engineering, tool wear
0 INTRODUCTION
Reliability of technological and production systems is considered of paramount importance regarding safe function and effective performance. The reliability of a machining system, e.g. a lathe with numerical control (NC), comprises the reliability of the following subsystems: machine tool, cutting tool and auxiliary accessories, which are connected in series. Nowadays studies on reliability of modern technological and machining systems: NC, CNC, DNC, FMS, RMS, IMS and so on, are certainly increasing.
The weakest link is the cutting tool, the excessive wear of which implies its ordinary replacement with obvious consequences on machining economy, and in this way the reliability of the tool expresses mostly the reliability of the whole system [1] to [4].
By the application of statistical methods, describing the breaking down of components of machining systems in the phase of effective performance, it is possible to determine the theoretical distribution, which fits best the experimental data.
The theoretical distribution models encountered more frequently are the following:
linear, exponential, hyper-exponential, normal (Gaussian), logarithmic normal, Weibull, Rayleigh, gamma, Erlang and Gumball or extreme (minimum or maximum) value of type I. The choice of theoretical distribution is checked through nonparameter tests: Pearson, Romanovski, Kolmogorov, Kolmogorov-Smirnoff and Mises.
The methodology followed in the present study assists stochastic modelling of tool life and can be seen in combination to tool life mathematical expressions (simple and generalized Taylor, Kronenberg etc.).
In this regard, the theoretical model of reliability of components in technical systems can be defined in one of the three following ways, as proposed by Dasic [4] and [5]:
•	on the basis of a previously defined theoretical distribution,
•	via choice of a theoretical distribution, which fits better the experimental data and
•	using distribution and complex [6] networks.
1 THE MACHINING CASE
The failures of cutting tools can occur at: cutting parts of the tool, tool holders, parts for
*Corr. Author's Address: High Technological School, Kosancićeva 36, 37000 Kruševac, Serbia,
dasicp@yahoo.com
mechanic fixation etc. However, it should be noted that the failures prevail at cutting parts of the tool. Under the term cutting tool failures are considered the facts, which have as a consequence the loosing of working capability. Besides, the cutting tool failures can be separated into:
•	failures with the possibility of renewing cutting tools working capacity. This can be realized by sharpening cutting parts of high-speed steel or introducing a new cutting blade, when coated or uncoated indexable hard metal inserts are used; or by means of repairing tool holder or other parts of the tool, and
•	failures without the possibility of renewing cutting tools working capacity i.e. failures, which cause rejection of cutting tools.
By observing the changes in time of some tool wear sensitive parameters during cutting process, for example flank wear width VB, a connection with deterioration of the process expressed by the reliability function R(T) can be obtained (Fig. 1).
The distribution may be given as the function of distribution F(T) related to the probability that the failure will appear until time instant t, i.e. the probability that the time of work without failures does not exceed a value t [1] to [3], [5] and [7] to [16]:
F(T) = P{T < t}
(1).
It is evident that at t=0, the value F(0)=0, and that at the value F(T) tends to be unity, i.e. F(&>)=1,
The probability of work without failures or the reliability function is a complement to the unreliability function and may be represented by the following equation:
R(T) = 1- F(T) = P{T > t}
(2).
It is evident that at t=0, the value R(0)=1, and that at the value R(T) tends to be zero, i.e. R(™)=0.
Fig. 1. Graphical illustration of the interaction between tool wear and reliability within a tool cycle
The failure frequency function is calculated according to the equation:
f(T) = F'(T) = -R (T)
(3).
And the failure intensity function is defined as:
m=
f(T) R(T)
(4).
In references [1] to [5], [8] to [13] and [17] to [22] the reliability of cutting tools is analysed during machining and the reliability in the case of using ceramic materials is studied in [23] to [25].
Metal cutting causes several types of wear which can be ascribed to a few mechanisms as illustrated in Figure 2 [26] and these mechanisms are defined by the international standard ISO 3685:1993 [27].
In references [28] and [29] typical tool wear features are illustrated in finish turning and VB and VB measures are defined.
max
The dependence of tool wear VB on machining time t is represented by a function:
VB = f(t)
(5),
which is of complex form and is derived experimentally.
In literature are presented different approximations of experimental tool wear curve with corresponding regression equations of different form.
Modelling of tool wear is performed mostly with parabolic functions [30] to [44].
In References [30] to [33], [36] to [38], [40], [42] and [43] proposed power and exponential functions are used to describe cutting tool wear.
Also, complex power-exponential functions have been proposed [30] to [33] and [36] and of polynomial m-th order form [30], [36], [45] and [46].
The necessary criteria for cutting tool failures are chosen:
•	for rough machining: on the basis of increased value of the cutting force or cutting temperature and
•	for finish machining: prescribed limits for quality of the machining products should not been exceeded, expressed by: roughness parameters, macro geometrical deviations, tolerance field etc.
Summing up, the description, stochastic or deterministic, of the reliability function can lead to more effective exploitation of cutting tools, as well as optimized selection of cutting tool geometry and material.
By applying a similar concept to soil processing systems and tools, the determination of the reliability function can lead to improvement of anti-wear performance characteristics. A tillage system has many in common with a machining system and a major problem related to use of tillage equipment is ploughshare wear, as it markedly affects tillage quality and agricultural production economy. A tillage operation system possesses a similar structure to a machining system and the weak link is again the tool, as its wear controls the reliability of the whole system to a great extent. Of
Chamfer
Thermal cracks in interrupted
cutting
Hi fib-speed steel
Carbide
Ceramic
land >
1.	Flank wear (wear
2.	Crater wear
3 Primary groove I outer diameter groove or wear notch I
Fig. 2. Tool wear mechanisms for different tool materials [26]
4.	Secondary groove (osidaiion wearl
5.	Outer metal chip notch
6.	Inner chip notch
course, in this case severe requirements for accuracy do not exist but excessive tool wear causes deterioration of technological and physiological parameters (power losses, insufficient cultivation depth etc.).
As aforementioned, by application of statistical methods describing the breaking down of cutting tools in the phase of their effective performance and exploitation, it is possible to determine the function of distribution that fits best the experimental data.
The present paper describes the introduction and relevant modelling of the reliability function in: firstly, machining by use of CBN and ceramic cutting tools, and secondly, in a tillage operation. In the first case a statistical approach is performed for finish turning of alloy steel 20CrMo5, which in the best way approximate experimental data on the basis of comparative analysis of the normal, logarithmic normal, Weibull and Gumbel extreme (minimum and maximum) distribution models. The selection of the suitable theoretical distributions is checked via the following non-parameter tests: Pearson, Romanovski, Kolmogorov, Kolmogorov-Smirnoff and Mises.
In the second case a deterministic analysis permits design optimization of the tillage tool. These two examples will be presented, separately, in the following.
Processing of experimental data by selecting a theoretical distribution is realized by DoRTSC2A software package [47].
1.1 Reliability Distribution of CBN Cutting Tools in Turning of Steel 20CrMo5
The cutting process in this case had the following characteristics:
•	cutting operation: external longitudinal finish turning, working diameter D = 50 mm and working length L = 78 mm;
•	material: steel 20CrMo5 (according to DIN standard) or C4721 (according to JUS standard) or 18CD4(S) (according to AFNOR standard), of 56-60 HRC hardness;
•	lathe: CNC lathe Max Muller MD5S with driving power 25 kW;
•	cutting tools: tool holder PTGNL2525M16 and the multi-bladed indexable inserts SPK-Feldmuhle TNMA160412T made of boron nitride WBN4 of special hexagonal
microstructure;
•	tool nose radius: r = 1.2 mm;
•	cutting conditions: cutting depth a = 0.2 mm, feed s = 0.12 mm/rev and cutting speed v = 90 m/min;
•	wet turning: processing with conventional coolant.
Failures of the CBN cutting tools have been observed during the turning operation, with the blunting of the cutting edge considered as a wear criterion (the maximum wear land width has been prescribed to VB = 0.4 mm).
The tests were replicated 11 times (sample size n = 11).
The time of work without failures t has been measured in [min] and has run from 18.08 min up to 23.60 min.
The resulting experimental distribution of the CBN cutting tools failures is positively skewed, i.e. asymmetrical on the left (ß = 0.7156) and close enough to normal (ß2 = 3.3329 or y2 = 0.3329).
For the scope of the present study the normal distribution is chosen owing to its symmetrical characteristics, simple form and unbiased character. The choice of the appropriate distribution model was checked by non-parameter tests: Pearson, Romanovski, Kolmogorov, Kolmogorov-Smirnoff and Mises and the best experimental data compliance was given by the normal distribution, the logarithmic normal distribution and the Guembel distribution of maximum values.
On the basis of comparative analysis [1], [5] and [47] of the distribution models: normal distribution, the logarithmic normal distribution and the Guembel distribution of maximum values, the chosen one is normal distribution, and its reliability function is as follows [3], [9] and [10]:
t 1.8111 ) The basic theoretical indicators of the reliability normal model are illustrated in Figure 3.
1.2 Reliability Distribution of Mixed Ceramic Cutting Tools for Turning of Steel 20CrMo5
The cutting process in this case had the following characteristics:
• cutting operation: external longitudinal finish turning;
Fig. 3. Graphic presentation of main indicators for the normal model of reliability of CBN cutting tools
in turning of 20CrMo5 steel
•	material: steel 20CrMo5 (according to DIN standard) or C4721 (according to JUS standard) or 18CD4(S) (according to AFNOR standard), of 56-60 HRC hardness;
•	lathe: universal lathe D-480;
•	cutting tools: tool holder CCLNR2525M16 and the multi-bladed indexable inserts CNGN 160812T02020 made of mixed ceramic SH1 from the firm SPK-Feldmuhle;
•	tool nose radius: r = 1.2 mm;
•	cutting conditions: cutting depth a = 0.5 mm, feed s = 0.16 mm/rev and rotation speed n = 280 rev/min, i.e. cutting speed v = 68.61 to 95 m/min;
•	mode of lubrication: dry cutting.
The trials were performed with sample size
n = 27.
The time for work free of failures has been measured from 12.05 min up to 19.41 min.
The resulting experimental distribution of the mixed ceramic cutting tools failures is positive, i.e. asymmetrical on the left (ß = 0.1518) and close enough to normal (ß2 = 2.398 or у = -0.662).
The normal distribution was chosen again owing to its symmetrical characteristics, simple form and unbiased character.
The selection of the appropriate distribution model was checked by non-parameter tests: Pearson, Romanovski, Kolmogorov, Kolmogorov-Smirnoff and Mises and the best experimental data fitting was given by the normal distribution, the logarithmic normal distribution, Weibull distribution and the Guembel distribution of maximum values.
R(T) [%]'	L i iS	' ЦТ)
F(T) [%]	g	
100 ■	25	2,5
80	■ 20-	■ 2,0-
60	■ 15-	■ 1,5-
40-	■ ю-	■ 1,0-
20	■ 5-	■ 0,5-
0	0	0
T [min]
Fig. 4. Graphical representation of the failure frequency function of the normal model of reliability for
turning with mixed ceramic tool
Based upon a comparative analysis [1], [5] and [47] of the theoretical distribution: normal distribution, the logarithmic normal distribution, Weibull distribution and the Guembel distribution of maximum values, the chosen one is normal distribution, and its reliability function is as follows [1]:
T-15.0383л
R(T) = 0.5 - Ф
1.7337
(7).
The basic theoretical indicators of the reliability normal model are illustrated in Figure 4.
2 TILLAGE EXAMPLE
Soil tillage operations consume large amounts of energy and cause significant wear to tillage tools. The latter results in deterioration of the overall performance of the plough i.e. higher energy losses demanding for higher fuel consumption, lower rates of work, decrease in tillage depth, time consuming changeover of the cutting edge and as a consequence, higher operation and production costs [48] and [49]. From the standpoint of agricultural production the proper soil manipulation directly affects the quantity of the production and finally the production cost of agricultural products giving also potential for savings.
Tillage quality expressed by a reliability function R = f(t, W) (t working time, W: wear land on the tool) will be directly affected by the wear law followed. If the latter is known experimentally or theoretically, an optimizing approach is possible regarding the appropriate selection of tool and soil processing conditions.
The basic expression combining reliability to wear and a work factor like plough area or time is as follows:
dR _ £ _dR dWj_ dr j _i dWj dT
(8),
where n is the number of worn elements that affect the working life of the machine part each time considered.
The expression for working life considering reliability is deduced from (8) as: Rd dR
(9),
Td _ I
R0 £
_1 д W1 dT
where Td is the working life and Ro, Rd the desired in the steady state and critical reliability values, accordingly.
Any change in the values of controlling parameters, such as tillage conditions (plough area, plough speed and tillage depth), the normal forces between the soil and the surfaces of the plough area, the proportion and mechanical properties of soil particles, the moisture content of the soil, and environmental effects and weather changes during the soil process will have an effect on the result of the integral (9), so we must consider both differentials of (9) separately.
Firstly, dWJdT expresses the wear rate for every element and depends on the material wear resistance and may be determined experimentally [50]. Secondly, the partial differential dR/dW may be determined analytically or experimentally on the basis of the wear behaviour; it is generally related to the element geometry and the material wear resistance.
Wear is related to the work factor by:
W _ f(T )
(10).
The wear gradient will be of the form:
dW
dT
1
(w)
(11).
The optimum design of the tool geometrical characteristics in view of wear would allow a reliable use for a longer working time and prolongation of its effective life.
In this way, the Equation (12) must be fulfilled according to the principle of extremes of functions
£	_ о
j_1 dWj dT
(12).
This equation denotes the conditions that rule the reliable performance of tillage machinery and considering the Equations (8) and (10), which provide the working time and wear rate, one can control the performance for varying process parameters. Also, it is evident that stabilization of the reliability function arises, when the latter is not affected by wear or if wear appears in many elements (surfaces).
Applying this concept to the anti wear behaviour of tillage tools an effective (minimum wear) tool geometric form can be selected. The
main points of a relevant analysis for a straight toothed harrow is described next; the detailed analysis is presented elsewhere [50].
A tool form that could be regarded optimum for given tooth material and wear resistance should allow maximum soil removal for the area Fd (Fig. 5). The working part of the tooth is described by a parabola. For any section, the relationship between the worn area F and the linear wear W is of the form:
F = Fd
1
W +—
H
1
1
yJ-
(13),
where x = F0/Fd = (0.1, 0.2 - 0.5), x = FJFi = 2 to 4, F0 is the smallest obelisk surface, Fd is the limiting allowable tooth working surface area, F is the tooth body cross section (16 *16 = 256 mm2), n = 0.5 to 2.5, W is the allowable tooth length regarding wear and H the tool length.
For 0.15 < x < 0.25, an empirically determined range combining the wear resistance of the tool material, the soil texture index and the geometry of the worn surface, a remarkable reduction in linear wear is achieved with regard to the maximum available material for re-sharpenings. According to this concept, to establish the maximum tooth working time possible for a reliable performance, the tooth must be of an "obelisk" shape with a limiting length corresponding to a final surface area equal to 20% of the limiting worn area. From the foregoing discussion, the introduction of the reliability function optimizes structural
parameters of tillage tools allowing the prolongation of their reliable working duration.
3 CONCLUSION
The use of reliability function in metal machining and soil processing, although obeying to different constraints, can lead to optimized anti-wear performance and better exploitation of tool useful life.
The two diverse examples discussed in this paper show that the determination of the reliability function in statistical or deterministic form permits another approach of tool wear in materials processing.
The normal distribution has been chosen for both types of tools considered. This fact reflects the stability of working without failures with these advanced cutting tool materials. All of the failure values are approximately symmetrically disposed around the mean value. A desirable state of surface finish can be maintained, as roughness of the workpiece increases regularly when approaching the end of tool cycle and this last interval can be beneficially exploited. In this regard, the latter could be considered in conjunction to the form of the tool reliability function and the degree of scatter and a supplementary criterion for predicting surface roughness, at any cutting operation, may be set.
It is evident that stabilization of the reliability function arises, when the latter is not affected by wear or if wear appears in many elements (surfaces).
Fig. 5. Harrow tooth modulation - Wear effects (ю /2: angle of the pyramidal modulation; H: tooth length; b: width of rectangular cross section; W: linear wear; f : cross section of the tooth holder; f : current tool cross section due to wear; fd: critical
tool cross section due to wear)
n
n
4 REFERENCES
[1]	Dasic, P. Determination of reliability of ceramic cutting tools on the basis of comparative analysis of different functions distribution. International Journal of Quality & Reliability Management (IJQRM), vol. 18, no. 4-5 (2001), p. 433-446.
[2]	Dasic, P. Examples of analysis of different functions of cutting tool failure distribution. Journal of Tribology in Industry, vol. XXI, no. 2 (1999), p. 59-67.
[3]	Dasic, P. Reliability analysis of the components of tribo-mechanical systems. International Journal of Applied Mechanics and Engineering - IJAME. vol. 7, Special issue (2002), p. 311-318.
[4]	Hitomi, K., Nakamura, N. and Inoue, S. Reliability analysis of cutting tools, Transactions of the ASME, vol. 10 (1979), p. 1-6.
[5]	Dasic, P. Algorithm approach to determination of reliability of components technical systems, Plenary and invitation paper. In: Proceedings of 5th International Conference Research and Development in Mechanical Industy - RaDMI 2005. Vrnjacka Banja, Serbia and Montenegro, 04.-07. September 2005, p. 34-45.
[6]	Wang, H. and Pham, H. Survey of reliability and availability evaluation of complex networks using Monte-Carlo techniques. Microelectronics and Reliability, vol. 37, ins. 2 (1997), p. 187-209.
[7]	Basu, A.P. Reliability and quality control. Amsterdam - New York - Oxford: Elsevier Science Publishers B. A., 1986.
[8]	Dasic, P. Analysis of the logarithmic normal model of the reliability of cutting tools made of cubic boron nitride for finishing turning of hardened steel. In: Proceedings of 7th International Conference on Accomplishments ofElectrical and Mechanical Industries - DEMI 2005. Banja Luka, Bosnia and Herzegovina, 27. and 28. May 2005, p. 209-216.
[9]	Dasic, P. Analysis of the normal model of the reliability of cutting tools made of cubic boron nitride. In: Proceedings of 1st DAAAM International Conference on Advanced Technologies for Developing Countries ATDC'02, Slavonski Brod, Croatia, 12. - 14. September 2002. p. 173-179.
[10]	Dasic, P. Reliability analysis of the cutting tools made of cubic boron nitride. In: Annals of DAAAM for 2002 & Proceedings of the 13th International DAAAM Symposium, Vienna, Austria, 23. - 26. October 2002. p. 125-126.
[11]	Dasic, P. Reliability analysis of the cutting tools made of multi-coatings hard metal. In:
Proceedings of 3rdInternational Conference, The Coatings in Manufacturing Engineering and EUREKA Partnering Event. Thessaloniki, Greece, 28. - 29. November 2002, p. 91-100.
[12]	Dasic, P., Djordjevic, A. Analysis of different functions of cutting tools failure distribution at the processing on the deep boring. Journal of the Balkan Tribological Association, vol. 9, no. 3 (2003), p. 370-380.
[13]	Dasic, P., Kuzin, V. V., Jecmenica, R. Reliability analysis of the ceramic cutting tools during turning. In: Proceedings of 5th International Conference on Tribology - BALKANTRIB 2005. Kragujevac, Serbia and Montenegro, 15.-18. June 2005, p. 273-276.
[14]	Dhillon, B.S. Quality control, reliability and engineering design. New York - Basel: Marcel Decker, 1985.
[15]	Kececioglu, D. Reliability Engineering Handbook - Volume 1. New Jersey: Prentice Hall, Englewood Cliffs, 1991.
[16]	Kececioglu, D. Reliability Engineering Handbook - Volume 2. New Jersey: Prentice Hall, Englewood Cliffs, 1991.
[17]	Klim, Z., Ennajimi, E., Balazinski, M., Fortin, C. Cutting tool reliability analysis for variable feed milling of 17-4PH stainless steel. Wear, vol. 195, no. 1-2 (1996), p. 206-213.
[18]	Lin, W.-S. Reliability study of cutting tool based on the reliability-dependent hazard rate function, Materials Science Forum, vol. 505-507, No. 2 (2005), p. 913-918.
[19]	Ning, F., et al. Study on reliability of ceramic cutting tools, Tool Engineering. Vol. 34, No. 8
(2000),	p. 3-5.
[20]	Sekulic, S., Nikolic, B. Influence of cutting conditions on reliability function parameters in turning. In: Proceedings of 6th International Conference Research and Development in Mechanical Industy - RaDMI 2006. Budva, Montenegro, 13.-17. September 2006, p. 345-349.
[21]	Vigneau, J., Bordel, P. Reliability of ceramic cutting tools. Ann. CIRP, vol. 37, no. 1 (1988), p. 101-104.
[22]	Wang, K.-S., Lin, W.-S., Hsu F.-S. A new approach for determining the reliability of a cutting tool. International Journal of Advanced Manufacturing Technology, vol. 17, No. 10
(2001),	p. 705-709.
[23]	Jacques, L. Ceramics reliability: Statistical analysis of multi-axial failure using the Weibull approach an the multiaxial elemental strength model. Journal of the American Ceramic Society, vol. 73, no. 8 (1990), p. 2204-2212.
[24]	Jacques, L. Statistical analysis of bending strengths for brittle solids: A multiaxial fracture problem. Journal of the American Ceramic Society, vol. 66, no. 3 (1983), p. 177-182.
[25]	Jacques, L. Statistical Approaches to Failure for Ceramic Reliability Assessment. Journal of the American Ceramic Society, vol. 71, no. 2 (1988), p. 106-112.
[26]	Kalpakjian, S. Manufacturing processes for engineering materials. 3rd Edition. Menlo Park (California): Addison Wesley Longman Inc., 1997. 950 p.
[27]	ISO 3685:1993 Tool life testing with singlepoint turning tools. Geneve: International Organization for Standardizations (ISO).
[28]	Chou, Y. K., Evans, C. J. Hard turning of M50 steel with different microstructures in continuous and intermittent cutting. Wear, vol. 255 (2003.), p. 1388-1394.
[29]	Chou, Y. K., Evans, C. J. Tool wear mechanism in continuous cutting of hardened tool steels. Wear, vol. 212 (1997), p. 59-65.
[30]	Dasic, P. Analysis of wear cutting tools by complex power-exponential function for finishing turning of hardened steel 20CrMo5 by mixed ceramic tools. The Annals of University "Dunarea de Jos " of Galati, Fascicle VIII Tribology, vol. 5 (2006), p. 54-60.
[31]	Dasic, P., Jecmenica, R., Radovanovic, M. Analysis of the wear tools for turning by modern cutting tools. Journal of Optimum Technologies, Technologic Systems and Materials in the Machines Building Field - MSTM, vol. 11, no. 1 (2005), p. 213-218.
[32]	Dasic, P., Nedic, B., Jecmenica, R. Analysis of function approximation of the wear tools for turning of the cast iron GG-25 by nitride ceramic cutting tools. Journal of Modelling and Optimization in the Machines Building Fields -MOCM, vol. 12, no. 2 (2006), p. 46-60.
[33]	Dasic, P., Radovanovic, M. Analysis of the wear cubic boron nitride for turning of the hardened steel. International Journal of Applied Mechanics and Engineering - IJAME, vol. 7, special issue (2002), p. 341-346.
[34]	Dawson, G. Machining hardened steel with polycrystalline cubic boron nitride cutting tools. Dissertation. Georgia (Atlanta): Georgia Institute of Technology, 2002. 264 p.
[35]	Dolinsek, S., Sustarcic, B., Kopac, J. Wear mechanisms of cutting tools in high-speed cutting processes. Wear, vol. 250, issues 1-12 (2001), p. 349-356.
[36]	Ekinovic, S. Cutting (in Bosnian language). Zenica: Dom štampe, 2001, p. 307.
[37]	Janac, A., Batora, B., Baranek, I., Lipa Z. Machining technology (in Slovakian language). Bratislava: Slovenska technicka univerzita (STU), 2004. 288 p.
[38]	Kalajdzic, M. Tehnology of metalworking I (in Serbian language). 9th Edition. Belgrade: Faculty of Mechanical Engineering, 2004. 407 p.
[39]	Kopac, J. Cutting tool wear during high-speed cutting. Strojniški vestnik - Journal ofMechanical Engineering, vol. 50, no. 4 (2004), p. 195-205.
[40]	Kopac, J. Cutting (in Slovenian language). Ljubljana: Univerza v Ljubljani, Fakulteta za strojništvo, 1991.
[41]	Poulachon, G., Moisan, A., Jawahir, I.S. Tool-wear mechanisms in hard turning with polycrystalline cubic boron nitride tools. Wear, vol. 250, issues 1-12 (2001), p. 576-586.
[42]	Radovanovic, M. Metalworking technology, Cutting work (in Serbian language). Niš: Faculty of Mechanical Engineering, 2002. - 328 p.
[43]	Stanic, J. Theory of metalworking I (in Serbian language). 2nd Edition. Belgrade: Faculty of Mechanical Engineering, 1989. 268 p.
[44]	Taylor, F. W. On the art of cutting metals, Transactions of the ASME, vol. 28 (1907.), p. 31-279.
[45]	Soković, M., Mikula, J., Dobrzanski, L. A., Kopac, J., Kosec, L., Panjan, P., Madejski, J., Piech, A. Cutting properties of the Al2O3 + SiC(w) based tool ceramic reinforced with the PVD and CVD wear resistant coatings. Journal of Materials Processing Technology, vol. 164165 (2005), p. 924-929.
[46]	Tadic, B. Contribution to analztical definition ofwear curves (in Serbian language). Journal of Tribology in Industry, vol. XII, no. 4 (1990), p. 117-120.
[47]	Dasic, P. DoRTSC2A: Methodology and Software package for Determination of Reliability of Technical System Components on the Basis of Choice ofTheoretical Distribution, which in the best way Approximate Experimental Data, on the Base of Comparative Analysis. ver. 3.0. Krusevac, 2002.
[48]	Bhakat, A.K., Mishra, A.K., Mishra, N.S., Jha, S. Metallurgical life cycle assessment through prediction of wear for agricultural grade steel. Wear, vol. 257 (2004), p. 338-346.
[49]	Natsis, A., Papadakis, G., Pitsilis, J. The influence of soil type, soil water and share sharpness of a mouldboard plough on energy consumption, rate of work and tillage quality. Journal of Agricultural Engineering Research, vol. 72 (1999), p. 171-176.
[50]	Natsis, A., Petropoulos, G. Optimal tribological design of soil processing tools. Journal of the Balkan Tribological Association, vol. 11, no. 4 (2005), p. 491-498.
Paper received: 20.2.2007 Paper accepted: 28.9.2007
Simulation Methods in Shipbuilding Process Design
Boris Ljubenkov1* - Goran Đukić1 - Marinko Kuzmanić2 'University of Zagreb, Faculty of Mechanical Engineering and Naval Architecture, Croatia 2ProMar Design, Kaštel Novi, Croatia
The paper aims at presenting the use of simulation methods in shipbuilding production process design. The basic design procedure principle is described, the main elements of the designing spiral, flow diagrams and interrelations between individual activities determined. The flow diagram indicates a point in which the simulation results affect the design flow.
A simulation program package was used to create a shipbuilding production process model. An example of possibilities and methods of presentation of the program package outputs is given, and the advantages ofusing simulation methods in the shipbuilding production process are described in a conclusion. © 2008 Journal of Mechanical Engineering. All rights reserved. Keywords: shipbuilding, process design, simulation methods
0 INTRODUCTION
The shipbuilding production is a complex and lengthy process, which demands careful planning and timely decision-making. Characteristic of an intermittent process like shipbuilding is a large number of working activities of different duration. Thus, it is necessary to ensure work-in-process storage areas, which demand an adequate space and a rather intensive use of transport devices. Position of machines, transport equipment and other devices within the process does not change. The major changes in the process are caused by diverse production programs. Products of different purpose and geometries pass through the same production process and incur different work-loads on workshops, equipment and work-in-process storage areas, which might cause interruptions in production [1].
Problems encountered in complex systems are efficiently resolved by simulation methods. Mathematical and statistical analyses do not play a major role in such systems, since they are not realistically described by mathematical equations. The literature does not offer an unequivocal definition of simulation. The term "simulation" would mean imitation, and the simulation procedure, according to one of definitions, is a set of activities ranging from real system modelling to experimenting with a model and analysis of results [2]. The simulation enables prediction of steps taken in real production process and
recognition of unfavourable situations, such as interruptions or bottlenecks in production. It is also possible to monitor the effects of parameter changes on overall process before the production in workshop starts, when error corrections are much more difficult and expensive.
Papers about simulation methods in the shipbuilding are scarce. In those very few authors' effort is to present and encourage the use of simulation method as a tool for easier decision making in shipbuilding production process. Medeiros et al. [3] and Williams et al. [4] are presenting their work on systems for modeling and visualization of production process and the importance of the internet as a communication tool. Kiran et al. [5] introduce hierarchical approach to modeling. Overall simulation model of the production process is divided in number of smaller models that are created separately. Smaller simulation models could be then controlled and analyzed easier. Their approach also includes methods that integrate these sub-models into an overall model in order to run different scenaria and identify global performance measures.
Research of McLean and Shao [6] presented overview of the generic simulations of shipbuilding operations. They are concerned to control and follow-up of labour resources and costs in complex system like shipbuilding. Dain et al. [7] present stochastic simulation model for determination, follow-up and control of costs, activity delay or labour resource shortage because they could cause
*Corr. Author's Address: University of Zagreb, Faculty of Mechanical Engineering and Naval Architecture, Ivana Lucica 5, 10000 Zagreb, Croatia, boris.ljubenkov@fsb.hr	131
huge disturbance in production. Nedess et al. [8] emphasize, among other things, that due to the fact that required space is a resource for time-critical processes with restricted capacities on shipyards and other shipbuilding companies. The space allocation has to be handled as one of the most important tasks in production planning, and this aspect was implemented in the simulation module for space allocation within the Simulation Toolkit Shipbuilding (STS). They also highlight the need for special simulation modules for shipbuilding industry. Steinhauer and Meyer-König [9] present importance of simulation and STS in sub-assembly and assembly stage planning.
1 SHIPBUILDING PROCESS DESIGN
A shipbuilding process design consists of elements each of which has its appropriate requirements and limitations. Elements of the shipbuilding process design are production concept, product analysis, material flows, building dynamics, arrangement of working areas, process flexibility, management, logistics, human resources and investment. Their interrelation is complex, and the ultimate solution is looked for in iteration of several steps [10]. Each next iteration step brings new solutions for the design elements, which affect other design elements. The final result is a compromise design of an overall shipbuilding process, which does not offer the best solution for each individual design element.
The shipbuilding process design procedure may be expanded by application of simulation program packages for modelling, visualisation and experimenting with a model, which enables prediction and avoiding of adverse effects the parameter changes have in a complex production process such as shipbuilding.
When the data obtained by analyses of product, production process and material flow are used, the simulation procedure enables creation of production process and final product model. Performing a large number of experiments with models in a short period of time will create a considerable database for quality and informed decision-making on solutions for the production process.
Figure 1 shows a block diagram with shipbuilding process design elements and their interdependence achieved by two-way relations.
Interdependence of elements is complex, and their relations in the design stage are open in both ways until the final solution has been found.
Simulation in the shipbuilding process design consists of a preparatory and executive process.
The preparatory procedure includes calculation of necessary model element and creation of the production process and final product models.
Executive part of the simulation procedure starts with the experiment scenario, which defines a research purpose, input data and analyses parameters, such as the equipment efficiency or experiment duration. The procedure continues with experiment with a model and analysis of output data, which results in conclusions. If the conclusions have determined that the research were to continue with a new experiment, the procedure is repeated, i.e. the changes are made in product analysis, production process and material flow. These changes affect the production process and final product models, which initiate a new experiment. When all the options are researched and final conclusion brought, the results are presented to responsible management, which makes decisions on all the shipbuilding design elements.
2 AN EXAMPLE OF SIMULATION METHOD
USE IN SHIPBUILDING WORKSHOP PRODUCTION PROCESS MANAGEMENT
The simulation method use in management of production processes will be shown using an example of situation analysis for a shipbuilding workshop [11]. The data obtained from the shipyard will be used to build a simulation model of the shipbuilding workshop, define the production program and carry out computer-supported experiments. Analysis of results of an individual experiment and comparison of results of a larger number of experiments creates a database facilitating quality decision-making on the production process changes. Taylor ED software is used for modelling, visualization and results analysis.
2.1 Simulation Program Package Taylor ED
Simulation software Taylor ED is object-oriented program package for modeling, visualization and analysis of the production
Fig. 1. Shipbuilding process design elements obtained by simulation procedure
process. In the program package whole production process and products of the workshop are created. A basic element of the simulation model is an atom. The atom represents single machine, tool, and transport device, necessary areas like interim stores, semi-products and products of the workshop. There are 80 defined and created atoms in the program package.
In the program package atoms are arranged on workshop layout. Their geometric characteristics and functions are defined, as well. Relations between atoms are defined to determine task scheduling in the workshop. Connecting entry and exit channels of atoms define the routing - the flow of the products through production process in a simulation model. Function of the atom is operation that has to be done. It could be set, for example, as duration of machine operation or speed, load and unload time of transport device.
2.2 Shipbuilding Workshop Simulation Model
The simulation program package models the production program with elements that represent equipment, devices and transport devices. Each
element is defined by its geometry (sizes of equipment and devices), position in space and operating characteristics (speed of transport devices or fabrication duration on equipment). Figure 2 presents a model layout.
Cut lengths of plates, profiles and bars are defined by their geometrical characteristics and ship structure analysis. Through cutting machines possibilities, duration of cutting is determined for each machine. Analysis of the cutting duration real values is basis for value distribution definition. Log-normal distribution is very close to real values and it is set in cutting machines atoms in the simulation model.
Basic areas are material stockyard and shipbuilding workshop. Plates and profiles are arranged on the material stockyard, which is equipped with a portal crane. Portal crane transports a material on the prefabrication line. After that plates and profiles are going to the shipbuilding workshop.
The shipbuilding workshop consists of 5 naves. The first nave is used for profile cutting on a robotised line, which consists of a line for the profile edge cleaning and a robot station for plasma
Fig. 2. Shipbuilding workshop simulation model
134
cutting. Profiles are transported in the nave with bridge crane, roller conveyors and transverse chain transport device.
The second nave is used for disposal of tools and as an auxiliary area off the main material flow. Being less important, this area is not modelled in the program package.
The third nave is used for shaping of earlier cut plates. The plates are shaped in press, and transported in the nave with bridge crane and trolley.
In the forth nave, the profiles are bent on bending machines and transported with bridge crane and trolley.
The plates are cut in the fifth nave. Three cutting machines are used, one for plasma and two for gas cutting (Numorex and Corta). The plates are transported with bridge crane and trolley.
2.3 Production Program
2.4 Experiments with Model
The modelling process is followed by launching of an experiment. The visualisation module used to monitor the experiment progress is started at the same time. It enables monitoring of what is happening with material, semi-products, equipment, devices and transport devices in any part of the process.
The paper describes two experiments with models. Having in mind complexity of the analysis, monitoring of process and results analysis will be limited to the profiles and plates cutting processes. The efficiency of equipment, devices and transport devices, situation at the interim store and quantity of fabricated material are monitored.
The first experiment monitors and evaluates situation in the basic simulation model for the shipbuilding workshop. Based on the analysis results, the changes and improvements are
To initiate the experiment with the simulation model, it is necessary to define the workshop production program, necessary material quantities, and sequence of fabrication. In this example, the production program includes plates, profiles and bars used for the hull building.
Quantity of material is determined from characteristic sections. The central part doublebottom sections and parts of the double-bottom with bilge of an oil product carrier were selected. Geometry and quantity of plates, profiles and bars was determined from the specification of material, same as the fabrication procedures. These data are used to define the scenario of the material input into the shipbuilding workshop. Figure 3 shows mass percentages of plates, profiles and bars in total quantity of material. In total quantity of material, plates account for 84%, and profiles and bars for 16%.
Fig. 4. Distribution of material according to the working site in the first experiment
Robot
Fig. 3. Mass percentages of plates, profiles and bars in total quantity of material
Fig. 5. Profiles cutting robot status during the first experiment
Simulation Methods in Shipbuilding Process Design
Numorex
Fig. 6. Status of cutting machines and plate edge preparation bench during the first experiment
proposed to be included in the new simulation model.
The second experiment monitors effects of the proposed changes on situation in the nave in which the profiles and plates are cut. The proposed improvements are estimated by comparison of obtained results.
Since an experiment with the simulation model in the program package is running continually and permanently, the experiment duration needs to be limited. The experiment duration is set at 1920 hours, i.e. according to the number of working days and effective working hours, at six months work in a shipyard.
2.4.1 The First Experiment with Model
Figure 4 shows distribution of material according to its fabrication site in the first experiment. Plates are commonly cut by plasma cutting machine because of its cutting speed. Profiles are cut by robot, and bars are usually cut by hand.
During the experiment, profiles and bars cutting robot operation was monitored, and the robot status is presented in Figure 5. The machine was busy 52% and idle 48% of time, and it cuts
8641 profiles and bars. Hence, the machine is not sufficiently busy and it could fabricate an additional quantity of profiles and bars.
In plates cutting workshop (nave 5), all the cutting machines and the bench on which semiautomatic cutting machines are used to prepare the plate edges for welding were monitored. For the experiment, 30% plates were planned to be forwarded after machine cutting to the benches for edge preparation for welding. Number of fabricated plates, profiles and bars are in Table 1. During the experiment, 6403 plates were fabricated at the cutting machines. Figure 6 shows status of machines and bench during the experiment. It is noticed that machines are blocked for quite a long time, which means that the plate is fabricated and waiting for transportation. The transportation in the workshop is done with a bridge crane, which is expected to perform numerous tasks. It is also noticed that in the interim store, near the edge preparation bench, a large quantity of material piles up and that it needs to be released of storage load.
Verification of the simulation model is performed and the model behaves with expectations. Validation of the model is established for 6 months period by comparison of experiments
136
Table 1. Number of fabricated plates, profiles and bars in first experiment
Cutting machine	Plates	Profiles and bars
Robot	-	8641
Plasma	3487	-
Numorex	1807	-
Corta	1109	-
Total	6403	8641
Semi-automatic cutting machine	667	-
results and real number of worked elements. Simulation model expresses real situation in the shipyard workshop very well.
According to the analysis results, the process improvements are proposed which are to be incorporated into the new simulation model and new experiment.
The proposals include:
-	Increase in profiles cutting machine efficiency by delivery of a larger quantity of bars for fabrication. In new experiment, the robot is fed with 60% bars compared to 22% from the first experiment.
-	Moving the plate edge prepared for welding to one of oxygen cutting machines so that number of plates that go to the edge preparation bench is reduced from 30% to 25%.
-	Fitting the plate-cutting workshop with a new bridge crane to share the transportation load with the existing crane.
2.4.2 The Second Experiment with Model
Introduction of changes into the model results in a new scenario of material input into the shipbuilding workshop. New distribution of
material according to the working place is shown in Figure 7.
During the experiment, same as in the previous experiment, operation of equipment, devices and transport devices, and quantities of material fabricated in the profiles and plates cutting workshops was monitored.
Figure 8 shows the profiles cutting robot status. The machine was busy 90% of time and it cut 14756 profiles and bars, which is an increase of 38% and 70% respectively compared to the first experiment.
Number of fabricated plates, profiles and bars is given in Table 2. The machines fabricated 7248 plates, which is 13% more compared to the previous experiment. Cutting machines status is shown in Figure 9. Blocked status was reduced by 25%, and efficiency increased because the crane load was reduced and transportation became more efficient.
According to a simulation, suggested improvements result with greater efficiency and performance of the machines in the workshop. Comparison of the data from table 1 and 2 shows that greater quantity of material could be processed on the cutting machines. This is, very important
Fig. 7. Distribution of material according to the working site in the second experiment
Fig. 8. Profiles cutting robot status during the second experiment
Plazma
Numorex
Corta
Semi-automatic cutting machine
Fig. 9. Status of cutting machines and plate edge preparation bench during the second experiment Table 2. Number offabricated plates, profiles and bars in second experiment
Cutting machine	Plates	Profiles and bars
Robot	-	14756
Plasma	4183	-
Numorex	1923	-
Corta	1142	-
Total	7248	8641
Semi-automatic cutting machine	1638	-
for plasma cutting machine. It could cut about 20% more plates with better transport organization. Improvements influence the quality of the material flows. It is seen that there are not huge quantities of material on work-in-process storage areas, which could cause problems with sorting and control.
The effects of proposed changes in workshops for plates and profiles cutting are positive, so it is necessary to make a technical and economic analysis that would estimate necessary investment and payback period.
3 CONCLUSION
Paper presents the usage of the simulation methods in the shipbuilding production process management. The entire workshop (machines,
tools, working areas, transport devices) and products (plates, bars and profiles) was modeled and analyzed by simulation.
Simulation model of the shipyard workshop was verified and validated, ensuring the confidence in suggested improvements of the technological processes.
Experiments with model are used to analyse different technological solutions for a production process, and to anticipate situations in case of investment into new equipment or when alternative production programs are considered which might affect particular stages in the production process. Differently defined experiment scenaria include changes in model parameters, and carrying out of a number of experiments in a short period of time enables a broadly based research, resulting in a
database that makes it possible to evaluate status of all the production process stages according to the criteria such as efficiency of equipment and devices, process duration or capacity of interim stores. The analysis conclusions are used for efficient management and decision-making process in a real system.
An advantage of simulation procedure use in production process management is that it enables creation of a realistic model of the production process and final product, and an approximation of situations encountered in a real production system to be used in production process control and management. Complex production processes, such shipbuilding, are too complicated to be resolved by mathematical analysis, and simulation methods have proven better for analysis and understanding of their behaviour. Their application is possible in already defined production processes as well as in those, which are in a design stage.
Application of the state-of-the-art program packages enables a more extensive use of simulation procedures in production process control and management than before. This is certainly due to advancement in computer technology and capacities, primarily in their operation speed and graphic software. Advancement in program packages resulted in their being user-friendlier; the users need not involve in programming, coding and testing of programs. Instead, they may focus on modelling, experiments with model and analysis of results.
An advantage the state-of-the-art program packages have is an option of 2D and 3D visualisation of production process model, which facilitates an experiment monitoring, detection of errors in model creation and pinpointing of interruptions or bottlenecks in the production process.
4 REFERENCES
[1] Ljubenkov, B. Simulation in management of the jack-up drilling rig building process, Master Thesis. University of Zagreb, Faculty of Mechanical Engineering and Naval Architecture, 2002. (in Croatian).
[2]	«iljak, V. Simulation by computers. Zagreb: Školska knjiga, 1982. (in Croatian).
[3]	Medeiros, D.J., Traband, M., Tribble, A., Lepro, R., Fast, K., Williams, D. Simulation based design for shipyard manufacturing process. Proceedings of the 2000 Winter Simulation Conference, Orlando, USA.
[4]	Williams, D., Finke, D.A., Medeiros, D.J., Traband, M.T. Discrete simulation development for a proposed shipyard steel processing facility. Proceedings of the 2001 Winter Simulation Conference, Arlington, USA.
[5]	Kiran, A., Cetinkaya, T., Cabrera, J. Hierarchical modelling of a shipyard integrated with an external scheduling application. Proceedings of the 2001 Winter Simulation Conference, Arlington, USA.
[6]	McLean, C., Shao, G. Simulation of shipbuilding operations. Proceedings of the 2001 Winter Simulation Conference, Arlington, USA.
[7]	Dain, O., Ginsberg, M., Keenan, E., Pyle, J., Smith, T., Stoneman, A., Pardoe, I. Stochastic shipyard simulation with SimYard. Proceedings of the 2006 Winter Simulation Conference, Monterey, USA.
[8]	Nedess, C., Friedewald, A., Wagner, L., Hubler, M. Simulation of material flow processes in the planning ofproduction spaces in shipbuilding. On-line: www.tudelft.nl.
[9]	Steinhauer, D., Meyer-König, S. Simulation aided production planning in block assembly. On-line: www.tudelft.nl.
[10]	Sladoljev, Z., Zaplatić,T., Ljubenkov,B. Some possibilities of shipbuilding process management by simulation method. The Tenth Congress of the International Maritime Association of the Mediterranean, Crete, Hellas, 13-17 May, 2002.
[11]	Kuzmanić, M. Technological process of shipbuilding workshop, Graduation project. University of Zagreb, Faculty of Mechanical Engineering and Naval Architecture, 2006. (in Croatian).
Paper received: 1.8.2006 Paper accepted: 28.9.2007
Investigation of Stress Concentration Factor for Finite-Width Orthotropic Rectangular Plates with a Circular Opening Using Three-Dimensional Finite Element Model
Kambiz Bakhshandeh* - Iraj Rajabi - Farhad Rahimi MUT University, Air Naval Research Center, Iran
In this paper, theoretical stress concentration factor for finite-width orthotropic plates with centered circular opening is investigated by use of three-dimensional finite element analysis. In the first part, the finite element (FE) model is validated by comparison analytic equations for infinite orthotropic and isotropic plates. With verified model, the stress concentration factor is calculated for orthotropic rectangular plates with centered circular opening subject to uniform tension. The effect of opening radius to width ofplate ratio for various orthotropic plates with different orthotropy ratia is investigated. This study shows the accuracy of analytical method is dependent on orthotropy ratio and a/W ratio. In high opening radius to plate width ratio, the orthotropy ratio has a little effect on stress concentration factor. © 2008 Journal of Mechanical Engineering. All rights reserved. Keywords: stress concentration factor, finite element methods, orthotropic plates
0 INTRODUCTION
Stress concentrations around cutouts have great practical importance because they are normally the cause of failure. For most materials, the failure strengths of the materials are strongly notch (or hole) sensitive. The net failure stress, taking into account the reduction in cross-section area, is typically much less than the ultimate tensile strength of the same materials without the notch or hole. Until now, the stress concentration factor for various isotropic structures is investigated. Heywood [1], Peterson [2] and Pilkey [3] have investigated various isotropic shape with wide range of holes. Howland [4] determined the solution of the problem of a long isotropic rectangular plate with a centered hole subject to a tension load. Peterson and Heywood introduced various equations for long length and finite width plate with different opening shapes.
Lekhnitskii [5] and Tan [6] introduced various formulations for stress concentration in orthotropic materials. Lekhnitskii derives equations for infinite plate with circular holes. For finite-width orthotropic plate Tan introduces various equations. He derives these equations by use of equilibrium relations, but he did not investigate the effect of orthotropy ratio. The orthotropy ratio is
defined as the ratio of two elastic moduli in both main directions. The analysis presented here suggested, that Tan's equations are applicable only for a specific range of orthotropy ratio.
In this paper, the effect of orthotropy ratio on the stress concentration factor of rectangular finite-width plates with circular opening is investigated. These plates are under action of unidirectional tension loads. In the first part, the finite element model is validated for infinite orthotropic and isotropic rectangular plate with circular opening. With validated model, the effect of orthotropy ratio on stress concentration factor is investigated.
1 THEORETICAL STRESS CONCENTRATION FACTOR FOR FINITE WIDTH PLATE
As is well known the theoretical stress concentration factor (TSCF) for normal stress is defined according to [2] and [7]:
snom
where the stress о represents the maximum stress
max A
to be expected in the member under the actual loads and the normal stress о is reference normal
*Corr. Author's Address: MUT University, Air Naval Research Center, Shiraz, Postal Code: 7194915685, Iran, 140	kambiz@bakhshandeh.ir
stress. In a rectangular plate under action of unidirectional tension load S (Fig. 1), a is
v ^ /? nom
calculated based on the whole cross-section of the plate (gross area) as [2]:
S 2 Wh
(2),
where:
S - traction force, uniform on edges of plate, W - plate half-width, h - plate thickness.
In this way, defined reference stress is in fact the far away field stress. TSCF calculated from Equations (1) and (2) is denoted by KTg and is very useful when investigating wide plates with small openings, i.e. when approaching infinite domain problems.
Some values of the TSCF have been compiled by some major references Peterson, Pilkey, Tan and Lekhnitskii, and are reproduced in a number of mechanical design texts used worldwide.
For a finite width plate, the stress concentration factor can be calculated by use of finite-width correction (FWC) factor [6] and [8]. FWC factor is a scale factor defined as the stress concentration factor of an infinite to a finite-width plate. Thus, the stress concentration of finite-width plate can be calculated with specification of FWC factor and stress concentration of infinite plate with similar opening and material (orthotropic or isotropic).
FWC factor for isotropic material is independent of material properties thus can be determined accurately by use of curve fitting technique. Whereas in orthotropic material, this factor must be obtained with stress analysis by use
Fig. 1. A finite-width plate containing a central elliptical opening
of elasticity equations or mainly using finite element methods [6].
According to investigation of Tan on rectangular orthotropic plate with elliptical opening under action of unidirectional axial load (Fig. 1), the FWC factor is introduced in the form of following equation:
(1-2Л)
K
Tg
Л2
K
Tg (1-Л)2 (1-Л)2
1+(Л2-1)
(1-Л)
W
-M
a ,, — M W
-1/2
1 + (Л2-1)
W
M
a
M
W
K %
1-1
Л
(3),
1 + (Л2-1)
W
M
-5/2
W
M
-7/2
1 + (Л2-1)
W
M
where:
1 - 8
M2 = -
1 - a
W
2 +
/ \3
1 - a
w
-1
-1
V
(4),
Л = a/b,
K%g - stress concentration factor in infinite plate, KTg - stress concentration factor in finite plate.
The magnification factor M, is only a function of a/W and is independent of material properties.
For a circular opening Л = a/b = 1 with using the L'Hospital's rule twice w.r.t Л [6], Equation (3) becomes:
11 a
Tg
.+1 (±M T (k;- 3)1 -faM 2{w }K T ' {W
2 л
KTg 2+(1 - a,
{ W J	(5).
This equation can be used for orthotropic and isotropic finite-width plates with circular opening under action of unidirectional tension load for a/W up to 0.9 [6].
For calculation of stress concentration factor in a finite-width plate besides of FWC
2
7
factor the stress concentration factor of infinite plate with similar opening must be calculated. For orthotropic plate, the stress concentration factor of finite-width plate is complicated with respect to isotropic plates. The Lekhnitskii [5] introduces stress concentration factor for orthotropic infinite plate, Kjgo as below:
KTg,o = — ( [- cos2 j + (k + n) sin2 jJ к cos2 в +
where:
Eli
+ [(1 + n)cos 2 j — к sin 2 j J s in 2 в — —n(1 + к + n) sin j cos j sin в cos в)
(6),
n =
v
pi. — u12 + -Ei-E22 2G12
к = Eii E2;
(7).
1 sin4 в
Ец
2v,-
G-
E11
sin2 в cos2 в +
cos4 в
E and E22 are elasticity moduli in main directions, Gl2 is modulus of rigidity in shear plane, u12 is Poisson ratio, в is polar angle measured from X axis (main direction 1) and ф is the angle of acting force with respect to main direction 1 (Fig. 2).
For an infinite plate with circular opening under action of unidirectional tension, with direction of the applied force aligned with the main direction 1, ф is 0 and polar angle в at points 'c' and 'd' on Figure 2 is 0 and 90 degree, respectively. The stress concentration factors for these two points are:
k;j = 1+n
Tg, o
K7l = -к
(8),
where:
КТО - stress concentration factor of point 'c' for orthotropic infinite plate,
KT О - stress concentration factor of point 'd' for orthotropic infinite plate.
In isotropic material E = E22 and G = E/ 2(l+v), and according to Equations (6) and (7), the stress concentration factor for point 'c' and 'd' are 3 and - l, respectively. Substituting, Kjg = 3 in the Equation (5), it is simplified to:
K;g = 3(1 - a/W)
К
(9).
v g 2 +(1 - a/W)3 This is familiar Heywood equation for isotropic infinite-plate with circular opening. Also, Howland determined the solution of the problem of a long rectangular plate with a central hole subject to tension in an open series form [4].
Another way to approaching the stress concentration factor and FWC factor is the application of the average normal stress across the net section (2(W-a)h) instead of the whole width of the plate (gross area). This stress concentration factor is known as net stress concentration factor, KTn and can be related to the previous stress concentration factor, K by:
KTn = Kg(1 - a/W)
(10).
In the common engineering sense it is more convenient to give stress concentration factors using reference stress based on the net area rather than the gross area [9]. Thus in this study the net stress concentration factor is used.
1.1 Finite Element Model Validation
The stress concentration of an isotropic and orthotropic plate with circular opening under action of unidirectional tension is calculated by use of a finite element software MSC.Nastran [10].
Due to symmetry and using the appropriate displacement boundary conditions only one quarter of the plate was used as occupied by the shaded area indicated in Figure 3. The utilized restrictions
Fig. 2. An infinite-plate containinga central circular opening
Fig. 3. Finite element model region
are zero displacement in the X direction for all the mesh nodal points located on the line X = 0, and zero displacement in the Y direction for all mesh nodal points, located on the line Y = 0.
The finite element model is discretized with twenty-node solid brick elements (Fig. 4).
Various mesh densities for each of plates are used to achieve a good mesh density for improving accuracy of results. The mesh density is changed from coarser to finer. On the average and roughly the fine mesh have ten times more nodes and elements than the coarse mesh. Densifying is iteratively repeated until the difference between the consecutive results becomes smaller than 0.35%. In Figure 5 the final discrepancy of the TSCF for various opening ratia is shown. This figure reflects the absolute value of the difference (percent) in the magnitude of the TSCF due to the mesh refinement. Between 67171 to 70501 elements are used in the study.
Because the isotropic material properties are not different in main directions, this is used for
Fig. 4. A typical finite element model
improving the finite element model accuracy besides of orthotropic material.
The finite element model is created for various opening ratia, on plates with dimensions: half-length L = 1000 mm, half width W = 100 mm and thickness 10 mm for isotropic and orthotopic material. The opening radius to half-length of plate ratio a/W is changed from 0.1 to 0.8.
The isotropic plate has elastic modulus E = 210 GPa, Poisson ratio v = 0.3 and mass density p = 7800 kg/m3.
The net stress concentration factor is calculated for isotropic plates by use of Equations (9) and (10) and compared with FE method results.
Table 1 shows the net stress concentration factors that are concluded by FE and analytical method for isotropic plates.
The error e in this table is defined as:
TSCFjheory — TSCFfem
TSCF
Theory
100
where TSCFTheory and TSCFfem are, respectively, the
E	13
1	!
X	E
га	a
2	.E
Opening radius to plate width [a/W]
Fig. 5. Convergence study results
Table 1. Net stress concentration factor for isotropic plate
	Net Stress Concentration Factor			EEOR***
a/W	Heywood*	Howland*	FE	e %
0.1	2.72	2.72	2.73	0.37
0.2	2.61	2.62	2.59	0.61
0.3	2.34	2.36	2.35	0.31
0.4	2.22	2.24	2.23	0.35
0.5	2.12	2.16	2.11	0.41
0.6	2.06	2.09	2.05	0.29
0.7	2.03	2.05	2.04	0.24
0.8	2.01	2.03	2.02	0.40
* Eq. (9) **Ref. [4]
*** Error (Eq. 11) of FE results with respect to Heywood results
stress concentration factors calculated with the analytic and FE method.
As it can be seen from Table 1 the agreement between the finite element method and theoretical solution for isotropic plate is excellent.
The distribution of net stress concentration factor in a vicinity of circular opening for isotropic plate with a/W = 0.1 is shown in Figure 6.
For orthotopic plate KTn is also calculated by use of FE and analytical method (Eqs. (5), (8) and (10)). The properties of the orthotopic material are:
Corresponding to a typical orthotropic material [11].
In Table 2, the calculated net stress concentration factor is compared. The error that is shown in this table is calculated by Equation (11).
Similar to isotropic plate, the FE method has a good accuracy for orthotropic materials.
The distribution of net stress concentration factor in a vicinity of circular opening for orthotropic plate with a/W = 0.1 is shown in Figure 7 .
Ell = 44.7 GPa, G23 = 3.45 GPa,
= 0.34,
E22 = E33 = 17.9 GPa, v12 = 0.25,
9 = 0°,
G12 = G13 = 8.96 GPa,
v13 = 0.25
Fig. 6. Net stress concentration factor distribution of isotropic plate with a/W = 0.1
Table 2. Net stress concentration factor for orthotopic plate
	Net Stress Concentration Factor		ERROR**
a/W	Theoretical*	FE	e %
0.1	3.43	3.42	0.29
0.2	3.15	3.16	0.44
0.3	2.94	2.93	0.32
0.4	2.76	2.77	0.36
0.5	2.61	2.62	0.38
0.6	2.47	2.48	0.39
0.7	2.34	2.33	0.43
0.8	2.22	2.21	0.45
* Eqs. (5), (8), (10) ** Eq. (11)
Fig. 7. Net stress concentration factor distribution of orthotropic plate with a/W = 0.1
2 ORTHOTROPIC RATIO EFFECT
INVESTIGATION PROCEDURE
In the previous section, the FE method was validated for an orthotropic and isotropic finite-width plate with different circular opening radius to plate width ratio. In this section, the effect of orthotropy ratio on the stress concentration factor is investigated. The orthotropy ratio Eu/E22 is defined as the ratio of elastic modulii in main directions.
The procedure employed in the research consisted in determining stress concentration factor for a group of plates with a constant length to width ratio and different a/W and orthotropy ratio. The a/W ratio is changed from 0.1 to 0.8 and for each ratio the stress concentration factor is calculated for
orthotropy ratio ranged from 2 to 50. In this research, the elastic modulus in width direction (E22) is assumed to be constant and has a value of 17.9 GPa.
3 RESULTS AND DISCUSSION
The stress concentration factor for various a/W ratio and orthotropy ratio is calculated by FE method and compared with analytical results (5), (8) and (10). The absolute difference of FE method and analytical method is calculated by Equation 11 and shown in Figure 8.
Tan mentioned that his Equation (5) could be used for wide range of a/W ratio even up to 0.9. According to Figure 8 this equation's accuracy is not good for all range of a/W ratio. As shown in Figure 8 the analytical and FE results for a/W < 0.4
Orthotropy ratio [EnIEì2]
Fig. 8. Difference of analytic and FE method results
have a little difference for various orthotropy ratio (difference is less than 4%), but with increasing the а/W ratio this difference increases. In low orthotropy ratio the difference between analytical and FE results have a little difference, but this difference increases with increasing orthotropy ratio. Thus, the accuracy of Tan's Equation (5) decreases with increasing the а/W ratio and orthotropy ratio.
The variation of the net stress concentration factor with respects to а/W ratia for various orthotropy ratios is shown in Figure 9. According to this figure the higher orthotropy ratio have a higher stress concentration factor. In addition, the variation of SCF respect to а/W ratio in high orthotropy ratio is higher than at low orthotropy ratio. The important observations regarding the

Fig. 9. Net stress concentration factor respect to orthotropy ratio and a/W ratio
graph of this figure are that all curves converge together with increasing the а/W ratio. Indeed, with increasing the а/W ratio, the effect of orthotropy ratio on net stress concentration factor is decreased.
According to Figure 9, the net stress concentration factor for low orthotropy ratio is near to each other and increasing this ratio the difference of а/W curves. The various а/W curves are converged to each other with low Eu/E22 ratio. In addition, with increasing the а/W ratio, the variation of orthotropy ratio has a little effect on the net stress concentration factor. This effect is also shown in Figure 9.
4 CONCLUSION
Formulation and comparative finite element results have been presented for a conforming finite element method, which may be used to model a plate with circular opening under action of tension load, and this comparison shows that appropriate finite element model has a good accuracy for stress concentration problems.
The stress concentration factor for low orthotropy ratio is near to each other and with increasing this ratio, the difference of а/W curves increases. The variation of net stress concentration factor, respect to orthotropy ratio increases with increasing the а/W ratio. In high а/W ratio the effect of orthotropic ratio is decreased, thus in high а/W ratio the elastic modulus changing in applied force direction has a little effect on the stress concentration factor. In other word: high orthotropy ratio and low opening ratio (а/W) increase net stress concentration factor.
The accuracy of analytical method is dependent on orthotropy ratio and a/W ratio and it is suitable for a/W < 0.4. Increasing the a/W ratio, this accuracy decreases further the influence of orthotropy ratio is small. At a/W ratio 0.8, discrepancy between analytical and numerical method is less than 6% , which is acceptable in engineering applications. Therefore, Tan's Equation (5) can be used instead of tedious and time consuming finite element analyses for all investigated ratia (up to 0.8).
5 REFERENCES
[1]	Heywood, R.B. Designing by photoelasticity. Chapman and Hall, 1952.
[2]	Peterson, R.E. Stress concentration factors. John Wiley & Sons, 1974.
[3]	Pilkey, W.D. Peterson s stress concentration factors. New York: John Wiley & Sons Inc., 2nd Edition, 1997.
[4]	Howland, R.C.J. On the stresses in the neighborhood of circular hole in a strip under tension. Phil. Trans. Roj. Soc.,Vol. 229, 1929, p. 49-86.
[5]	Lekhnitskii, S.G. Anisotropic plates. Gordon and Breach, 1968.
[6]	Tan, S.C. Finite-width correction factors for anisotropic plate containing a central opening.
J. of Composite Materials, Vol. 22, 1988, p. 1080-18.
[7]	Wang, Q.Z. Simple formulae for the stress concentration factor for two - and three-Dimensional Holes in Finite Domain. J. of strain analysis, 37, 2002, p. 259-264.
[8]	Tan, S.C. Laminated composites containing an elliptical opening. II. Experiment and Model Modification. J. of Composite Materials, Vol. 21, 1987, p. 949-20.
[9]	Hwai-Chung,W., Bin, M. On stress concentrations for isotropic/orthotropic plates and cylinders with circular hole. J. of Composites, Part B: engineering, 34, 2003, p. 127-134.
[10]	Manual: MSC.Nastran for Windows: Quick Start Guide. MacNeal-Schwendler Corporation, 2004.
[11]	Reddy, J.N. Mechanics oflaminated composite plates and shells: Theory and analysis. 2nd Edition. New York: CRC Press, 2004.
This paper presents the development of a robot that mimics the movement ofa live snake. A prototype robot comprising six links was constructed. Torque actuators between the links modify the robot's shape. Anisotropic friction between the links and the ground generates the force that propels the robot. A control variable that determines actuator angles is used to achieve a wave-like body motion. The corresponding signal is transmitted over a radio link. Measurements of average velocity and trajectory of the robot were performed with different control parameters. Basic properties of the robot's movement are presented. © 2008 Journal of Mechanical Engineering. All rights reserved. Keywords: snake-like robots, control, serpentine windings, friction forces
Strojniški vestnik - Journal of Mechanical Engineering 54(2008)2, 148-153	Paper received: 12.3.2007
UDK - UDC 007.52	Paper accepted: 27.6.2007
Development of a Snake-like Robot
Jure Bezgovšek* - Igor Grabec - Peter Mužic - Edvard Govekar University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
0 INTRODUCTION
In relation to the development of modern endoscopes, there is increasing interest in snake-like robots achieving motion in a desired direction by undulation of the body as shown in Fig. 1. Such robots make feasible new solutions to complex problems, such as the non-destructive testing of mechanical parts in technical devices, flexible machining of geometrically complicated cavities and various surgical operations in medicine [1].
0.1 Characteristics of Snake-like Robots
A limited number of degrees of freedom in existing endoscopes lead to restrictions in their application where movement in a geometrically complex space is required. Consequently, the development of a robot mechanism with many of degrees of freedom has been initiated recently in order to enable movement in diverse environments. Imitation of biological organisms with flexible bodies, such as highly developed snakes, is a possible course in the development of such robot mechanisms.
Snakes can be found in diverse natural habitats [2]: they live in deserts, tropical forests; they are capable of swimming in rivers and oceans. During evolution snakes have lost their limbs and developed an extended spine in order to facilitate movement in constrained environments, such as narrow underground tunnels, or grass and bush terrain. These characteristics lead us to imitate snakes in the development of a robot designed for performance that requires high flexibility.
The body of a snake consists of a spinal column, muscles and a sensory-neural network. The spinal column is a passive element that transmits forces and determines the dimensions of the body. The same property is achieved by the mechanical structure of a snake-like robot that includes links and that permits variation of the robot shape. Muscles can be replaced by controllable torque actuators causing undulation of the robot. The actuators should be positioned at the joints in order to change the angles between adjacent links. Adaptation of a body to complex forms in the environment and creeping locomotion is controlled in snakes by a complicated sensory-neural network. In our development of this snake-like robot, this task is performed by the controller of the torque actuators.
Our goal was to develop an appropriate control strategy by which desired robot propulsion could be achieved. The structure of a snake-like robot is presented in Figure 2.
0.2 Propulsion Principle
Different types of snake locomotion for various purposes exist [2]. For example, one type of locomotion is suitable for soft terrain, while another works on hard terrain. Locomotion in narrow channels is also possible. Well-known locomotion types are serpentine, rectilinear, lateral rolling, concertina, etc. [3] and [4]. Our research is focused on serpentine locomotion, because it appears to be the simplest. Its main characteristic is the wave-like undulation of the snake body. To generate propulsion with this kind of movement, adequate friction between the body and
*Corr. Author's Address: University of Ljubljana, Faculty of Mechanical Engineering, Aškerčeva 6, SI-1000 Ljubljana, 148	Slovenia, jure.bezgovsek@fs.uni-lj.si
\Л
Fig. 1. Undulating movement
the ground must exist. The friction should be anisotropic with regard to the direction of velocity and orientation of a particular link [1]. This means that the transverse frictional force on a particular link should be much greater than in the direction of the link. The basic problem is then to find the undulation of the body that would result in the movement of the robot along a prescribed path.
1 CONSTRUCTION OF THE SNAKE-LIKE ROBOT
The aforementioned propulsion principle is applicable to achieve robotic movement in two or three dimensions. In our research we have considered just two-dimensional cases on a flat plane. Hence the joints axes between the links are parallel to each other and normal to the plane. The robot was composed of six links connected by hinges (Fig. 3). Between the links there are five torque actuators.
For our research it is not necessary for the robot to withstand high mechanical loads and the path following tolerance is not very strict. Therefore, the robot needs not to be built of highly accurate and rigid parts with tight tolerances. Consequently, we have assembled the robot from simple flexible plastic
Fig. 3. Snake-like robot made up of six links connected by hinges
elements [5]. These elements also enabled the simple changing and improvement of the structure. The final version of the robot is presented in Figures 4 and 5.
As torque actuators we used DC electric servo motors, of the type usually used in remote control model aircraft. These motors are convenient for developing robot prototypes because they have a built-in feedback loop that simplifies control of their position angles. Furthermore, these motors are designed to work with radio control and receiving units permitting remote control. The servo motors are positioned on the turning axis of the joints between adjacent links. The motors do not need any additional gears, since they generate enough torque for the change of the angles between links. The motors receive control signals from the radio receiver unit placed on the head of the robot.
To assure anisotropic friction, we need suitable contact of the robot with the ground. Different materials have already been tested for this purpose by other researchers [2]. Passive wheels have been found to be the most appropriate, since they have a very high coefficient of friction in the transverse direction and almost no friction in the direction of rolling. Consequently, we have also utilized passive wheels. Initially, we used plastic wheels, but these
Fig. 4. Final version of a snake-like robot
Receiver
Servo motors
Fig. 5. Drawing of a complete experimental snake-like robot
ф5 (t ) = R%2 • sin (2nv • t -(s - 1)ф) + Rtr (2).
lacked sufficient transverse friction for good propulsion. To increase friction, we replaced the plastic wheels with rubber wheels.
As a power supply for the radio receiver and servo motors we utilized rechargeable electrical accumulators. It is important to achieve uniform and sufficient friction along the robot. For this reason a pair of accumulators is placed on each link as extra weight.
2 CONTROL OF THE SNAKE-LIKE ROBOT
The propulsion of the robot results from friction of the wheels due to undulation of the body on the surface. This undulation is achieved by controlling the angles between links by servo motors [4]. A control variable ф(И) determines the values of the angles between the links. To control the movement of the robot along a prescribed planar path this control variable was determined by a sum of two control signals. The first component causes undulation of the robot body, while the second one causes an average bending that corresponds to the curvature of the prescribed path. The corresponding control variable ф/t) is described by the following function:
ф5 (t) = R • sin (• t-(s - 1)ф) + Rtr (1).
Here R, v, ф and Rr are settable parameters. R is the relative amplitude of oscillation of the servo motors with regard to the maximum possible amplitude, v is the frequency of oscillation of the servo motors and ф is the phase delay between servo motors. Rr is the parameter that determines the average bending of the robot body and leads to steering along the prescribed path as shown in Figure 6. Index s denotes the index of the link and tr denotes path curvature.
In our testing, two control strategies were examined: in the first the relative amplitude of oscillation was kept constant (1), while in the second it was modulated with the position of the link [1]:
Here ^=s/S is the normalized index of a link, where S is the number of links, and additional parameter describes the modulation of the undulation amplitude that increases along the body from the head towards the tail of the robot.
The robot is controlled by a personal computer. In our case the control signal is transmitted from the computer to the servo motors over a radio unit designed for remote control of model aircraft [6]. Since the robot has no direct mechanical connection to the controller, the only restriction regarding its movement is the range of the radio signal. The range greatly depends on the ambient obstacles (walls and cabinets in the laboratory), but in open space it should be at least 80 m.
The structure of this control system is presented in the diagram in Fig. 7. This is an open-loop, realtime control. The user enters control parameters either through a keyboard or joystick to the control algorithm that runs on a personal computer, which generates the required control signal in real-time.
Mechanism of the robot
Prescribed path
Fig. 6. Scheme of a snake-like robot tracing a prescribed path
Radio transmission of control signal to the robot
----- Robot
Fig. 7. Flow diagram of control system
The computer is connected to the radio control unit through the parallel port [7]. The radio unit transmits the signal using modulation of high frequency radio waves. The radio receiver on the robot extracts the control signal from the received radio wave and transforms it to match the characteristics of the servo motors. Then it routes the correct signals to particular motors. A weakness of this type of control is that it does not provide for a feedback connection by which the controller could receive information about the state of the robot. We assume that using a feedback would lead to increased robot movement accuracy.
3.1	Movement on Various Surfaces
Properties of the ground surface directly affect the movement of the robot. They are described by the coefficient of friction k and the roughness R . The higher the coefficient of friction between the surface and the passive wheels, the more effective is the creeping of the robot. The best conditions are met when there is no slip of the wheels in the transverse direction.
To analyze the influence of surface on the effectiveness of creeping we conducted a set of experiments on surfaces of linoleum (Ra ~ 0.5 |im), spongy plastic (Ra ~ 10 |im) and smooth plastic (Ra ~ 0.1 |im). We found that propulsion of the robot is most effective on the spongy plastic material, which has the highest roughness Ra. The results on linoleum were also satisfactory, while creeping on smooth plastic was rather ineffective.
3.2	Measurements of Average Velocity
Measurements of the average velocity of the robot were made on the linoleum surface. We adapted the non-modulated control variable (1), where four control parameters R, v, p and Rr are present. For velocity measurements we only used creeping in a straight line. In this case Rtr = 0 and the control variable is:
ф (t) = R■ sin(inv ■ t-(s- 1)p) (3).
We measured the average velocity at seven different combinations of the remaining three parameters R, v and ( (Table 1). The highest velocity of approximately 0.26 m/s was observed at the following combination of parameter values: R = 100%, v = 0.8 Hz and p = 1 rad. Based on the combinations of parameters used, we can roughly estimate the influence of control parameters on movement. Results of the measurements, given in Table 2, indicate that in the tested range of parameters velocity increases with increasing amplitude and increasing frequency. The velocity also increases with the phase delay increasing from p = 1 rad (Fig. 8) to p= 1.5 rad (Fig. 9).
3 EXPERIMENTAL RESULTS
3.3 Measurements of Head Trajectory
Measurements of average velocity and head trajectory were made on the robot creeping on various surfaces.
Measurements of the head trajectory of the robot were also made on the linoleum surface. The trajectories were recorded from the vertical
Table 1. Combinations of control parameters	Table 2. Results of velocity measurements
used in average velocity measurements
Parameter set	R	V	ф
no.	[%]	[Hz]	[rad]
1	100	0.4	1.5
2	100	0.6	1.5
3	100	0.8	1.5
4	100	0.6	1
5	60	0.6	1
6	60	0.6	1.5
7	100	0.8	1
Rank	v [m/s]	Parameter set no.
1.	0.258	7
2.	0.202	4
3.	0.197	5
4.	0.176	3
5.	0.126	2
6.	0.099	6
7.	0.088	1
Fig. 8. Shape of a snake-like robot at phase delay ф = 1 rad
perspective by a digital camera. During experiments we applied the same non-modulated control variable (3) as for the velocity measurements and the same seven combinations of control parameters.
Figure 10 shows the head trajectory for three different combinations of control parameters, including the parameters at which the highest velocity was achieved. The trajectories are rather smooth and resemble a sine curve determined by the control of the body undulation. Since there are no peaks in the curves we conclude that there is no abrupt slipping of the body. Beside this, we observed rather good tracing of a straight line. Similar observations have been found with following a curved prescribed path, when the radius of curvature is much larger than the length of a single link.
4 CONCLUSIONS
We have developed a snake-like robot capable of creeping along a prescribed planar path on various
Fig. 9. Shape of a snake-like robot at phase delay ф = 1.5 rad
surfaces. The robot is controlled by a personal computer, where a user enters control parameters by keyboard or joystick. Application of a joystick offers convenient and continuous manual control of the robot's movement.
In our experiments, velocity and head trajectory of the robot were measured. We tested seven different sets of control parameters and roughly estimated the influence of the control parameters on the movement. The highest achieved velocity of the robot was approximately 0.26 m/s. A surface with a high coefficient of friction has proven to be the most suitable for effective propulsion.
At this stage we have only used open-loop control. The next research step is to upgrade the open-loop into a closed-loop control, where the controller constantly receives feedback information from the robot. The feedback information can contain position coordinates, indication of obstacles, etc. By applying the closed-loop control of the robot, increased tracking accuracy and reliability of the robot can be expected.
Parameter sest no. 1
Parameter set no. 4
Parameter set no. 7
x [cm]
Fig. 10. Head trajectories measured at different combinations of control parameters. At parameter combination 7, the robot achieved the highest velocity.
Some work has already been done on a development of a control system with a feedback loop that includes an infrared optical tracking sensor, where we achieved some promising preliminary results. Our long term goals are to fully implement closed-loop control and to utilize an artificial neural network for autonomous optimization of the robot movement.
5 REFERENCES
[1]Grabec,	I. Control of a creeping snake-like robot. 7th Int. Workshop on Advanced Motion Control, Maribor, 2002, Proc. IEEE Catalog Number: 02TH8623, p. 526-531.
[2]	Miller, S.P.G. Snake robots for search and rescue, Neurotechnology for Biomimetic Robots, edited by J. Ayers, J. L. Davis, A. Rudolph. The MIT Press, 2002.
[3]	Mori, M., Hirose, S. Three-dimensional serpentine motion and lateral rolling by Active Cord Mechanism ACM-R3. Proc. 2002 IEEE/ RSJ Intl. Conf. on Intelligent Robots and Systems. EPFL, Lausanne, Switzerland, 2002.
[4]	Ma, S., Araya, H., Li, L. Development of a creeping locomotion snake robot. International Journal of Robotics and Automation, Vol. 17, No. 4, 2002, p. 146-153.
[5]	Bezgovšek J. Development of a snake-like robot. Degree project, University of Ljubljana, Faculty of Mechanical Engineering, 2005. (In Slovenian).
[6]	PC-to-R/C Interface. On-line: http:// www.mh.ttu.ee/risto/rc/electronics/pctorc.htm
[7]	Oberoi, D. Centronics port generates narrow pulse widths, EDN, 2.8.2001.
Corrigendum
Determination of the Friction Coefficient of Groove Forms Using the Stress-Function Method
C. Erdem Imrak - Ismail Gerdemeli
Strojniški vestnik - Journal of Mechanical Engineering, vol. 53 (2007), no. 1, p.13-17.
As was pointed out to us, the graph plot errors exist in Figure 3. The values of the shape factor for the coefficient of friction should be corrected as illustrated below. The other results of the paper are unchanged.
Instructions for Authors
From 2008 the Journal of Mechanical Engineering is to be published in English only, but with separate Slovene abstracts. The authors are entirely responsible for the correctness of the language. If a reviewer indicates that the language in the paper is poor, the editor will require the author to correct the text with the help of a native English speaker before the paper is reviewed again.
Papers can be submitted electronically to the journal's e-mail address (info@sv-jme.eu) or by post.
Papers submitted for publication should comprise the following:
-	Title, Abstract, Keywords,
-	Main body of text,
-	Tables and Figures (graphs, drawings or photographs) with captions,
-	List of References,
-	Information about the authors, the corresponding author and a full set of addresses.
For papers from abroad (in the case that none of the authors is a Slovene) the editor will obtain a Slovenian translation of the Abstract.
Papers should be short, about 8 to 12 pages of A4 format, or at most, 7000 words. Longer papers will be accepted if there is a special reason, which should be stated by the author in the accompanying letter. Short papers should be limited to less than 3000 words.
THE FORMAT OF THE PAPER
The paper should be written in the following format:
-	A Title, which adequately describes the content of the paper.
-	An Abstract, which should not exceed 250 words. The Abstract should state the principal objectives and the scope of the investigation, the methodology employed, summarize the results and state the principal conclusions.
-	An Introduction, which should provide a review of recent literature and sufficient background information to allow the results of the paper to be understood and evaluated.
-	A Theory and the experimental methods used.
-	An Experimental section, which should provide details of the experimental set-up and the
methods used for obtaining the results.
-	A Results section, which should clearly and concisely present the data using figures and tables where appropriate.
-	A Discussion section, which should describe the relationships and generalisations shown by the results and discuss the significance of the results, making comparisons with previously published work.
-	Because of the nature of some studies it may be appropriate to combine the Results and Discussion sections into a single section to improve the clarity and make it easier for the reader.
-	Conclusions, which should present one or more conclusions that have been drawn from the results and subsequent discussion and do not duplicate the Abstract.
-	References, which must be numbered consecutively in the text using square brackets [1] and collected together in a reference list at the end of the paper.
THE LAYOUT OF THE TEXT
Texts should be written in Microsoft Word format. The paper must be submitted in an electronic version, by e-mail or by post on a CD.
Do not use a LaTeX text editor, since this is not compatible with the publishing procedure of the Journal of Mechanical Engineering.
Equations should be on a separate line in the main body of the text and marked on the right-hand side of the page with numbers in round brackets.
Units and abbreviations
Only standard SI symbols and abbreviations should be used in the text, tables and figures. 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. m/s, K, min, mm, etc.).
All abbreviations should be spelt out in full on first appearance, e.g., variable time geometry (VTG).
The meaning of symbols and units belonging to symbols should be explained in each
case or quoted in special table at the end of the paper before the References.
Figures
Figures must be cited in consecutive numerical order in the text and referred to in both the text and the caption as Fig. 1, Fig. 2, etc. Pictures should be saved in a resolution good enough for printing, and in any common format, e.g., BMP, GIF or JPG. However, graphs and line drawings should be prepared as vector images, e.g., CDR, AI.
All Figures should be prepared in black and white, without borders and on a white background. All the figures should be sent separately in their original formats.
When labeling axes, physical quantities, e.g., t, v, m, etc., should be used whenever possible. Multi-curve graphs should have the individual curves marked with a symbol. The meaning of the symbol should be explained in the figure caption.
In the case that the author wishes, for whatever reason, to publish Figures in colour, the author must pay the resulting costs.
Tables
Tables must be cited in consecutive numerical order in the text and referred to in both the text and the caption as Table 1, Table 2, etc. In addition to the physical quantity, e.g., t (in italics), units (normal text), should be added in square brackets. Each column should have the title line. Tables should not duplicate information that is already noted anywhere in the paper.
Acknowledgement
An acknowledgement for cooperation or help can be included before the References. The author should state the name of the research (co)financer.
The list of references
All references should be collected at the end of the paper in the following styles for journals, proceedings and books, respectively: [1] Wagner, A., Bajsić, I., Fajdiga, M. Measurement of the surface-temperature field in a fog lamp
using resistance-based temperature detectors. Strojniški vestnik - Journal of Mechanical Engineering, February 2004, vol. 50, no. 2, p. 72-79.
[2]	Boguslawski L. Influence of pressure fluctuations distribution on local heat transfer on flat surface impinged by turbulent free jet. Proceedings of International Thermal Science SeminarII, Bled, June 13.-16., 2004.
[3]	Muhs, D., et al. Roloff/Matek mechanical parts, 16th ed. Wiesbaden: Vieweg Verlag, 2003. 791 p. (In German). ISBN 3-528-07028-5
ACCEPTANCE OF PAPERS AND COPYRIGHT
The Editorial Board reserves the right to decide whether a paper is acceptable for publication, obtain professional reviews for submitted papers, and if necessary, require changes to the content, length or language.
The corresponding author must, in the name of all authors, also enclose a written statement that the paper is original unpublished work, and not under consideration for publication elsewhere.
On publication, copyright of the paper shall pass to the Journal of Mechanical Engineering. The JME must be stated as a source in all later publications.
Submitted materials will not be returned to the author. Unpublished materials are not preserved and will not be sent anywhere without the author's consent.
The paper, prepared for publication, will be sent to the author in PDF format. The author should check for any necessary corrections, which should be the minimum required. With this the author confirms that the paper is ready for publication.
PUBLICATION FEE
For all papers the authors will be asked to pay a publication fee prior to the paper appearing in the journal. However, this fee only needs to be paid after the paper is accepted for publication by the Editorial Board. The fee is €180.00 (for all papers with a maximum of 6 pages), €220.00 (for all papers with a maximum of 10 pages) and €20.00 for each additional page. The publication fee includes 25 off-prints of each paper, which will be sent to the corresponding author.
Vsebina
Strojniški vestnik - Journal of Mechanical Engineering letnik 54, (2008), številka 2 Ljubljana, februar 2008 ISSN 0039-2480
Izhaja mesečno
Uvodnik
Alujevič, A.	SI 16
Povzetki razprav
Rihtaršič, J., Subelj, M., Hočevar, M., Duhovnik, J.: Analiza toka skozi centrifugalno puhalo sesalne
enote	SI 17
Čas, J., Klobučar, R., Hercog, D., Safarič, R.: Daljinsko vodenje SCARA robota na osnovi nevronskega
krmilja	SI 18
Imrak, C.E.: Uporaba umetnih nevronskih mrež pri dvojnem/trojnem skupinskem sistemu krmiljenja SI 19 dvigal
Lemeš, S., Zaimović-Uzunović, N.: Uporaba analize uklona za napoved valovitosti pri stopnjevanem
obdelovanju pločevine	SI 20
Dašić, P., Natsis, A., Petropoulos, G.: Model zanesljivosti rezalnih orodij: Primer v proizvodnem
inženirstvu in agrotehniki	SI 21
Ljubenkov, B., Đukić, G., Kuzmanić, M.: Simulacijske metode pri načrtovanju postopka gradnje ladij SI 22 Bakhshandeh, K., Rajabi, I., Rahimi, F.: Raziskave faktorja koncentracije napetosti ortotropnih pravokotnih plošč s končno širino ter krožno odprtino z uporabo trirazsežnega modela končnih elementov	SI 23
Bezgovšek, J., Grabec, I., Mužič, P., Govekar, E.: Razvoj kačastega robota	SI 24
Poročila
Eksergija: Staro je novo	SI 25
Osebne vesti
Doktorat, magisteriji in diplome
SI 26
Uvodnik
V	tem zvezku našega Strojniškega vestnika smo pripravili osem znanstvenih člankov, od tega je enega sprejel Izdajateljski svet, ostale pa še Uredniški odbor v prejšnji širši sestavi leta 2007. Trenutno imamo na zalogi precejšnje število že sprejetih prispevkov, še več pa jih je v recenzijskem postopku na tujem. Če bomo uspeli objavljati v zastavljenem obsegu 8 člankov v vsaki od naslednjih letošnjih številk, to enostavno pomeni, da vse prispevke 2008 lahko uvrstimo šele v 55 letnik. Ker bo naslednje leto izjemno s 13 polnimi lunami (31. decembra. 2009), ko se bom upokojil kot univerzitetni profesor, to enostavno pomeni, da bo moral moj novi uredniški namestnik profesor Vincenc Butala prevzeti opravila v čim krajšem možnem času.
V	letu 2008 se prav tako spominjamo stoletnice rojstva zaslužnega profesorja Bojana Krauta, ki je bil ustanovitelj in prvi urednik sedaj
našega Strojniškega vestnika, tam v letu 1955, zelo znan pa je tudi po svojem žepnem Strojniškem priročniku (prvič izdanem že 1954 pri TZ Litostroju). Profesor B. Kraut je bil rojen v Kamniku 12. aprila 1908, diplomiral pa je 1932 na TF v Zagrebu. V letih 1936 do 1946je deloval Tovarni železniških vagonov Slavonski Brod. Od 1948 do 1950 je bil tehnični direktor Litostroja, istočasno pa tudi že izredni profesor na Oddelku za strojništvo v Ljubljani, kjer je na novi Fakulteti za strojništvo leta 1962 postal redni profesor za tehnologijo kovin in tirna vozila, do upokojitve 1972. Leta 1978 mu je Univerza v Ljubljani podelila naziv zaslužnega profesorja, Univerza v Mariboru pa častni doktorat leta 1984. Umrl je 22. avgusta 1991 v starosti 83 let. Strojniški vestnik je le ena od njegovih duhovnih dediščin.
Urednik Prof. dr. Andro Alujevič
Prof. Bojan Kraut
Prejeto: 3.1.2008 Sprejeto: 12.3.2008
Analiza toka skozi centrifugalno puhalo sesalne enote
Janez Rihtaršič* - Matjaž Subelj - Marko Hočevar - Jože Duhovnik Univerza v Ljubljani, Fakulteta za strojništvo
V prispevku je z eksperimentom in numerično simulacijo prikazana vizualizacija toka tekočine skozi centrifugalno puhalo sesalne enote. Za eksperiment sta bila izbrana dva obratovalna režima.
Eksperimentalna vizualizacija toka je bila izvedena z dodajanjem pasivnih zrn tekočine. Z uporabo teorije dinamične podobnosto smo znižali vrtilno frekvenco in namesto zraka je bila kot delovno sredstvo uporabljena voda.
Primerjalne numerične simulacije smo izvedli z uporabo stisljive in nestisljive tekočine. Pri eksperimentu je tok tekočine skozi kanale puhala predstavljen s potmi pasivnih zrn in z relativnimi hitrostmi, medtem ko je pri numerični simulaciji predstavljen s tokovnicami in z relativnimi hitrostmi. © 2008 Strojniški vestnik. Vse pravice pridržane.
Ključne besede: vizualizacija toka, centrifugalni rotorji, dinamična podobnost, numerične simulacije
Sl. 1. Mrežni model puhala difuzorja in povratnih kanalov sesalne enote.
*Naslov odgovornega avtorja: Univerza v Ljubljani, Fakulteta za strojništvo, Aškerčeva 6, 1000 Ljubljana, janez.rihtarsic@fs.uni-lj.si
Prejeto: 20.3.2007 Sprejeto: 19.12.2007
Daljinsko vodenje SCARA robota na osnovi nevronskega
krmilja
Jure Čas* - Rok Klobučar - Darko Hercog - Riko Safarič Univerza v Mariboru, Fakulteta za elektrotehniko, računalništvo in informatiko
Opisana je izdelava eksperimenta za daljinsko vodenje SCARA robota (Selective Compliance Assembly Robot Arm - SCARA) z uporabo krmilja, ki deluje na osnovi nevronske mreže. Daljinsko vodenje SCARA robota je izvedeno preko medmrežne povezave. Izdelana uporaba se na Univerzi v Mariboru kot oddaljen preizkus uporablja v studijske namene. Osnovana je na programskih orodjih kot sta MATLAB/ Simulink in LabVIEW. MATLAB/Simulink je skupaj s knjižnico DSP2 Library for Simulink uporabljen za razvoj krmiljnega algoritma in generiranje izvrsne kode. Izvrsna koda se s pomočjo prej omenjene knjižnice naloži in izvaja na digitalnem signalnem procesorju (DSP), ki preko svojih analognih in digitalnih vhodov ter izhodov izvaja krmiljenje SCARA robota in omogoča podatkovno povezavo z laboratorijskim strežnikom. Virtualni instrument (VI), ki je izdelan v programu LabVIEW, je uporabljen kot uporabniški vmesnik, ki s tehnologijo "Remote Panels " omogoča nastavljanje parametrov krmilja ter prikaz želenih podatkov v obliki numeričnih kazalnikov in grafov tudi v medmrežnem brskljalniku. Glavna prednost nevronskega krmilja je zmožnost samostojnega učenja dinamike robota mehanizma in s tem posledično sposobnost natančnejše regulacije. V primeru, ko se pojavi dodatno trenje ali nepričakovana ovira, uporabnik oddaljenega preizkusa ne potrebuje podatkov o nastalih spremembah dinamike robota, ker so spremembe dinamike ocenjene neodvisno od uporabnika preko algoritma za učenje nevronske mreže. © 2008 Strojniški vestnik. Vse pravice pridržane.
Ključne besede: robotika, nevronalne mreže, kontrolerji, oddaljeni preizkusi, LabVIEW, MATLAB, Simulink
lega, hitrost...
Lab PC (strežnik) ^
primerjalna lega
uporabnik
daljinsko vodenje
- "L_______
Sl. 1. Shematski prikaz razvitega preizkusa za daljinsko vodenje SCARA robota preko medmre'ne povezave
Sl. 2. Na DSP robotski karti se izvaja razvita programska oprema
*Naslov odgovornega avtorja: Univerza v Mariboru, Fakulteta za elektrotehniko, računalništvo in informatiko, SI 18	Smetanova 17, 2000 Maribor, jure.cas@uni-mb.si
Prejeto: 12.1.2007 Sprejeto: 19.12.2007
Uporaba umetnih nevronskih mrež pri dvojnem/trojnem skupinskem sistemu krmiljenja dvigal
C. Erdem Imrak Tehnična univerza Istanbul, Turčija
Z umetnimi nevronskimi mrežami lahko, v primerjavi z običajnimi krmiljnimi sistemi dvigal, dobimo boljšo rešitev problema porazdelitve klicev potnikov. Zato v prispevku obravnavamo uporabo nevronskih mrež pri skupinskem sistemu krmiljenja. Predstavljena je smiselnost uvajanja umetnih nevronskih mrež. Skupinski sistem krmiljenja dvigal z nevronskimi mrežami lahko predvidi naslednje nadstropje postanka, z upoštevanjem priučenega z obdelavo vzorca sprememb v zahtevah potnikov. V prispevku obravnavamo uporabo umetnih nevronskih mrež za porazdelitev najprimernejših kabin v nadstropja z upoštevanjem zahtev potnikov. Umetne nevronske mreže smo uporabili pri dvojnem/trojnem skupinskem sistemu krmiljenja za zmanjšanje čakalne dobe potnikov. Za učenje nevronskih mrež smo uporabili povratni algoritem. Predstavili smo analizo prometa z dvigali in simulacijo rezultatov ter ju primerjali z običajnim sistemom krmiljenja dvigal.
© 2008 Strojniški vestnik. Vse pravice pridržane.
Ključne besede: dvigala, kontrolni sistemi, nevronalne mreže, povratno učenje
Sl. 2. Značilen model nevrona
*Naslov odgovornega avtorja: Tehnična univerza Istanbul, Fakulteta za strojništvo, Gümü§suyu, TR-34437, Istanbul, Turčija, imrak@itu.edu.tr	SI 19
Prejeto: 19.3.2007 Sprejeto: 28.9.2007
Uporaba analize uklona za napoved valovitosti pri stopnjevanem obdelovanju pločevine
Samir Lemeš* - Nermina Zaimović-Uzunović Univerza v Zenici, Fakulteta za strojništvo, Bosna in Hercegovina
Predstavljena je zamisel uporabe analize uklona za napoved valovitosti izdelka pri stopnjevanem oblikovanju pločevine. Pričnemo s predpostavko, da lahko z analizo uklona dobimo geometrijsko obliko pločevinastega izdelka po prvem opravilu stopnjevane deformacije. Ta oblika in vrednost kritične sile uklona se izračunavata z metodo končnih elementov v elastičnem področju. Analiza je opravljena za različne vrednosti premera orodja. Rezultati so preverjeni z meritvami geometrijske oblike izdelka po prvem opravilu stopnjevane deformacije orodja. Analiza je pokazala, da z manjšim gornjim orodjem dobimo enako kakovost izelka, izraženo s številom valov. © 2008 Strojniški vestnik. Vse pravice pridržane.
Ključne besede: preoblikovanje pločevine, stopnjevano preoblikovanje, uklon, metode končnih elementov
Sl. 1. Stopnjevano preoblikovanje majhnih	Sl. 2. Stiskalnica za stopnjevano preoblikovanje
izdelkov	P2MF 200x4 - Sertom, Milan, Italy
*Naslov odgovornega avtorja: Univerza v Zenici, Fakulteta za strojništvo, Fakultetska 1, 72000 Zenica, Bosna in SI 20	Hercegovina, slemes@mf.unze.ba
Prejeto: 13.3.2007 Sprejeto: 19.12.2007
Model zanesljivosti rezalnih orodij : Primer v proizvodnem inženirstvu in agrotehniki
• •
• •
P. Dašić1* - A. Natsis2 - G. Petropoulos3 'Visoka tehnična šola Kruševac, Srbija 2Agronomska univerza v Atenah, Oddelek za agrotehniko, Grčija 3Univerza Thessaly, Oddelek za strojništvo in industrijski inženiring, Grčija
Zanesljivost obdelovalnih sistemov je v veliki meri odvisna od zanesljivosti značilk rezalnih orodij kot najšibkejših členov zaradi njihove obrabe. Sistemi za obdelavo zemlje imajo precej skupnega z sistemi za obdelavo kovin - pomemben problem, ki se nanaša na uporabo obdelovalne opreme je obraba lemeža, ki značilno vpliva na kakovost obdelave zemlje in gospodarnost agrotehnične proizvodnosti.
V prispevku sta za določanje in modeliranje zanesljivosti R(T) podana primera obdelave kovin in obdelave zemlje: struženje legiranega jekla 20CrMo5 z orodji iz kubičnega borovega nitrida (CBN) in mešane rezalne keramike ter obdelava zemlje z uporabo preprosto nazobčane brane. V primeru obdelave kovin, skozi primerjalno analizo teoretičnih modelov porazdelitve, ki se dobro prilegajo eksperimentalnim podatkom, je tehtno izbran Gausov model za ponazoritev funkcije zanesljivosti obeh rezalnih materialov. V primeru obdelave zemlje je, z izvedenim postopkom optimizacije, ki temelji na zanesljivosti R(T) sistema za obdelavo, dobljena natančna oblika zoba orodja, ki omogoča najdaljši čas delovanja zoba zaradi zanesljive značilke.
© 2008 Strojniški vestnik. Vse pravice pridržane.
Ključne besede: obdelovalni sistemi, rezilna orodja, zanesljivost sistemov, poljedelska tehnika obraba orodij
Toplotne razpoke pri prekinjevanem Faza	rezanju
Hitrorezno jeklo
Karbidna trdina
Rezalna keramika
1.	Obraba proste ploskve
2.	Kotanjasta obraba
3.	Primarna brazda (zaradi dolbljenja
4.	Sekundarna brazda (oksidacijska obraba)
5.	Zunanja zareza odrezka
6.	Notranja zareza odrezka
zunanjega premera ali obrabne zareze)
Sl. 2. Mehanizem obrabe orodja za različne rezalne materiale [26]
*Naslov odgovornega avtorja: Visoka tehnična šola, Kosančićeva 36, 37000 Kruševac, Srbija, dasicp@ptt.yu
Prejeto: 20.2.2007 Sprejeto: 28.9.2007
Simulacijske metode pri načrtovanju postopka gradnje ladij
Boris Ljubenkov1* - Goran Đukić1 - Marinko Kuzmanić2 'Univerza v Zagrebu, Fakulteta za strojništvo in ladjedelništvo, Hrvaška 2ProMar Design, Kaštel Novi, Hrvaška
Namen tega prispevka je prikaz uporabe simulacijskih metod pri načrtovanju postopka gradnje ladij. Opisan je osnovno načelo postopka načrtovanja, glavni elementi spirale načrtovanja, tokovni diagrami in določeni medsebojni odnosi med posameznimi dejavnostmi. Tokovni diagram označuje točko, pri kateri rezultati simulacije vplivajo na načrtovani tok.
Za pripravo modela postopka gradnje ladij smo uporabili simulacijski programski paket. Prikazali smo primer možnosti in metod prikaza programskega paketa. V sklepu smo predstavili prednosti uporabe simulacijskih metod v postopku izgradnje ladij. © 2008 Strojniški vestnik. Vse pravice pridržane.
Ključne besede: ladjedelništvo, proizvodni postopki, simulacijske metode
Numorex	Plasma
i \
Corta
polavtomatski obdelovalni stroj
mostni žerjav
Sl. 2. Simulacijski model delavnice za gradnjo ladij
*Naslov odgovornega avtorja: Univerza v Zagrebu, Fakulteta za strojništvo in ladjedelništvo, Ivana Lučića 5, SI 22	10000 Zagreb, Hrvaška, boris.ljubenkov@fsb.hr
Prejeto: 1.8.2006 Sprejeto: 28.9.2007
Raziskave faktorja koncentracije napetosti ortotropnih pravokotnih plošč s končno širino ter krožno odprtino z uporabo trirazsežnega modela končnih elementov
Kambiz Bakhshandeh* - Iraj Rajabi - Farhad Rahimi Tehnična univerza Malek-Ashtar, Iran
V prispevku je predstavljena raziskava teoretičnega faktorja koncentracije napetosti pravokotnih plošč s končno širino ter središčno okroglo odprtino z uporabo trirazsezne analize končnih elementov. V prvem delu smo, s primerjavo analitičnih enačb za neskončne ortotropne in izotropne plošče, preverili model končnih elementov. S preverjenim modelom smo izračunali faktor koncentracije napetosti za otrhotropic pravokotne plošče s središčno okroglo odprtino v odvisnosti od nespremenljive natezne napetosti. Raziskali smo vpliv razmerja polmera odprtine in širine plošče za različne ortotropne plošče z različnimi ortotropnimi razmerji. Raziskava je pokazala, da je točnost analitične metode odvisna od ortotropnega razmerja in razmerja a/W. Pri velikem razmerju polmera odprtine in širine plošče, ima ortotropno razmerje le majhnen vpliv na faktor koncentracije napetosti. © 2008 Strojniški vestnik. Vse pravice pridržane.
Ključne besede: faktor koncentracije napetosti, metode končnih elementov, ortotropne plošče
2L
							
-	2 W	la	У				_„
					\	X	
				\J			
				2b			—
Sl. 1. Plošča končne širine s središčno eliptično odprtino
Sl. 2. Neskončna plošča s središčno krozno odprtino
*Naslov odgovornega avtorja: Tehnična univerza Malek-Ashtar, Letalsko-pomorski raziskovalni center, Shiraz, 7194915685, Iran, kambiz@bakhshandeh.ir	SI 23
Prejeto: 12.3.2007 Sprejeto: 27.6.2007
Razvoj kačastega robota
Jure Bezgovšek* - Igor Grabec - Peter Mužič - Edvard Govekar University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
Opisan je razvoj robota, ki posnema gibanje kač. Sestavili smo prototip robota iz šestih členov. Med posameznimi členi so generatorji navora, ki povzročajo spremembo oblike robota. Pri ustreznem spreminjanju oblike robota se zaradi anizotropnega zakona trenja med členi in podlago ustvarja pogonska sila, ki potiska robota naprej. Valovito spreminjanje oblike telesa določa krmilna spremenljivka, ki spreminja položajne kote generatorjev navora. Ustrezni krmilni signal potuje do robota preko brezžične radijske povezave. S prototipnim robotom smo opravili meritve povprečne hitrosti premikanja in tirnice robota pri različnih krmilnih parametrih. Podane so osnovne ugotovitve o gibanju kačastega robota. © 2008 Strojniški vestnik. Vse pravice pridržane.
Ključne besede: kačasti roboti, vodenje, serpentinasto vijuganje, silea trenja
Sl. 3. Kačasti robot sestavljen iz šestih med seboj vrtljivo povezanih členov
Sl. 4. Končna izvedba kačastega robota
*Naslov odgovornega avtorja: Univerza v Ljubljani, Fakulteta za strojništvo, Aškerčeva 6, 1000 Ljubljana, SI 24	jure.bezgovsek@fs.uni-lj.si
Poročila
Eksergija: Staro je novo
Na letnem srečanju združenja ASHRAE (American Society of Heating, Refrigerating and Air-Conditioning Engineers) v Long Beachu v Kaliforniji so člani Prve delovne skupine -Eksergijska analiza trajnostnih zgradb knjižnici združenja podarili izvirni zgodovinski članek, v katerem je bil prvič na svetu predstavljen termin "eksergija". Naslov članka, ki ga je napisal dr. Zoran Rant, profesor termodinamike na Fakulteti za strojništvo Univerze v Ljubljani, je Vrednost in obračunavanje energije. Članek je bil objavljen v slovenščini v Strojniškem vestniku - Journal of Mechanical Engineering, leto 1 (1955), št. 1.
Izvod revije, skupaj s prevodom izvirnega članka v angleščino je združenju podarilo slovensko društvo SITHOK (Slovensko društvo inženirjev za tehnologije hlajenja, ogrevanja in klimatizacije) kot priznanj e za prizadevanj e ASHRAE pri uveljavljanju trajnosti z eksergijo. Prav tako je bilo darilo priznanje članom Prve delovne skupine za prispevek k svetovni literaturi s področja eksergije in trajnostnega razvoja.
Terry Townsed, dr. Peter Novak in Jeff Littleton pri predaji prve številke Strojniškega vestnika
Dogodek je organiziral dr. Peter Novak, upokojeni profesor Univerze v Ljubljani in član Prve delovne skupine. Darilo sta sprejela član predsedstva Terry Townsed in izvršni podpresednik ASHRAE Jeff Littleton.
prevod iz ASHRAE Insights, leto 22, št. 10
Osebne vesti
Doktorat, magisteriji in diplome
DOKTORAT
Na Fakulteti za strojništvo Univerze v Ljubljani je z uspehom zagovarjal svojo doktorsko disertacijo:
dne 7. januarja 2008: Alira Srdoc, z naslovom: "Model vodenja kakovosti na osnovi baze podatkov in znanja" (mentorja: prof. dr. Alojz Sluga in akad. prof. dr. Ivan Bratko).
S tem je navedena kandidatka dosegla akademsko stopnjo doktorja znanosti.
MAGISTERIJI
Na Fakulteti za strojništvo Univerze v Ljubljani sta z uspehom zagovarjala svoji magistrski deli:
dne 22. januarja 2008: Stanislav Meglic, z naslovom: "Načrtovanje in vodenje dobavnih verig podjetja" (mentor: prof. dr. Marko Starbek) in Aleš Silic, z naslovom: "Integracija potisnega in vlečnega sistema proizvodnje" (mentor: prof. dr. Marko Starbek).
Na Fakulteti za strojništvo Univerze v Mariboru so z uspehom zagovarjali svoja magistrska dela:
dne 7. januarja 2008: Saša Ajd, z naslovom: "Modeliranje onesnaževanja podtalnice" (mentorja: prof. dr. Leopold Škerget in prof. dr. Matjaž Hriberšek) in Borut Križe, z naslovom: "Raziskava neposrednega legiranja in rafinacije pri proizvodnji nodulatorjev in cepiv ter vpliv spremenjene
tehnologije na okolje" (mentor: prof. dr. Alojz Križman);
dne 22. januarja 2008: Jasmin Kalju, z naslovom: "Ergonomski in estetski vidiki razvoja izdelkov" (mentor: prof. dr. Bojan Dolšak)
DIPLOMIRALI SO
Na Fakulteti za strojništvo Univerze v Ljubljani sta pridobila naziv univerzitetni diplomirani inženir strojništva:
dne 31. januarja 2008: Primož LAVRIN, Anton SVET.
Na Fakulteti za strojništvo Univerze v Mariboru sta pridobila naziv univerzitetni diplomirani inženir strojništva:
dne 31. januarja 2008: Tine JURŠE, Matej KUHAR.
*
Na Fakulteti za strojništvo Univerze v Ljubljani sta pridobila naziv diplomirani inženir strojništva:
dne 11. januarja 2008: Dejan ANTOLOVIČ, Janez VILFAN.
Na Fakulteti za strojništvo Univerze v Mariboru sta pridobila naziv diplomirani inženir strojništva:
dne 31. januarja 2008: Andrej KNAP, Jožef
VEK.