Modelling of ECM and EDM processes
Modelling of ECM and EDM processes
S.Hinduja(1)a,*,M.Kunieda(1)b
a School of Mechanical,Aerospace and Civil Engineering,The University of Manchester,UK
b Department of Precision Engineering,The University of Tokyo,Japan
1.Introduction
Because of their ability to machine tough,hard and heat-
resistant materials with complicated shapes,electro-chemical and
electrical discharge machining(ECM and EDM)processes were
?rst applied for the machining of aerospace alloys and die and
mould making in the1950s.Although both processes are
categorized as electrical machining to differentiate them from
conventional mechanical processes and both have similar machine
tool structures and applications,their principle and machining
characteristics are signi?cantly different.
The encouraging trend of ECM was not sustained for long
because of inherent dif?culties such as:
(i)controlling/predicting the process due to electro-chemical,
hydrodynamic and thermal factors;
(ii)predicting the equilibrium shape of the workpiece;and
(iii)determining the tool shape for a given workpiece geometry.
In the case of EDM,the uncertainty of tool electrode wear was
an early dif?culty.
Researchers have tried to overcome these dif?culties by
developing analytical and numerical models of the processes
but progress has been slow.The reasons for this are as follows.
(i)ECM and EDM are unlike conventional machining processes
such as turning or milling.In the latter,standard tool shapes are
usually used and some of the process parameters such as feed
rate and cutting velocity can be changed without affecting the
shape of the?nal workpiece.But in ECM and EDM,the tool
shape is unique to the workpiece geometry and a change in
feed rate results in a different workpiece shape.
(ii)Conventional processes usually require only one or two
physical phenomena to be modelled whereas ECM and EDM
require several,such as?uid?ow,gas evolution,chemical
reactions,heat generation at the electrodes and in the
electrolyte,and mass transport of the species.
In spite of the above complexities,in the last50years,
empirical,analytical and numerical models have been developed.
This paper reviews the development of these models and their
application in industry,and discusses some of the computing
issues.
2.Overview of modelling the ECM/EDM processes
Both ECM and EDM have several common modelling require-
ments.In both,the primary goal is to predict the shape of the
workpiece.EDM is predominantly a thermal process and therefore,
in its case,thermal modelling is more important.In the case of
ECM,however,determination of the current density distribution is
of primary interest.Once the distribution is known,other
parameters such as workpiece dissolution rate can be computed.
ECM has the advantage that there are no residual stresses but
they have to be evaluated in the case of EDM.
In the case of EDM,removal occurs at the discharge location
only.Even under the same pulse conditions,the thermo-
hydrodynamic and electromagnetic behaviours of the anode,
cathode,and working?uid materials bring about signi?cantly
different results in material removal.The material removal in
consecutive discharges is a cumulative result of single pulse
discharges.However,results of multiple discharges cannot be
obtained from a linear superposition of the results of a single
discharge,because the medium in the gap is composed of the
dielectric liquid,gas bubbles generated due to discharge and solid
debris particles.The composition,pressure,gap width,and
temperature vary both temporally and spatially,which makes
EDM simulation signi?cantly dif?cult.
CIRP Annals-Manufacturing Technology62(2013)775–797
A R T I C L E I N F O
Keywords:
ECM
EDM
Modelling
A B S T R A C T
The modelling of ECM and EDM processes requires not one but several models to simulate the different
phenomena that occur during machining.This paper reviews the models that have been developed to
simulate each of these phenomena,e.g.potential models to calculate the current density distribution in
ECM,thermal models for the plasma arc in EDM,moving boundary models to simulate the anodic
dissolution in ECM and probabilistic models to determine the discharge location in EDM.In addition to
discussing the relative merits of the techniques deployed in these models,the paper describes some
salient applications and concludes with desirable future enhancements to these models.
?2013CIRP.
*Corresponding author.
E-mail address:sri.hinduja@https://www.sodocs.net/doc/8f4222461.html,(S.Hinduja).
Contents lists available at SciVerse ScienceDirect
CIRP Annals-Manufacturing Technology
journal homepage:https://www.sodocs.net/doc/8f4222461.html,/cirp/default.asp
0007-8506/$–see front matter?2013CIRP.
https://www.sodocs.net/doc/8f4222461.html,/10.1016/j.cirp.2013.05.011
Workpiece shape change,in the case of ECM,is dependent primarily on the current density on the anode surface.The electrolyte’s electrical conductivity depends,amongst other factors,on the electrolyte temperature,which varies spatially, and the hydrogen bubbles released at the cathode surface.
In EDM,wear of the tool electrode necessitates the considera-tion of geometrical change of the tool electrode as well.To design the tool electrode shape,not only the gap width but also the tool electrode wear has to be considered.In the case of ECM,there is no tool wear except when reversed polarity is used to remove any solid debris adhering to the tool.
In continuous ECM,the tool is fed towards the anode at a constant value and in pulsed ECM,the cathode or anode is moved into a?xed position for each cycle.But in EDM,because the tool electrode is fed by a servo feed control,where the feed rate is based on the averaged gap voltage,modelling of the feed control is signi?cantly important to obtain the gap width distribution. Moreover,modelling of the feed control is also necessary to estimate the processing time,which normally results in a much greater error when compared with other cutting processes.
3.Modelling the ECM processes
Although the fundamental principles involved in ECM are well known[26,104],it is dif?cult to model the process because there are several physical and chemical phenomena,some of which occur simultaneously.The phenomena that need to be modelled are discussed below and summarized in Table1.
(i)Electro-chemical reactions.The chemical reactions occurring
at the anode and cathode cause ions,oxygen and hydrogen to be released from the electrodes and the electrolyte;this mass transfer as well as the spatial and temporal concentrations of the species has to be determined,leading eventually to the current distributions within the electrolyte and on the electrode surfaces.
(ii)Electrolyte?ow.The?ow of electrolyte through the gap may be laminar or turbulent.The chemical reactions at the anode and cathode result in hydrogen bubbles being released at the cathode and oxygen bubbles at the anode.The presence of
these bubbles causes the?ow of electrolyte to become two-phase;these bubbles affect the electrolyte’s conductivity,thus increasing the complexity of the process.
(iii)Thermal effects.The electro-chemical reactions cause heat to be generated in the double layer and in the bulk of the electrolyte(Joule heating).This heat energy causes the electrolyte temperature to increase,resulting in a further change to the electrolyte’s electrical conductivity.
(iv)Anodic dissolution.The electrochemical reactions occurring cause dissolution of the workpiece,resulting in its shape changing with time.This temporal change to the workpiece is modelled as a slow moving boundary problem.
The modelling of these phenomena requires the development of various models(see Fig.1)which are inter-dependent.
3.1.Multi-ion and potential models
The aim of these models is to determine the distribution of the current density on the electrode surfaces and in the electrolyte. With the multi-ion model,the total current i is given by the net?ux of charged species[117]:
i?
X
k
z k FˉN k(1)
where F is Faraday’s constant,z k the charge,andˉN k the?ux density of species k is given by:
ˉN
k
?àz k u k Fc k r VàD k r c ktc kˉv(2) where c is the molar concentration,D the diffusion co-ef?cient,u the mechanical mobility,V the potential,andˉv the velocity of the electrolyte.Two other fundamental equations are the electro-neutrality condition
X
k
z k c k?0(3) and the conservation of charge.
@c k
@t?àrá
ˉN
k
tR k(4) R k is the production rate of a species in the bulk of the electrolyte.Since in ECM the reactions occur on the electrode surfaces,R k is zero.
If the electrolyte velocity distribution is determined?rst,then Eqs.(1)–(4)can be solved to give(i)the voltage distribution in the electrolyte,and(ii)the concentrations of each species(c k)and hence the total current.Researchers have developed models which take all three components of Eq.(2)into account,i.e.migration, diffusion and convection.Qui and Power[137]developed a two-dimensional boundary element model which they applied to a parallel cell reactor and predicted the change in the cathode shape. Bortels et al.developed a more comprehensive model but it was based on two-dimensional?nite elements and they referred to their model as a multi-ion transport and reaction(MITReM)model
Table1
Modelling requirements for ECM and EDM.
ECM EDM
Aim of modelling Predict workpiece shape,
cycle times,optimum
process parameters
Predict workpiece shape,
tool wear,gap width,
wire vibration,
temperature distribution
Physics of the process Electro-chemical
modelling(potential and
current density
distributions,ion mobility,
mass transport)thermal
modelling,?uid?ow
modelling
Thermal modelling
(heat-affected zone,
residual stresses)?uid
?ow electromagnetic
modelling
Geometry considerations Modelling workpiece
shape(ef?ciency,
Faraday’s law)
Modelling both tool and
workpiece shapes
Power source and polarity Continuous and pulsed
dc;Tool(cathode)
Workpiece(anode)
Pulsed dc and ac
Machine control Constant feed rate
(continuous ECM)or
static(pulsed ECM)
Adaptively controlled
Adaptive control model None Pulse conditions,feed,
jumping
Inverse problem Predict tool shape Predict tool shape
Information required in database Electrical and thermal
properties of electrolyte,
valency,over-potential
and current ef?ciency
Discharge delay time,
energy distribution,
plasma diameter,
removal per pulse,
electrical and
thermophysical properties
of electrodes
Fig.1.Modelling the ECM process(adapted from[125]).
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797
776
[9].These models are applicable to both electro-chemical machining and deposition.
Generally,in ECM the following assumptions can be made. (i)The time-dependent concentration term is zero.
(ii)The electrolyte is continuously refreshed and therefore the concentration gradient terms can be ignored.
Because of these assumptions,it can be shown that the above equations reduce to Laplace’s equation which,in the context of ECM,is often referred to as the potential problem.
k e ráer VT?0(5) where the electrical conductivity k e is given by
F2
X
k z2
k
u k c k(6)
One can use the multi-ion mass transport model for ECM to calculate the concentration of the species and the current density within the electrolyte;this gives a greater insight into the process [30,31]but at the expense of greatly-increased processing time.On the other hand,the potential model is popular with researchers because of its simplicity and the relative speed with which the results can be computed.
3.1.1.Empirical models
Over the years,researchers have used various techniques to predict the workpiece(anode)shape.The advent of computers in the1970s enabled researchers to use numerical techniques to solve for the potential distribution in the inter-electrode gap but prior to this,researchers developed empirical and mathematical models.
The earliest empirical equations to predict the overcut and side gaps in drilling were developed by Koenig and Pahl[75].Koenig and Degenhardt[73]enhanced these equations for completely insulated and partially bare cylindrical tools.They also developed nomograms from which one could determine the side gap for a given tool tip radius,land width and equilibrium gap.Ippolito and Fasalio[59],using a multiple regression technique,also developed equations to predict the over-cut.The main disadvantage with these empirical equations is that they are valid only for the range over which the experiments were conducted.
3.1.2.Mathematical models
Of all the analytical techniques,the simplest and effective approach is the cos u method?rst pioneered by Tipton[159].For workpiece shape prediction,u is the angle between the normal to the tool surface and the feed direction.Kubeth[85]and Kawafune et al.[68]also suggested the same method but used the complementary angle and referred to it as the sin u method.In the cos u method,the tool is subdivided into several straight-line segments and the corresponding workpiece segment is deter-mined to be at a distance of h e=cos u where h e is the frontal equilibrium https://www.sodocs.net/doc/8f4222461.html,wrence suggested something similar but referred to it as the?at inclined cathode theory[97,98].However, according to Jain et al.[61],the cos u method gives rise to dif?culties when u>458,and the tool has sharp corners or has radiused edges because it assumes parallel?ux lines and neglects the effect of stray currents.If the tool has a complex pro?le,then this method is not to be recommended for obtaining the?nal equilibrium workpiece shape;however,most researchers,when using an iterative numerical technique to determine the workpiece shape,use the cos u method to determine the starting shape.
A mathematical technique,which can predict the workpiece shape albeit under idealized machining conditions,is that due to Collett et al.[17].They used a conformal mapping technique and were able to analytically determine the workpiece shape for plane-faced tools with and without insulation on the vertical face. Hewson-Browne[50]extended the work of Collett et al.by considering a tool with rounded corners and partial insulation.This method can only be applied to a limited number of geometries that can be transformed to the complex plane;moreover it is restricted to two-dimensional problems and transient workpiece shapes cannot be predicted.However,it can provide exact workpiece shapes for certain tool geometries and these shapes can be,and have been,used by researchers to benchmark their solutions.
A more recent analytical-cum-computational method has been suggested by Hocheng et al.[51]to predict the shape of an EC drilled hole.At each time step,they consider a point on the workpiece surface as being in?uenced by all the point sources on the tool and the amount of dissolution at an anodic point is obtained by integrating over the entire width of the tool.
3.1.3.Numerical models
As processing power began to become readily available in the late1970s and early1980s,researchers started to use numerical methods such as the?nite difference and?nite element methods to determine the potential distribution within the inter-electrode gap.The boundary element method(BEM)remained in the shadow of the other two methods and it was not until the early1980s that it was used to model the ECM process.
The?nite difference method(FDM)was the?rst of the numerical methods to be used;Tipton[159]is generally credited for having pioneered the application of this method to ECM.Tipton tried different methods of moving the workpiece;instead of moving a point along a direction normal to the surface,he considered only the vertical component of the cut vector.With this method the?nal equilibrium shape is not affected but the intermediate transient shapes do not correspond to reality. Klingert et al.[71]determined the primary and secondary current distributions with the activation over-potential modelled as a linear/logarithmic current density-dependent https://www.sodocs.net/doc/8f4222461.html,wrence [98]predicted the shape he would obtain with a semi-cylindrical tool and veri?ed his predictions with both experimental and analytical results.In their FDM model,Koenig and Humbs[74] made the current ef?ciency dependent on the current density and temperature whereas Dabrowski and Kozak[19]combined their FDM model with the analogue method.Subsequently Kozak developed FDM-based models suitable not only for EC drilling but also for die sinking[79].This model considered the change in the electrolyte’s electrical conductivity due to temperature increase and void fraction.Like Klingert et al.[71],Prentice and Tobias[135] also developed their FDM model to include primary,secondary and tertiary current distributions but they applied their model to electroplating.
It is generally accepted that the standard FDM cannot deal accurately with curved boundaries;it assumes a linear variation of the governing variable.Both these factors make it necessary to deploy thousands of grid-points in the inter-electrode gap, although Nanayakkara[112]and Nanayakkara and Larsson [113]tried to reduce the number of grid-points by assuming a quadratic variation for each grid-point.These disadvantages can be overcome by using the?nite element method(FEM)because higher-order iso-parametric elements can be used,the sides of which can be curved.Therefore,no approximations have to be made when representing curved tool and workpiece shapes;also, these elements allow a quadratic or cubic variation of the potential, thus obviating the need for a large number of elements.
The FEM was?rst used in ECM by Pandey and Jain;in[63]they modelled the inter-electrode gap with a very coarse linear triangular mesh and in[62]they repeated the same but this time they obtained the temperature distribution within the gap.The geometry of the reactor cell was limited to two parallel plates. Alkire et al.[3]highlighted the dif?culties in obtaining accurate distributions of the potentials in the vicinity of a singularity. Assuming ideal machining conditions,Hardisty et al.predicted the workpiece shape generated by a parabolic tool and veri?ed it with an analytical solution[44].
Although the use of?nite elements is a considerable improve-ment over?nite differences,it has one serious drawback when
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797777
modelling the ECM process.Because the workpiece shape changes after every time step,one has to re-mesh the entire domain after every few time steps.In the1980s,when the availability of processing power was still limited,Brookes[11]found the re-meshing of the inter-electrode gap very time-consuming.Instead, he used the same mesh but after every time step,he deleted some elements,modi?ed others and improved the shape of the elements close to the workpiece boundary.Unfortunately,he found that modifying the mesh was almost as time-consuming as re-meshing the inter-electrode gap.
Another disadvantage of the FEM is that the values of the voltage gradient,denoted by@V=@n or q,at nodes on the
workpiece surface are not directly calculated.They are subse-quently determined from the computed values of potentials at the nodes in an element and this gives rise to discontinuity of@V=@n between adjacent elements.
The BEM overcomes these problems to some extent.Although Christiansen and Rasmussen[15]were the?rst to apply the integral equation method to modelling ECM,it was Hansen who employed it more generally for different axi-symmetric con?g-urations[42].Narayanan et al.also developed a general2-D model but were able to demonstrate the accuracy of their model by comparing it with known analytic solutions[115].Since the BEM reduces the dimensionality of a problem by one,it means that a3-D model of the ECM process requires the surfaces of the tool and workpiece to be discretised with2-D triangular elements.Fig.2 shows a boundary element model consisting of a workpiece(grey surfaces)in which a slot is to be machined by a tool(golden surfaces)with a rectangular cross-section.The?gure shows the tool about to start machining.Virtual surfaces(green)are introduced to form a closed shell.
It is interesting to compare the accuracy that is obtainable with the different numerical methods.Narayanan et https://www.sodocs.net/doc/8f4222461.html,ing a plane-faced rectangular tool,with and without insulation,did one such comparison[115].They compared the results obtained from FD,FE and BE models with those obtained by Collet et al.[17]and Hewson-Browne[50]for a rectangular tool(Fig.3)under ideal machining conditions(Table2).Assuming that the results from the mathematical methods are exact,it is clear that the results obtained by PERA[132]and Lawrence[97]using?nite differences are considerably in error whereas Brookes[11]using?nite elements,predicted the ratios h c/h g and h o/h g to be0.83and 1.191,respectively,for the tool with no insulation(see Fig.3).The BEM,using quadratic elements,calculated these ratios as0.8043 and1.1586respectively.The FE value for h1/h g is considerably far away from the theoretical value of1.0as calculated by Hewson-Browne;the BE value of1.079is considerably closer,demonstrat-ing the superiority of the BEM.It should be noted that these results were obtained in the mid-1980s and the meshes deployed then were not as?ne as those one would use today.However,the example serves to illustrate the relative accuracy that can be obtained with the different methods.In summary,the accuracy depends upon the type of iso-parametric element used(higher-order elements give more accurate results)the time step(if a large time step is used,oscillations are induced in the workpiece shape which become dif?cult to suppress)and the mesh density.
The BEM has one major disadvantage.The electrolyte’s electrical conductivity(k e)varies in the inter-electrode gap because of Joule heating and the release of hydrogen bubbles. This variation can be modelled to some extent by subdividing the gap into several zones and assuming a constant value of k e for each zone.Subdivision into zones makes the use of the BEM less appealing.If one can assume that k e is constant in the entire gap, then the BEM is the ideal choice;otherwise,if a more accurate model of the process is required wherein the velocity,temperature and k e variations are taken into account,then the FEM should be used.
3.2.Workpiece shape change model
ECM is a slow-moving boundary problem and for modelling purposes,the total machining time is divided into several time steps.The rate of dissolution at a point on the workpiece is governed by Faraday’s law
dh
dt
?
h M
r zF J(7) where h is the current ef?ciency,M the atomic mass,h the inter-electrode gap,and J the normal current density which is given by J?k e q(8) Eq.(7)is usually solved using the marker method although more recently the level set method(LSM)has been deployed.With the former,using the explicit Euler integration scheme,the gap at node i on the workpiece surface,after the k th time step,is given by
h kt1
i
?h k itD tá
h M
r zF Jàf cos
u
(9)
where D t is the time increment and f the feed rate.Eq.(9)requires not only the value of J but also its direction.At a node lying on the workpiece boundary,one of two cases may arise.In the?rst case, the elements meeting at this node may share a common tangent. For example,in two-dimensional models,if the workpiece surface is modelled as a B-spline curve using higher-order elements,then slope continuity is ensured.Alternatively quadratic elements can be used and slope continuity enforced[115];in such cases,a single value of the voltage gradient q is computed by the BE or FE program,and its direction is also uniquely de?ned,i.e.normal to
Fig.2.A BE model for a milled component[129].
Fig.3.Equilibrium frontal and side gaps[115].Table2
Comparison of results for a rectangular tool[115].
Tool with no insulation Tool with insulation
h c=h g h o=h g h c=h g h o=h g h1=h g
Collett[17]0.80 1.1590.80 1.159
Hewson-Browne[50] 1.0
Pera[132] 1.0 1.7 1.0 1.7
BEM[115]0.8043 1.15860.74090.7548 1.079 FE[11]0.83 1.1910.7750.875 1.25 FDM[97]– 1.3––
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797 778
the tangent.In the second case,there is a geometric discontinuity (e.g.sharp corner)at the node and the voltage gradient q becomes multi-valued.One solution to this problem,which is popular with researchers,is the method suggested by Alacron et al.[2].According to this method,a single value of q is forced on the solution and its direction is determined by assuming a linear variation of the potential over the element edges meeting at the node.
In3-D BE models,linear triangular elements are used resulting in the workpiece surface becoming faceted.Therefore,there may be several elements meeting at a node,with each element having its own normaleˉn eT.In such cases,the unit normal(~n nd)at that node is given by the average of the element normal vectors[135],
~n nd ?
1
e
P e?1
N e
n e!
P e?1
N e
n e!
2
64
3
75(10)
where N e is the number of elements meeting at node i.In some cases,the above does not yield a valid normal direction.According to Purcar[136],for a normal to be valid,it should be visible from all the faces.He deployed a method?rst suggested by Kallinderis and Ward[67]wherein the two faces subtending the most acute angle at a node is determined,followed by the plane bisecting them.The normal is forced to lie on this plane and its direction is chosen so that it makes equal angles with the remaining faces meeting at node i.
Another approach to solve the sharp corner problem with BE models is by the use of discontinuous or partially discontinuous elements[131].With these elements the nodes are not located at the ends of an edge but along the edges or even within the element. Yet another technique to cope with the corner problem without changing the element type is to use the double node technique introduced by Brebbia[10].There are disadvantages with the use of double nodes.The elements meeting at the double node are displaced in different directions,creating a virtual gap between the nodes which have to be bridged somehow.Also,the use of double nodes increases the number of equations,and hence the computing time.However,in the case of ECM,double nodes are ideal to represent sharp corners on the tool surface[129](e.g.at the junction of the end and side faces)because there is no relative movement between nodes on the tool surface(at the end of every time step,all the nodes on the tool undergo the same rigid body movement).
Another modelling problem occurs when the cathode(or anode)is adjacent to an insulating surface and the angle between the two surfaces is greater than p/2(see Fig.4(a)).For example,during the n+1th time step the cathode/anode will be moved due to deposition/ dissolution resulting in the two surfaces becoming disconnected (Fig.4(a)).In such cases,Deconinck suggests either an additional virtual element to connect the end of the current electrode pro?le with the insulated surface(Fig.4(b))or extending the electrode pro?le until it meets the insulated surface(Fig.4(c))[28].
The accuracy of the computed workpiece shape depends also on the magnitude of the time step.For a simple cell geometry consisting of two parallel plates,Hardisty et al.predicted the number of time steps that would be required to reach equilibrium; they tested their improved algorithm for small and large initial starting gaps[43].Narayanan[114]showed that too large a time step induced oscillations in the computed workpiece surface that were dif?cult to suppress.Deconinck[27]showed that the error depends upon Wagner’s number W and the magnitude of the inter-electrode gap h.Based on a one-dimensional analysis,Purcar[136] suggested that the time step is given by
D t?
??????????e
r
á
J av
max
(11)
where J av and J max are the average and maximum values of the current density respectively,e the permissible error,m the number of time steps and b=1/(2(W+h min)3)where h min is the minimum gap.
When the explicit scheme is used,Volgin and Lubiynov found that,irrespective of which of the four variations(right,left,central and upwind)of the explicit difference method is used,numerical instability sets in when sharp corners are encountered[168].A sharp corner results in the normal being de?ned ambiguously and subsequent iterations cause a swallow tail or self-intersection in the workpiece pro?le to be formed.Special topological routines have to be developed to detect and eliminate them.Fig.5shows a workpiece containing a recess and as it grows,its shape self-intersects[136].
This problem can be avoided by using the level set method,which was?rst pioneered by Sethian[144]and later described in detail in Ref.[145].This method has been applied to other dynamic problems such as?ame propagation and waves.In this method,a scalar implicit function?is used to represent the moving front and also its evolution.For example,the following implicit function would be used to represent the evolution of the workpiece front in ECM @?
@tv r??0(12) where v is the velocity with which the dissolution of the workpiece occurs.The LSM was?rst applied to electro-chemical machining by Volgin and Lubiynov[168]who used the stationary formulation of the LSM.This method is more computing intensive than the marker method,although the computing time can be reduced to some extent if the velocities are calculated only for a band of nodes on either side of the workpiece surface.It has been applied in2-D models but has yet to be applied to3D ECM models.
3.3.Thermal models
Most of the early potential models assumed a constant value of conductivity when computing the current density distribution.In reality,the increased temperature of the electrolyte and the presence of hydrogen gas in the inter-electrode gap affect k e.There are two empirical equations for k e and the one suggested by Thorpe and Zerkle[158]is given by:
k e?k e0e1àaTm?1tgeTàT0T (13a)
where a is the void fraction,T the electrolyte’s current temperature,g the conductivity constant of the electrolyte and
Fig.4.Junction between electrode and insulated surface[27].
Fig.5.Self-intersection of the work-piece pro?le[136].
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797779
m a constant varying between1.5and2for the heterogeneous mixture of liquid and gas.Suf?x0refers to the initial properties of the electrolyte at the time of entering the gap.The second empirical equation is due to Riggs[139]and Riggs et al.[140]: k e?k18e0:024eTà18Te1àaT1:5(13b)
where k18is the conductivity at188C.Needless to say,use of the above equations requires knowledge of the temperature distribu-tion within the electrolyte which can be determined by solving the following convection-diffusion equation,
r C p @T
@tt
r c pˉv r T?ráek t r TTtP bulk(14)
where r is the density,c p the speci?c heat,ˉv the velocity,k t the thermal conductivity and P bulk the heat load.Researchers avoid using this equation in its entirety because of the following dif?culties:
(i)a priori knowledge of the velocity distribution of the electrolyte
within the gap is required;and
(ii)the thermal conductivity of the electrolyte in the equation is the sum of the molecular and turbulent conductivities;the latter is determined from the k–v Reynolds averaged turbulence model.
The heat load usually consists of P bulk,the Joule heat generated in the electrolyte,and P dl,the heat generated in the double layer. The latter can be neglected as most of it goes to heat the anode and cathode surfaces.
Most of the early researchers used a one-dimensional approach to calculate the increase in electrolyte temperature.In this approach,the temperature at a section along the gap is obtained by performing a simple energy https://www.sodocs.net/doc/8f4222461.html,ing this one-dimen-sional approach,Clark and McGeough[16]compared their predictions with experimentally measured temperature values. Although a good correlation was found between the two sets of temperatures at the exit of the gap,there was considerable discrepancy between the two sets at intermediate sections.The measured temperatures of the electrolyte were much higher.This is understandable because the predicted temperatures did not take into the account the additional resistivity caused by the hydrogen bubbles.Loutrel and Cook[100]observed through a microscope that the average void fraction of the?owing bubble layer is0.5; they,as in[16],used a one-dimensional model but took into account the increased resistivity due to the bubbles and were able to make more accurate predictions of the workpiece shape.
Instead of having to assume a value for the void fraction,many researchers have developed a two-phase numerical model for the electrolyte?ow which they solved using?nite differences.Some of these two-phase models incorporate one-dimensional?ow by averaging the variables across the width of the gap at each section [37,52].Using a two-dimensional?ow model,Hourng and Chang were able to compute the temperature increase along the?ow stream under equilibrium conditions taking into account the void and other process parameters[53].
Jain and Pandey were one of the?rst to use a two-dimensional FE model to predict the current densities and hence calculate the increase in temperature of the electrolyte[62].Their FE results consistently underestimated the temperature increase,probably because of the very coarse mesh used.
With the advent of pulsed electro-chemical machining(PECM) in the late1980s,the interest in electro-chemical machining has been revived because PECM offers,by virtue of smaller gap sizes (10–100m m),better accuracy and surface?nish.The determina-tion of the increase in electrolyte temperature depends largely on the Strouhal number(S h),which during the off-time is given by
S h?
L
vt o
where L is the length of the electrode and t o is the pulse-off time.
When S h(1,there is suf?cient time for the electrolyte to be
replenished and by-products to be?ushed away;hence the system
can be considered to have completely recovered during the rest
period.Under these circumstances,it is suf?cient to model the
PECM process by considering the effect of just one pulse(or group
of pulses)applied during the on-time.Kozak,Rajurkar and Ross
adopted this approach and they derived the characteristic PECM
equations using a simple energy balance and then solved them
using?nite differences[82].In a more re?ned mathematical model
for PECM,Kozak,Rajurkar and Wei suggest a critical upper limit for
the pulse-on time for different process parameter settings[83].
Fig.6shows the variation of the charge density with gap size,for
different pulse-on timeset pT.One can operate below this critical
time limit by ensuring that for a given current density,the
corresponding pulse time and gap size are to the left of the dotted
line.
Smets et al.[148]developed an analytical model that was based
on the complete convection-diffusion equation;they studied the
temperature evolution in PECM in a rectangular gap taking the
electrodes into consideration.They showed that the temperature
transient curve during the off-period comprises two parts(see
Fig.7);in the?rst part which shows a steep fall in temperature,
most of the heat is transported away by the electrolyte due to
convection.In the second part,the fall in temperature is more
gradual as heat?ows from the electrode across the thermal
boundary layer and into the electrolyte.The second part has
therefore a much higher time constant than the?rst.
Even when the Strouhal number is small,the temperature
evolution is transient in nature taking several cycles to reach
steady state.Since the time scale for the pulse-on time is far
smaller than the time required for the temperature to reach steady
state,a complete transient may become computationally very
expensive as several thousand time steps will be required.Smets
et al.overcame these dif?culties by developing several analytical
and numerical models to calculate the averaged and pulsed
temperature history[149,150].Their quasi-steady state short-cut
(QSSSC)model is capable of calculating the average temperature
until a certain time period and then switching over to smaller time
steps to determine the cyclic temperature variation over the next
few cycles.
Fig.6.Variation of current density with gap size in PECM[83].
Fig.7.Time-temperature evolution in PECM for a point located on the anode and at
the gap exit(S h=0.006,pulse period=0.1s)[148].
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797
780
3.4.Over-potential
The over-potential of a cell has three components,i.e. activation,concentration and resistance over-potentials.Whilst these over-potentials are dependent,amongst other factors,on the current density and?ow,it is dif?cult to characterize each of them https://www.sodocs.net/doc/8f4222461.html,ually,for the sake of simplicity,the three over-potentials are combined together and represented by a single over-potential(V pol)and only its relationship with the current density is taken into consideration either by a Tafel-like logarithmic relationship
V pol?atb ln J(15a) or by a simple linear relationship
V pol?atbJ(15b) where a and b are constants or by using the Butler–Volmer equation[117]
J?J o?e a a FV pol=RTàe a c FV pol=RT (15c) where R is the gas constant,a a and a c are the transfer co-ef?cients, and J o the exchange current density.The values of J o,a a and a c have to be experimentally determined.
Assuming a linear relationship and a simpli?ed analysis,Altena showed that the current density is given by[4]:
J?eVàaTk e
k e(16)
where V is the applied voltage.The product b k e expresses the effect of concentration on the current density and becomes important when the gap sizes are small as is the case in pulsed ECM.In the case of continuous ECM where the gap sizes can be as much as0.5mm,h is much greater than b k e.Hence,serious errors are not introduced if V pol is assumed to be independent of J. However,when the inter-electrode gap is in the range of100m m or less,the effect of the concentration over-potential becomes more pronounced and neglecting the term bák e will introduce errors.
Consideration of over-potential using the logarithmic or linear relationship presents an awkward problem because both V pol and J are not known.It causes the problem to become non-linear,thus requiring an iterative solution at every time-step.Danson et al.[20] and Adey[1]used the following iterative method.
V kt1 pol ?V kà1
pol
tCeV k polàV kà1
pol
T(17)
where C is a damping factor and k the current iteration.Prentice [134,135]also used a similar iterative technique but made the damping factor a variable dependent on the normalized changes in V pol and J,and also on Wagner’s number.When the iterative technique of Danson et al.was tried for a stepped tool with linear and logarithmic over-potential relationships,it took more than twelve iterations per time step to converge.A much faster method based on Newton-Raphson’s method was developed by Narayanan and who was able to achieve convergence within three iterations [114].
One of the dif?culties that modellers face in taking polarization into account is the lack of experimental data for the constants in Eq.(15).Altena[4]experimentally determined the polarization voltage data for different amounts of concentration of NaNO3 electrolyte;one set of his results is shown in Fig.8which clearly exhibits a linear relationship between V pol and J.
To investigate the effect of over-potential,Narayanan assumed the following normalized linear over-potential relationships for the tool and workpiece surfaces and an equivalent logarithmic relationship[114].
V polàcathode?0:2J(19a) V polàanode?1:0à0:2J(19b)
The computed workpiece shapes are shown in Fig.9,from which it is clear that polarization decreases the equilibrium gap. Also the computed workpiece shapes with linear and logarithmic over-potential relationships are different only in the vertical section of the workpiece,a section where the current density values are relatively small.
3.5.Current ef?ciency
One of the parameters required for modelling is the current ef?ciency(h);its value depends upon whether the electrolyte is passivating(e.g.NaNO3)or not(e.g.NaCl).If it is a non-passivating electrolyte,then h can be assumed to be a constant.Otherwise,it is a function of the current density,pulse time and electrolyte concentration.It is in?uenced to a lesser extent by the build up of the anions on the anode surface[23].
Instead of h,Kozak et al.advocate the use of the electro-chemical machinability coef?cient(k v)which is given by
k v?
h:k c
r(20)
where k c is the electrochemical equivalent of the workpiece material[81].Kozak et al.argue that k v should be determined experimentally because,during alloy dissolution,the electroche-mical equivalent of each constituent is different from that when the constituent materials are dissolved individually.Instead of having to determine values for two parameters(i.e.h and k c) experimentally,they showed that only one is necessary,i.e.k v which is given by f=J o where J o is the mean current density;they then went on to determine experimentally the value of k v for passivating and activating electrolytes when machining different alloys.For example,using a13%water solution of NaNO3and machining an alloy NC10,they obtained the relationship k v=1.64–2.13eà0.034J which,of course,can be easily programmed into a system.
Altena investigated the effect of current density,pulse time and concentration on the current ef?ciency.Fig.10shows one set of
Fig.8.Variation of polarization voltage with current[4].
Fig.9.Workpiece shapes for different over-potentials[114].
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797781
results that he obtained using a water solution containing250g/l of NaNO3[4].These results were captured by Altena into a single hyperbolic tangent equation:
h?aetanhebt ptcTTáJtedát pteTtf(21a) where a–f are constants and t p is the pulse on-time.In the case of continuous ECM,the relationship is much simpler,i.e.
h?aáetanhebáJtcTT(21b) Instead of having to experimentally determine the value of h, Van Damme et al.,in a seminal work,developed a FE model to predict the current ef?ciency when machining steel using a solution of NaNO3as the electrolyte[165,166].They predicted the current ef?ciency by calculating three variables at each node,i.e. the potential,and ion concentrations c Me z+and c OHà.These variables were obtained by solving three differential equations,the ?rst of which was the potential equation.Eq.(2)(but without the migration term)was solved twice,once for each of the two ionic species.The other highlight of this work is the introduction of a water depletion function which was used as a weighting factor for the assumed current density function.Their results for different pulse-on times and for a particular concentration are shown in Fig.11along with the experimental results obtained by Altena[4]. Although there is some discrepancy between the two sets of results,the proposed procedure is very encouraging and the way forward for modelling current ef?ciency.
3.6.Flow models
Flow is an important process parameter because in addition to forming an electrical bridge between the two electrodes,it transports away most of the heat energy(and solid debris) generated during machining.How effectively it does this depends upon the?ow patterns generated,whether there are regions within the gap which are starved of the electrolyte and/or regions where the electrolyte forms eddies.
Single phase?ow is obtained by solving the incompressible Navier–Stokes equations.
rˉv?0(22a) r
@ˉv
@ttˉvárˉv
?àr ptmDˉv(22b)
where p is the pressure,ˉv is the velocity and m is the dynamic viscosity.
Jain et al.were among the early researchers to develop a numerical model which accounted for the void fraction,and temperature and pressure changes in the inter-electrode gap[64]. However,the model was very speci?c in its application as it required the anode and cathode to be cylindrical which made it possible to solve some of the equations analytically.Hourng and Chang[52]initially generalized the problem by developing an integrated model in which the bubbly two-phase?ow was?rst computed using Navier Stokes equations using a one-dimensional model.The computed velocities were then fed into a thermal energy model to calculate the temperature distribution in the electrolyte.The temperatures were used to update the values of the electrolyte’s electrical conductivity,followed by a re-computa-tion of current distribution in the gap.This iterative process was repeated until convergence of the anode shape was obtained. Hourng and Chang were able to predict the spatial variation of the process variables within the inter-electrode gap.
However,one-dimensional?ow models have limited use especially when the anode and cathode shapes contain sharp bends because they cannot identify any eddies or separation of the electrolyte.In a subsequent paper,Hourng and Chang[53] enhanced their model by considering the?ow to be two-dimensional.They were able to demonstrate that better workpiece accuracy was achievable with this model than that with a similar 1D model.A very convincing example is shown in Fig.12;it clearly demonstrates that?ow,even if it is single phase,should be modelled at least in two dimensions[29].This?gure clearly shows that re-circulations are formed,on the upstream and downstream sides,near the two concave vertices of the cathode.These re-circulations reduce the ef?ciency of heat removal and actually result in an unsymmetrical workpiece shape(see Section5.4).This example also shows that a visualization of the?ow in the gap could be of invaluable assistance in the design of the cathode shape.
But?ow in ECM should be modelled as two-phase because suf?cient quantities of hydrogen are released at the cathode. Oxygen is also released at the anode but the quantity is not signi?cant.Since a considerable portion of the inter-electrode gap, especially the downstream part,is occupied by hydrogen gas, considering?ow as single phase would be an approximation.
Because modelling two-phase?ow in two-dimensions is not straightforward,early researchers simpli?ed the problem by considering the two-phase?ow only as a one-dimensional model. This made it possible for researchers,notably Thorpe and Zerkel [158]to study the dynamics of the system of the process analytically.Loutrel and Cook[100]were among the?rst to build a two-phase1D numerical model of the process;in this model they assumed a linear increase of the void fraction from gap entrance to exit.This causes the gap to decrease towards the exit.The use of a
Fig.10.Current ef?ciency for different current densities and pulse-on times for concentration of250g NaNO3/l.[4].
Fig.11.Predicted current ef?ciency for different current densities and pulse-on times[166].
Fig.12.Electrolyte?ow in the inter-electrode gap[29].
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797 782
two-phase two-dimensional model enabled Chang and Hourng [14]not to make this assumption.Instead the rate of dissolution of hydrogen gas was incorporated into the model using Faraday’s laws of electrolysis and by so doing they were able to obtain a much better correlation with experimental results.3.7.Tool shape prediction model
Predicting the tool (cathode)shape for a given workpiece shape is often referred to as the inverse problem.The inverse problem is more important than the direct problem because,in practice,several trial tools have to be manufactured resulting in high development costs.Yet tool design has received far less attention from researchers probably because the problem is ill-posed,i.e.in some cases,for a given workpiece shape,there is no corresponding tool shape and sometimes for a given workpiece shape,there are several solutions.The inverse problem is also a free boundary problem because one of the boundaries,i.e.tool boundary,is unde?ned.There are three main approaches to solve this inverse problem,i.e.analytical,transformation to the complex plane and embedded.The analytical method is due to Tipton who suggested that the tool surface distance is proportional to cos u ,where u ,in this particular case,is the angle between the normal to the anode and the tool feed direction [160].Since this method gives rise to dif?culties,it is now normally used by researchers as a starting approximation to the required cathode shape.
As early as 1968,Krylov [84],followed by Nilson and Tsuei [118,119],showed that it is possible to calculate the cathode shape directly if:
(i)the workpiece shape is de?ned using an analytic function,e.g.Fourier series;and
(ii)the shape can be transformed from the physical plane to the
https://www.sodocs.net/doc/8f4222461.html,cey [94]demonstrated this transformation technique by calculating the cathode surface for a concave hyperbolic anode surface and for different values of k e eV =f T.
However,there are dif?culties with this approach.Only the simplest shapes can be transformed to the complex plane;also a corner on the anode surface gives rise to a singularity,the consequence of which is that the direct solution can be obtained only for small values of k e eV =f Ti.e.high feed rates,low anode voltage or low electrolyte conductivity.
Hunt transformed the potential equation to a second-order differential equation which he integrated and then,using vertical lines or a multi-grid,searched for points at which Laplace’s equation was satis?ed [54].The curve joining these points de?ned the cathode shape.Hunt also suggested the embedded technique [55]which removes the ill-posedness associated with the inverse problem.He did this by searching through several direct solutions obtained by successive modi?cations to the cathode shape until a shape was found which matched,or was close to,the required anode shape.This embedded approach,in some form or another,is still the most popular technique deployed by researchers.
Narayanan et https://www.sodocs.net/doc/8f4222461.html,ing the BEM utilized the fact that the conditions on the anode surface are over-speci?ed [116].Given a workpiece shape and the feed rate,the required value of the
voltage gradient q r is equal to the dissolution rate,i.e.f cos u =˙M
where ˙M
is the dissolution rate.They considered each ?ux line independently and suggested three different formulations to calculate the geometrical error at the termination point of the ?ux line on the cathode surface.For example,in the third formulation,the most promising of the three,the geometrical error at the end point of a ?ux line is given by
D error ?
l 2
eq r àq w T
eV w àV t Ttl eq r àq w T
where l is the length of the ?ux line,q w the calculated voltage gradient on the workpiece surface,and V w and V t the voltages on the tool and workpiece surfaces respectively.They repeatedly modi?ed the cathode shape until the error was acceptable.In many ways,this is similar to Hunt’s embedded technique [55]but it does not suffer from the disadvantage that analytical functions have to be used for representing the anode and cathode surfaces.Fig.13shows that even after 50iterations,the converged calculated cathode shape does not have the sharp corner that is present in the ‘‘exact’’tool shape,i.e.the tool which was used to generate the workpiece shape in the ?rst place.But this does not matter because the calculated converged cathode shape yielded the same work-piece shape as the exact tool.Bhattacharya et al.[6]used ?nite elements to determine the cathode shape iteratively but used a much simpler criterion:the difference between the required and calculated anode shapes.Das and Mitra [22]viewed it as a non-linear optimization problem in which they minimized
X
N i ?1
eq r àq w T2where N is the number of nodes on the anode surface.Zhou and Derby also considered it as an optimization problem but they minimized the difference between the required and calculated shape with respect to the coef?cients in the analytical expression de?ning the cathode shape [179].Like other researchers,they showed that the number of co-ef?cients in the analytical expression affects the accuracy to which the solution converges.
Chang and Hourng developed a comprehensive model for predicting the tool shape [14].Using Hunt’s embedded technique [55]and representing the tool and workpiece shapes as analytical functions,they computed the cathode shape considering the temperature increase due to Joule heating,void fraction and ?ow in two dimensions.
The only work to-date in predicting the inverse problem in 3-D is by Sun et al.[155]who claim to have obtained the tool shape in one iteration by calculating the theoretical value of the equilibrium gap using the cos u method.They were probably able to do this because the anode surface,in their case,was a turbine blade with a very gentle curvature.
4.Modelling of EDM processes 4.1.Generalized model for EDM processes
4.1.1.Discharge location
The generalized model for EDM processes is shown in Fig.14.In EDM,removal occurs at the discharge spot,where a tiny crater is generated.For each pulse,discharge occurs only at a single location where the dielectric strength is lowest.Hence,the EDM simulation
Fig.13.Determining the tool shape [116].
S.Hinduja,M.Kunieda /CIRP Annals -Manufacturing Technology 62(2013)775–797783
starts from the determination of discharge location.If discharge always occurs where the gap width is narrowest,the tool electrode shape can be replicated on the workpiece with a high accuracy.However,discharge does not occur deterministically at the narrowest gap.As shown in Fig.15,Schumacher [143]indicated that the gap space is not ?lled with the dielectric liquid,but a large fraction of the gap is occupied by bubbles and debris particles in consecutive pulse discharges.Hence,it is necessary to determine the discharge location taking into account the in?uence of debris particles and bubbles.
4.1.2.Simulation of EDM arc plasma
Simulation of arc plasma is a necessary pre-requisite for determining the temperature distribution and material removal in electrodes as it provides the necessary boundary conditions,i.e.heat input;diameter of heat source;pressure,and ?ow velocity.The heat input,which is the energy received from the arc plasma by the electrode,and diameter of the arc plasma exert a signi?cant in?uence on the temperature distribution in the electrode.As shown in Fig.15,discharge occurs in a narrow gap ?lled with dielectric liquid between parallel plane electrodes within a short duration of nano to micro seconds order.Plasma constituents are variable due to mass transfer from the evaporated dielectric and electrode materials.The conservation equations of mass,momen-tum and energy,Ohm’s law and Maxwell’s equations have to be solved considering moving boundaries.Moreover,considering that an analysis of arc plasma needs the temperature distribution and topography of both the anode and cathode surfaces as boundary conditions,the arc plasma simulation must be performed iteratively together with simulations of temperature distribution and material removal.Thus,a precise analysis of the EDM plasma is dif?cult.
4.1.3.Simulation of temperature distribution and material removal due to single discharge
Most of the simulations of temperature distribution and removal were performed assuming that the heat input and diameter of arc plasma are known.Recent developments of high speed video cameras have made it possible to observe and obtain the diameter of the plasma.The inverse problem solution based on the analysis of electrode temperature distribution combined with temperature measurement gives the heat input from the arc plasma.However,it is still dif?cult to calculate the removal amount because not all the molten region is removed [126,167,178].Furthermore,since the removal under the heat ?ux affects the temperature distribution in the discharge crater,simulations of the temperature distribution and material removal should be solved as a coupled problem.
4.1.4.Simulation of ?ow ?eld in gap
To determine the discharge location and to obtain the gap width distribution,the ?ow of debris particles and bubbles should be calculated.These ?ows are generated by the explosion of the bubble at the discharge spot,?ushing ?ow of dielectric liquid,and the jumping motion of the tool electrode.
4.1.
5.Simulation of geometry
The load of removal,which is the volume of workpiece removed per unit area of the tool,is not uniform on the tool electrode surface because it depends upon the depth of cut [18]and the surface inclination as shown in Fig.16(a).Hence,wear does not occur uniformly over the working surface.In addition,the load of removal depends also on the curvature of the tool electrode surface.A convex/concave unit area faces a greater/smaller workpiece surface than a ?at unit area.This non-uniform load removal can easily be reproduced by repeating the removal of both anode and cathode at the discharge location for each discharge,provided the location is determined correctly as described later.Provided the volume of removal per discharge in the simulation is correct,the obtained geometries should agree with the machining results.In reality however,the simulation accuracy is not satisfactorily high.This is because the bubbles existing prior to the discharge result in different plasma growth and discharge crater formation.Debris particles decrease the dielectric break-down strength of the gap,thereby increasing the gap width,which affects the plasma diameter,varying the removal amount [126,143].
4.1.6.Machine control
In EDM,discharge ignition is delayed with increasing gap width.On the other hand,if the gap width is too small,discharge occurs immediately.If the feed rate is too high,the tool electrode makes contact with the workpiece,leading to a short circuit.Hence,the tool electrode is fed based on the servo feed control as shown in Fig.17[151].It is therefore necessary to simulate the discharge delay time,which varies depending on factors such as the gap width,debris particle concentration and working surface area.
4.2.4, 4.2.5, 4.2.6
4.2.7, 4.2.8, 4.2.9, 4.2.10
4.3.1, 4.3.2, 4.3.3, 4.3.4
4.2.1, 4.2.2, 4.2.3
4.4
Fig.14.Generalized model for EDM processes.
Fig.15.Schematic view of discharge gap due to multiple discharge.
(a) Sinking EDM
(b) Wire EDM
Fig.16.Non-uniform load removal and gap width due to probabilistic nature of discharge.
S.Hinduja,M.Kunieda /CIRP Annals -Manufacturing Technology 62(2013)775–797
784
4.1.7.Simpli?ed and partial models
The idea that EDM can be simulated accurately by repeating all the steps in the generalized model is not realistic or possible.It is rather useful to simulate a part of the generalized model for understanding the gap phenomena,optimizing the machine control,evaluating the machining accuracy and for process planning.
The simulation method of time-sequential repetition of the ?owchart which models a single pulse discharge needs a long calculation time,longer than the actual machining.Hence, Tricarico et al.[164]repeatedly calculated the thickness of the removal layer on the tool electrode and workpiece for every small feed step taking into account the dependence of the material removal rate on the local gap width using surface models for the tool electrode and workpiece geometries as shown in Fig.18.In this simulation however,it is dif?cult to reproduce the gap phenomena for each discharge.
4.1.8.Geometry considerations
Geometrical simulation in EDM needs consideration of tool electrode wear.In the case of surface models shown in Fig.18, distance between adjacent nodal points decreases on the tool electrode due to wear,whilst nodal points on the lateral surface of the machine hole move further apart.Calculation of the direction of the normal at each nodal point is susceptible to signi?cant error. Hence,the simulation must be interrupted frequently for re-meshing to avoid topological instability[108].In contrast,the use of voxel models as shown in Fig.19,makes re-meshing unnecessary;also the singularity problem at sharp edges can be eliminated[109].4.2.Modelling of a single discharge
4.2.1.Probability of discharge
To ignite an electric discharge in clean oil at an open voltage of 100V,the gap width must be less than several microns.In the case of EDM however,since the dielectric liquid is contaminated with electrically conductive debris particles,whose average diameter is even more than one third of the gap width[176],discharge can occur at gap widths of tens of micrometres or more [8,102,143,154].This fact sometimes reduces the replicating accuracy due to the uneven distribution of the debris particles. On the other hand,the extended gap width is favourable for easy gap control because it is dif?cult to keep the gap width constant at several microns.
In order to study the in?uence of debris particles,a debris particle with a diameter of5m m was placed in a gap of20m m as shown in Fig.20[89].If it were true that the discharge occurs at a point where the gap is shortest,then the discharge should occur at the point where the debris particle is placed.However,experi-mental results showed that in most cases,the discharge crater was not generated at the point where the debris particle was placed. This is because the probability of discharge in a certain area is obtained from the product of the probability of discharge per unit area and its surface area(area effect).Since the projected area of the debris particle in the direction normal to the electrode surface is negligible compared with the area of the electrodes (50mm?50mm),the probability that discharge does not occur on the debris is greater than the probability that the discharge occurs on the debris.This result indicates that the discharge location is determined in a probabilistic way and not deterministic.
Another example which indicates the importance of the probabilistic model is already shown in Fig.16.In sinking EDM, when the tool electrode is fed by a distance during a certain time, the duty of removal volume per unit area on the tool electrode is greater at a curved surface than that at a?at surface.The duty on an inclined surface is lighter than that on a horizontal surface. Therefore,the discharge frequency per unit area must be higher on a convex surface,and lower on an inclined surface compared with that on?at and horizontal surfaces,respectively.To satisfy these conditions,the gap width must be smaller on a convex surface,and greater on an inclined surface than the corresponding gaps on a horizontal surface.In wire EDM,the gap width distribution around the wire is not uniform as shown in Fig.16(b),which would not have happened if the discharge location were to be determined deterministically from the gap width distribution.
4.2.2.Discharge delay time
The discharge delay time depends upon the local values of the gap width and concentration of debris particles at the discharge location.Therefore,the discharge delay time can be used as a measure to determine the discharge location.Moreover in reality, the tool electrode feed is controlled based on the discharge delay time.Hence,modelling of discharge delay time is the key to EDM simulation.
Morimoto and Kunieda[109]quanti?ed empirically the average discharge delay time t d,ave[m s]as a function of the gap width,gap[mm],concentration of debris particles,conc[mm3/
Fig.17.Principle of servo feed control.
Fig.18.Surface model for geometrical simulation[164].
Fig.19.Voxel models for sinking EDM[109].
Fig.20.Probabilistic nature of discharge location.
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797785
mm 3],machining area,area [mm 2],and debris diameter,r [mm]under the following machining conditions:copper anode,carbon steel cathode,discharge current of 3A,discharge duration of 300m s and open circuit voltage of 120V.t d ;ave ?8:2?10
12
á
ga p 8:8?r 2:9
area 1:2?conc 1:6
(25)
The machining area is a factor because the probability of discharge within a speci?c area is proportional to the ratio of its surface area to the whole discharge area as described in the previous section.t d,ave was obtained experimentally based on the fact that the discharge delay time t d conforms to the exponential distribution obtained from the Laue plot method [5,8].The Laue plot shows the percentage of electric insulation that does not break down until time t after the supply of a pulse voltage [58].Suppose a single pulse discharge was generated N times and the number of discharges which did not occur until time t was n ,then n as a percentage of N can be expressed by:
n ?exp àt d ;ave
(26)where t d,ave the average value of t d ,is given by [107]:t d ;ave ?
S t d
N
(27)
The Laue plot enables easy evaluation of the discharge delay time,which normally has a large scatter,from the slope t d,ave as shown in Fig.21.
4.2.3.Theory of determination of discharge location
To determine the discharge location,the discharge delay time of each voxel of the tool electrode is calculated probabilistically using Eqs.(25)and (26).In Eq.(25),area is set to the area of each voxel.conc and gap are the concentration of debris particles and gap width at each voxel whose values are recalculated for every discharge cycle.conc is de?ned as the volume of debris particles divided by the volume of the gap voxels close to the voxel of the tool electrode.The last variable r ,debris diameter,can be obtained from machining experiments.Then,using Eq.(26),the Laue plot can be drawn for every voxel as shown by the green line in Fig.21.Note that the slope of the Laue plot for each voxel is signi?cantly small when compared to that obtained from actual machining,because the area of each voxel is nearly equal to the discharge crater size which is much smaller in size than the tool electrode surface.Thus t d,ave calculated for each voxel is much longer than that in actual machining.It must be emphasized that t d of each voxel should be determined probabilistically.Substituting a random number from 1to 100%for n /N along the vertical axis in Fig.21,t d can be calculated based on the Laue plot of each voxel as indicated by the blue arrow.This operation is veri?ed by the theory of statistics which shows that,to generate a value of a random variable X having a distribution function U =F (X ),it
suf?ces to generate a value of a random variable U that is uniformly
distributed [107].Since the number of voxels existing over the machining surface is several thousands to several tens of thousands,t d may even be as small as a few micro-seconds.Therefore,if the voxel with the minimum t d is determined as the discharge location,the same orders of t d as measured in actual machining processes can be reproduced.Consequently,the greater the machining area,the shorter will be the discharge delay time.Thus,the area effect is successfully reproduced in the simulation.4.2.4.EDM arc plasma
Since erosion in the EDM process is made by accumulation of removal due to single discharges,the simulation model for heat conduction in electrodes is important.With the advent of powerful computers and numerical analysis methods,it is nowadays not dif?cult to take into account the time-dependent radius of the circular heat source and time-dependent heat ?ux,or temperature dependence of thermo-physical properties of electrode materials.However,boundary conditions such as arc plasma diameter,which is time-dependent,and ratios of power distributed into the anode and cathode have not yet been obtained theoretically,although they exert a signi?cant in?uence on the calculation results.This is because simulation of EDM arc plasma is extremely dif?cult.
Hayakawa et al.[46,47]?rst conducted magneto-hydrody-namic analysis of a DC arc in air between parallel plane electrodes under the gap conditions used in EDM.They assumed that the species in high-temperature air which includes copper electrode vapour are N 2,O 2,NO,N,O,Cu,NO +,N +,O +,Cu +,N 2+,O 2+,Cu 2+and electrons.Considering the temperature dependence of the thermo-physical properties of the plasma,the electromagnetic ?eld,temperature,pressure and velocity distributions were calculated for the regions including both the electrodes and discharge gap.The conservation equations of mass,momentum and energy,Ohm’s law and Maxwell’s equations were solved.The energy equation included Joule heating,conduction,convection,and radiation terms.Fig.22shows an example of the temperature ?elds.It was found that most of the discharge power is distributed in the electrodes,and heat transfer due to convection and radiation is negligible.Although this result is in agreement with the experimental results of a DC arc [48],the arc which they analysed was not in dielectric liquid but in air.The arc was not in transient but in steady state,and removal of the electrodes was not taken into consideration.Furthermore,the arc was assumed to be in thermo-equilibrium,and the equations of motion of the three species:electrons;ions and neutral particles,were not solved separately.Thus,the gap phenomena were symmetrical between the anode and cathode thereby making the fractions of energy distributed to the anode and cathode equal.Consequently,it was dif?cult to obtain the difference in the energy distribution between the anode and cathode.
4.2.
5.Energy distribution
When copper was used for both the anode and cathode,the anode removal amount was greater than that at the cathode for a
Fig.21.Evaluation of discharge delay time using Laue plot [109].
Fig.22.Magnetohydrodynamics analysis of EDM arc coupled with heat conduction analysis in electrodes [46].
S.Hinduja,M.Kunieda /CIRP Annals -Manufacturing Technology 62(2013)775–797
786
discharge duration shorter than20m s,while it was smaller with a discharge duration longer than20m s[170].Motoki and Hashi-guchi[110],and Van Dijck[167]explained that such a phenomenon is caused by the variation of the energies distributed into the anode and cathode with the discharge duration,on the basis of the T–F electron emission theory[99].Other papers suggest that this is due to the small mass of electrons compared to ions,resulting in quicker impingement of electrons into the anode surface than the slow ions,which arrive at the cathode surface later [24].However,there is no scienti?c evidence that the delay of ions can be in the order of nanoseconds or more.
On the other hand,Koenig et al.[77]measured the energy distribution by measuring the temperature of the electrodes and dielectric?uid in consecutive pulse discharges.Xia et al.[170] measured the energy distribution in a single discharge by comparing the measured temperatures of the foil electrodes with the calculated results obtained under the assumed ratio of the energy distributed in electrodes,using a FD model and the experimental setup shown in Fig.23.When the calculated and measured temperatures were in agreement,they found that the energy distributed to the anode and cathode is about40–48%and 25–34%respectively.They reported that the energy distributed to the anode is always greater than that to the cathode and is independent of the discharge duration.Recently,Zahiruddin et al. [178]measured the energy distribution in micro EDM.Fig.24 shows the overall results of energy distribution ratio versus discharge durations.The energy distribution into the anode is consistently greater than that into the cathode regardless of the discharge duration.
Thus,it is obvious that the difference in the volumes removed between the anode and cathode cannot be fully explained by the energy distribution.Motoki et al.[111],Ikai et al.[56]and Mohri et al.[105]explained that the amount of electrode wear is reduced due to the protective effects of carbon layer which is formed on the anode surface.Since a thicker carbon layer is generated when the discharge duration is long,tool(anode)wear is smaller than that on the workpiece(cathode)even though the energy distribution to the anode is more.4.2.6.Plasma diameter
Another important boundary condition is the plasma diameter. There are many papers in which a point heat source is assumed at the cathode spot based on the researches on gas and vacuum discharges ignited with gap widths which are signi?cantly greater than the plasma diameter.DiBitonto et al.[34]and Patel et al.[128] obtained the energy distribution ratios into the anode and cathode as8%and18%,respectively,assuming the point cathode source model.However,there have been no photographs taken to verify the point source model in the narrow EDM gap yet.Another simpli?ed model,which is widely used without evidence,assumes that the plasma diameter is equal to the diameter of discharge craters[180].
To obtain an evidence-based diameter,Snoeys and Van Dijck [152,153,167]systematically analysed the temperature distribu-tion assuming a circular heat source with time-dependent radius and time-dependent heat?ux on a semi-in?nite cylinder.The heat source growth function was obtained from an iterative calculation of the heat conduction equation by assuming that the temperature at the centre of the heat source corresponds to the metal boiling temperature at a pressure equal to an average pressure in the gas bubble calculated from the thermo-dynamical model.Good agreement was found between the calculated melting point isothermals and those measured from pictures of the cross-sections of the craters cut perpendicular to the surface.
Recent developments of high-speed video cameras have facilitated measurement of the arc plasma diameter in EDM [33,78].The radial temperature distribution in the arc plasma was measured by optical emission spectroscopy and the plasma area where discharge current?ows was determined based on the fact that the degree of ionization depends on the temperature[78].The measured diameter of arc plasma generated in air was0.5mm,?ve times larger than the crater diameter as shown in Fig.25.The plasma completed expanding within2m s after dielectric break-down,and thereafter,its diameter remained constant during discharge.Although the expansion rate was lower when the arc plasma was generated in a dielectric liquid,its diameter was still greater than the crater diameter.
4.2.7.Simulation of removal
4.2.7.1.Temperature rise in electrode.Given the energy distribution ratio,and time-dependent radius of the circular heat source as boundary conditions,the following heat diffusion equation can be solved to obtain the temperature distribution in the electrodes. @
@x k
@T
@x
t
@
@y k
@T
@y
t
@
@z k
@T
@z
t˙q?r c p
@T
@t(28)
where˙q is the rate at which energy is generated per unit volume of the electrode due to Joule heating.By integrating the solution for a point source which is liberated instantaneously at a given point and time with regard to appropriate space and time variables,one can obtain solutions for instantaneous and continuous sources of any spatial con?guration[12].The analytical solutions of the temperature distributions at time t produced by stationary Gaussian energy distribution sources in a semi-in?nite solid are given by Pittaway[133].
Fig.23.Method to obtain energy distributed into electrode in single pulse discharge [170].
Fig.24.Energy distribution into anode and cathode[178].
Fig.25.High-speed video frames of EDM arc plasma(i e:23A,t e:80m s,u i:280V) [78].
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797787
Nowadays,with the development of powerful computers and numerical methods,it is no longer dif?cult to take into account the temperature dependence of thermophysical properties of elec-trode materials and latent heat of melting and vaporization. However,it should be emphasized that correct temperature distributions can be obtained only if correct boundary conditions are used,i.e.energy distribution and plasma diameter.Although there are many papers dealing with temperature analysis,in most cases there is insuf?cient discussion about correct boundary conditions.
Zingerman[181]states that Joule heating does not play a substantial role in the machining of metallic materials,because electrical potential drops in metallic workpieces are negligibly small.According to the calculated results of Rich[138],however, for high-resistivity metals,e.g.,Hg,Sb,Fe,and Bi,Joule heating is comparable to the energy input from the arc plasma.Saeki et al. [141]found that the removal of high-electric-resistivity material
(Si3N4–30wt%SiC)in a single discharge can be greater than that of steel due to Joule heating.
4.2.7.2.Mechanism of removal.Calculation of the temperature distribution in the electrode only is insuf?cient to obtain the volume of removal due to a single discharge.Simulations of crater formation and eruption of molten pool are indispensable.There is no doubt that the material can be ejected when its temperature exceeds the boiling point.However,there are few papers dealing with the removal model of the melting area.Zolotykh[182],Snoeys and Van Dijck[152,153,167],and Tao et al.[157]reported that most of the metal removal occurs due to boiling of the superheated molten mass in the crater at the end of the discharge because boiling of that superheated metal is prevented by the bubble pressure during the discharge.Van Dijck calculated the volume of the region inside the normal boiling point isothermal surface,at the end of the discharge.The calculated volume agreed well with the measured material removed per pulse.It was also found that the material removal ef?ciency,which was de?ned as the ratio of the ejected to melted volume,was only1–10%.However, observation of the gap of single pulse discharge in dielectric liquid using X-ray showed that85%of material removal occurs during the discharge duration[38].High-speed camera images of?ying debris scattered from discharge points showed that material removal also occurs during the discharge[45].Furthermore,Yoshida et al.[176] found that the metal removed per pulse in air is almost equal to that in liquid when the discharge duration is longer than100m s, indicating that metal removal can occur without a sharp drop in bubble pressure.Thus,development of computation models which can simulate the material removal phenomena is eagerly awaited.
4.2.7.3.Advanced analysis methods for simulation of crater formation and material removal.In some cases,molecular dynamics(MD) simulation can be very helpful to analyse the phenomena which are dif?cult to model using conventional thermo-hydrodynamics simulations.Shimada et al.[146]simulated the formation of a thin electrode in a single discharge.Yang et al.[173,174]analysed the forming mechanism of a discharge crater ignited in vacuum as shown in Fig.26.Both the space and time domains which can be handled are unrealistically small compared with those used in EDM.Nevertheless,Fig.26clearly demonstrates the evidence of the superheating theory,indicating that removal can occur at any time when the pressure inside the molten pool exceeds the pressure outside[174].Tao et al.[157]simulated realistic crater morphology,such as build-up crest and bulged bottom,using FLUENT,a commercial computational?uid dynamics(CFD) software.To incorporate both molten and solid materials together in the same analysis domain,the volume of fraction(VOF)method was used,where both liquid and solid phases are modelled using liquid type cells.
4.2.8.Bubble
Electrode materials and dielectric liquid evaporate,molecules are dissociated,and atoms ionized,resulting in a rapid expansion of a bubble.Fig.27shows side views of a bubble generated by a single discharge in deionized water and EDM oil.The bubble and dielectric liquid are analogous to the spring and mass oscillation system respectively[36,57].Starting from the initial condition, where the bubble is compressed in a small volume,the dielectric liquid is accelerated radially.At the moment the pressure inside the bubble equals the atmospheric pressure,the kinetic energy peaks.Hence,the bubble continues expanding.The diameter of the bubble peaks when all the kinetic energy is transferred to the potential energy of the bubble.The diameter of the bubble reaches several millimetres,several tens of times greater than the gap width.Thereafter the bubble starts contracting until it is compressed to its initial diameter.In reality,the viscosity of the liquid causes damping in the oscillation.After discharge,ions and electrons are recombined,and the evaporated atoms and molecules are solidi?ed or condensed to form debris particles or dielectric liquid,but gases such as hydrogen and methane which are generated by the dissociation of the working oil are left to form
a bubble.
4.2.9.Reaction force acting on electrodes
Fig.28shows a model to calculate the bubble oscillation in the gap between parallel plane electrodes[36].The reaction force acting on the tool electrode can be calculated by integrating both the pressure in the bubble and that in the dielectric liquid over the working surface[86].The reaction force in Fig.29was measured by Kunieda et al.[91]using the Split Hopkinson Bar method[172].At the initial state in which the bubble is compressed,the force reaches its peak value.With the expansion of the bubble,the bubble pressure decreases and so does the force.The force even becomes negative because the bubble continues to expand even after the bubble pressure falls below the atmospheric pressure due to the inertia of the dielectric liquid.Since the natural frequencies of EDM machine structures are much lower than the frequency components included in the reaction force waveform,the gap width cannot respond to the change in force.Furthermore,the reaction force in a series of pulse discharges decreases with time while the working gap is?lled with bubbles generated in
Fig.26.Molecular dynamics simulation of material removal[173].
Fig.27.Side view of bubble oscillation(i e:20A,t e:100m s,gap width:0.1mm, anode:Cu f5mm,cathode:Cu f5mm).
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797 788
consecutive discharges[86,91].Hence in?uence of the reaction force caused by bubble pressure is negligibly small in sinking EDM. However,it causes vibration and de?ection of electrodes in the case of wire EDM[106,123]and micro EDM[72].
4.2.10.Residual stress
The temperature rise due to discharge results in a compressive thermal stress on the electrode surface until the stress is released by material yielding at a high temperature.The stress becomes nearly zero when the material becomes molten.After discharge, due to the temperature decrease,the molten material re-solidi?es and starts shrinking,generating a tensile residual stress on the surface.This model was analysed using FEM[21,175]and MD simulation[174].
4.3.Modelling of gap phenomena due to multiple discharges
4.3.1.Difference from single discharge
Generation of bubbles causes?uid?ow as described in the next section.In addition,pressure in the bubble generated by a discharge is signi?cantly in?uenced by bubbles which are already present[91],resulting in different plasma growth and discharge crater formation.Observation of gap phenomena through a transparent electrode clearly showed that discharge occurs even in bubbles,and that craters formed in bubbles are different from those in the liquid[70].Since debris particles weaken the dielectric breakdown strength of the gap,the gap width is increased.Thus, the plasma diameter is increased,thereby varying the removal amount.Therefore,diameters of craters and debris particles in consecutive discharges are not equal for each discharge and show different intermediate values between single pulse discharges in the gap?lled with gas and liquid[156].
Furthermore,the uneven concentration of debris particles results in a non-uniform gap width,deteriorating the machining accuracy.The cooling effect of the liquid may also be affected by the bubbles.Concentration of discharge locations may result in locally elevated temperature and accumulation of debris particles, leading to unstable discharge.Consequently,results of multiple discharges are not equal to those obtained from a superposition of the results obtained from a single pulse discharge[126].Therefore, simulation of gap phenomena considering multiple discharges is necessary.
4.3.2.Bubbles in consecutive discharges
The working gap is mostly occupied by bubbles although both electrodes are submerged in dielectric liquid[70].The?ow?eld in the gap was calculated by superimposing the dielectric?ow caused by the bubble generated at a discharge point and the radial?ow of dielectric?ushing supplied from the centre hole[156].Assuming that discharge occurs at random on the electrode surface,an animation of the dielectric?ow?eld was made as shown in Fig.30.
4.3.3.Debris movement
Fig.31shows an electrostatic?eld distorted by a negatively charged debris particle[92].The measured gap width is occasionally larger than100m m,whereas the measured average diameter of debris particles is25m m or less.This result can be explained in terms of the electrophoresis movement of debris particles in the direction perpendicular to the electrode surfaces, which was observed by Suda et al.[154],Bommeli et al.[8],and Schumacher[143].They also observed that some particles are linked in series to form chains parallel to the electric?elds.The speed of the electrophoresis-induced motion was calculated to be 0.136m/s,high enough for a particle to change its position during the discharge delay time in the narrow discharge gap[92].
4.3.4.Flushing
To maintain stable machining and to obtain narrower and more uniform gap width,it is critical to?ush bubbles and debris particles and cool the working gap in order to prevent the localization and concentration of discharge locations.Pressure or suction?ushing through holes in the electrode or workpiece remains one of the most ef?cient?ushing methods especially if these holes have to be provided anyway or do not harm the workpiece.Both in pressure and suction?ushing,one can observe lower electrode wear and larger gap width at the outlet point in comparison with the inlet point[76].Normally,the tool electrode is lifted periodically to replace the contaminated dielectric?uid in the gap with fresh?uid.
Fig.29.Calculated force acting on electrodes due to bubble oscillation(i e:45A,t e: 150m s,gap width:150m m)[91].
Fig.30.Bubbles in gap after300times of discharge(110ms after start of machining, i e:19A,t e:250m s,t o:50m s,diameter of copper electrode:30mm)[156].
Fig.31.Calculated equipotential lines generated by negatively charged particle [92].
Fig.28.Pressure distribution in bubble generated by single pulse discharge in gap
?lled with dielectric liquid.
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797789
Cetin et al.[13]calculated the three-dimensional?uid?ow in the gap considering the suspended debris particles and obtained the relationship between the tool lifting velocity and height and the ?ushing capability.The working gap can be?ushed by a fresh dielectric?uid jetted from nozzles placed adjacent to the discharge gap.Masuzawa et al.[103]however demonstrated through hydrodynamic analysis that jetting of dielectric?uid merely from one direction causes increased density of debris particles in the downstream region,resulting in uneven distribution of gap width deteriorating the machining accuracy.
Okada et al.[124]analysed the three-dimensional?uid?ow around the wire electrode when the workpiece is being cut by WEDM as shown in Fig.32.They investigated the in?uence of the position of the dielectric jet nozzles relative to the upper and lower surfaces of the workpiece and?ushing?ow rate on the ability to ?ush the debris particles out of the gap.
From these?uid dynamics analyses,debris concentration can be calculated quantitatively.However,the relationship between the debris concentration and gap width cannot be determined unambiguously.The gap width cannot be de?ned uniquely by the debris concentration because the gap width is time dependent.The gap width increases with the passage of time even if the tool electrode is not fed,because discharge can still be ignited even if the gap width is increased,although the probability of discharge decreases.Hence,a probabilistic approach should be used to obtain the gap width distribution from the debris concentration.
4.4.Machine control
Given the shape of the tool electrode,the workpiece shape can be obtained precisely,only if the gap width distribution is known. Since the gap width is determined by the servo feed control,a precise geometrical simulation cannot be performed without an accurate model which can reproduce the feed control of the tool electrode.Since t d can be obtained for each pulse discharge using the method described in Section4.2.2,the moving average of the gap voltage can be calculated from the t d values sampled in a certain period of time.Thus the tool electrode can be fed forward or backward depending on the difference between the average gap voltage and servo reference voltage as in actual servo-control.
4.5.Empirical data required
Although the EDM process can be modelled as described in the previous sections,certain data which can only be obtained empirically is required.For example,a pre-requisite for the determination of the discharge location is the evaluation of the average discharge delay time,as expressed by Eq.(25).This equation has to be determined experimentally for every machining condition:open voltage;dielectric liquid;and electrode materials.
Similarly,the simulation of arc plasma is not possible and therefore the energy distribution to the electrodes and plasma diameter are not https://www.sodocs.net/doc/8f4222461.html,ually researchers assume values for these two parameters;instead they should be obtained from experiments and then used as boundary conditions for determin-ing the temperature distribution in electrodes.
Another importation parameter which has to be empirically determined is the removal ef?ciency.Its value has an important bearing,during the geometric simulation,on the removal volumes of anode and cathode per single pulse discharge.Even if the pulse conditions,materials used for the electrodes,and dielectric liquid are the same,the removal ef?ciency will be different depending on the working surface area and gap width[164]as shown in Fig.33. Other data which must be supplied are the volume of the bubble generated by a single discharge,debris particle size,and thermo-physical properties of the electrode materials.
5.Modelling examples in ECM and EDM
5.1.EC milling
The advantage of EC milling is that the required shape is obtained by moving a tool of simple shape along the three orthogonal axes.However,calculating the current density distribution becomes more dif?cult as the inter-electrode gap is now three-dimensional.Kozak et https://www.sodocs.net/doc/8f4222461.html,ed an analytical method in which they assumed a linear variation of the current through the thickness of the inter-electrode gap[80].This assumption helped them to calculate the current density and simulate the machining of a free-form surface with a spherical tool.However,their model did not consider stray machining which is a very important consideration especially in EC milling.Domain-based methods like ?nite elements do not lend themselves readily because the mesh would have to be re-generated with3-D elements after each time step.
Pattavanitch simulated and actually milled several components varying in complexity from a simple slot[130]to the diamond-shaped pocket with a human-shaped protrusion shown in Fig.34 [129].He subsequently modelled the machining of these components with the BEM.Whilst the accuracy obtained from the models was reasonably good,he faced several dif?culties in preparing data,two of which are mentioned below.
(i)Since the pocket was machined in several axial passes,the
shape of the human?gure had to be expanded and the boundary of the pocket contracted by the cumulative overcut.
Commercial CAM programs do not have this facility and special routines had to be developed.
(ii)A preliminary investigation showed that the time step should not exceed0.1s.Since the workpiece was machined with a feed
Fig.32.Simulation of?iud?ow in wire EDMed kerf[124].
Fig.33.Material removal rate V w and volumetric tool wear rate V E versus gap width [164].
Fig.34.Pocket with a human-shaped protrusion[129].
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797 790
rate of9mm/min,it meant that in one time step,only
0.015mm of the tool path could be considered.The tool path
was608mm long and therefore it had to be sub-divided into well over40,000segments,with every few segments requiring
a mesh regeneration.Tool paths were generated using
MASTERCAM,a commercial software normally used for generating tool paths for turned and milled components.
5.2.Graphical and mathematical simulation
The accuracy of FE and BE models depends,in addition to the mesh density and quality,element type and time step,on the input values such as polarization voltage and current ef?ciency.Unless the latter are accurately speci?ed,there will always be discrepancy between the experimental and theoretically predicted results.To ensure that potential models generate accurate results,one needs to know,for different electrolyte concentrations and pulse times, the relationships between:
(i)the polarization voltage and current density,and
(ii)the current ef?ciency and current density.
Altena derived the above relationships experimentally[4] which have now been used by several modellers.In addition to the above,Altena determined other relationships such as that between the process voltage and frontal gap.With the help of these relationships,accurate graphical and mathematical models for the process were built by DeSilva et al.[35].These models were used not to predict the current density but instead to determine the tool dimensions and the side gaps so that the slots were not unduly tapered.For example,to machine the300m m wide slots in the shaver head in Fig.35,the width of each of the90radially-oriented lamellae was determined as120m m and it was found that the actual width of the machined slots was measured as285m m[125].
5.3.EC turning
Although researchers have EC machined axi-symmetric com-ponents,very few have attempted to model the process.Ma et al. [101]proposed a mathematical model based on Faraday’s and Ohm’s laws to calculate the amount of material removed in pulsed EC-turning.Their model is a function of the cutting velocity,width of the tool and pulse-time.However,their model is of limited use as it is one-dimensional in nature.Pattavanitch et al.[130]have shown that a3D model of EC turning is feasible;they reported the development of a BE model to simulate the machining of the disc shown in Fig.36(a).Fig.36(b)shows part of the BE model-it does not show the bounding virtual surfaces.Since the length of the ring was equal to the depth of the tool,it was not necessary to feed the tool in the longitudinal direction.Also no radial feed was involved. The purpose of modelling was to determine the time required to reduce the diameter of the thin disc by51m m which was determined as90s and subsequently veri?ed experimentally.
5.4.Multi-ion and transport model for ECM
The importance of modelling the electrolyte?ow and thermal effects was recently demonstrated by Deconinck et al.[30,31].Fig.37shows the anode shape(black curve)obtained using their multi-ion transport and reaction model which takes into account the electrolyte?ow and heat energy generated in the electrolyte and double layer.The shape of the anode curve is asymmetrical (black curve)and Deconinck et al.[32]attributed this asymmetry to the increased temperature of the electrolyte in the upstream region.The temperature increase was due to the formation of eddies,preventing the heat being ef?ciently being transported away.Deconinck et al.also suggested that water depletion may be another factor causing the asymmetry.Water depletion makes oxygen evolution dif?cult,causing a greater portion of the current to be spent on metal dissolution,thus increasing the current ef?ciency.It demonstrates that modelling of?ow and internal energy has a signi?cant effect on the computed anode shape.
https://www.sodocs.net/doc/8f4222461.html,puting time in ECM
No review on modelling would be complete without a few words about computing time which can become prohibitive especially since several coupled models are involved.The computing time depends,amongst other factors,on the ef?ciency of the algorithm to solve the resulting simultaneous equations.For this,the FAST Multipole Method suggested by Greengard and Rohklin[39]should be used to make the stiffness matrix sparse and then a GMRES technique used to solve the resulting equations [25].
In the case of EC milling,nodes lying on the workpiece surface should be classi?ed as being either‘‘active’’or‘‘passive’’[129].The current density on the workpiece surface has positive values only in and around the projected area of the tool(yellow nodes)and a little distance away from the tool,its magnitude becomes zero.In Fig.38,the?ux density at the black nodes(which are on the workpiece surface but are far away from the projected area),is virtually zero;such nodes are said to be‘‘passive’’since both the potential and current density are known.Such nodes should not be considered and this will decrease the computing time substantially.
Fig.35.Machining slots in the Philips shaver head[125].
Fig.36.BE model for a thin ring[130].
Fig.37.Temperature distribution using multi-ion model[32].
Fig.38.Active and passive nodes[129].
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797791
5.6.Sinking EDM
Kunieda et al.[88,109]developed a sinking EDM simulation using the discharge location search algorithm shown in Fig.39 which was obtained by simplifying the generalized model(Fig.14). In the?rst step,the discharge location was determined probabil-istically as described in Section4.2.3.In the second step,voxels of both the tool electrode and workpiece at the discharge location were converted to gap voxels,reproducing the material removal. Since the simulations of arc plasma,temperature distribution,and removal were not performed,the removal volumes of tool electrode and workpiece per discharge were obtained from experiments.In the third step,debris particles newly generated by the discharge were located evenly on the periphery of the bubble generated at the same time,and debris particles present prior to the discharge were re-located due to the dielectric liquid ?ow caused by the bubble expansion.In the servo feeding step,the average gap voltage was calculated from the moving average of the discharge delay time of consecutive pulse discharges.Then the feeding speed of tool electrode was determined depending on the difference between the average gap voltage and servo reference voltage as is done on actual EDM machines.These steps were repeated until the machining depth reached the pre-set value.
A cylindrical steel workpiece was machined using a copper tool electrode with the same diameter.Dielectric liquid was supplied into the gap from the centre of the tool electrode.Fig.40shows distribution of debris particles simulated for two different?ushing ?ow rates.Fig.41compares the calculated gap width distributions with experimental results.The gap width increases in the radial direction,because debris particles are transferred with the?ushing ?ow of dielectric liquid.The discharge location search algorithm combined with heat transfer analysis realizes the simulation of surface temperature distribution in consecutive discharges [60,66,171].
5.7.Reverse simulation in EDM
Contrary to the forward simulation where the tool electrode is given,it is practically more important to develop a reverse simulation for obtaining the tool electrode shape with which the target workpiece shape can be machined precisely[87](see Fig.42).The reverse simulation was conducted using the same algorithm as that developed for the forward simulation in Fig.39 assuming that the tool electrode is machined by the workpiece having the same initial shape as the target workpiece shape.The data of the removal volumes per discharge obtained from experiments were switched between tool electrode and workpiece in the reverse simulation.
https://www.sodocs.net/doc/8f4222461.html,ling EDM
The machining of a groove by EDM was simulated considering wear of the tool electrode[177].In-process measurement of the tool electrode wear enabled the adaptive compensation of the depth of cut to machine stepped grooves to a high accuracy[7]. Using the discharge location search algorithm,a three-dimensional geometric simulation of micro-EDM milling was developed as shown in Fig.43[49].
5.9.Wire EDM
5.9.1.Forces acting on wire electrode
In WEDM,there are four kinds of forces acting on the wire electrode[32,127]:discharge reaction force;electrostatic force; hydrodynamic force and electromagnetic force(see Fig.44).The discharge reaction force is caused by the rapid expansion of a bubble at the discharge spot during the discharge.The electrostatic force acts mostly when an open voltage is applied between the wire and workpiece during ignition delay time.The electromag-netic force acts on the wire during the discharge and can be calculated from the area integral of the vector product between the current density and magnetic?ux density in the wire.The
Tool electrode voxels Fig.39.Discharge location searching algorithm[109].
Fig.40.Simulation of debris movement under?ushing?ow[109].
Fig.41.Gap width distribution in radial position.
Fig.43.Simulation of micro EDM milling[49].
Fig.42.Reverse simulation.
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797 792
hydrodynamic force is the drag force generated by the ?ow of
dielectric ?uid.These forces cause vibration and de?ection of the wire,thereby lowering the machining accuracy,speed,and stability [32,69].
Obara et al.[123]and Mohri et al.[106]obtained the reaction force from solutions of the inverse problem in which wire vibrations calculated using assumed values of the force were compared with measured ones.They used thin workpieces to exclude the in?uences of the electrostatic and electromagnetic forces.Obara et al.[121]measured the change in the resultant force involving the reaction force,electromagnetic force,and electrostatic force with various discharge frequencies.They then obtained the electrostatic force by extrapolating the resultant force to the limit of zero discharge frequency,because both the reaction force and electromagnetic force are zero when the discharge frequency is zero.As for the hydrodynamic force,Kuriyama et al.[93]conducted a CFD analysis of the force and investigated the in?uence of the jet ?ushing conditions on the wire de?ection.
Regarding the electromagnetic force,Tomura and Kunieda [161]developed a FE model to calculate the electromagnetic ?eld and the electromagnetic forces as shown in Fig.45.Tomura and Kunieda [162]found that,using a workpiece 40mm thick,the reaction force is larger than the electrostatic force with large discharge energy,while their magnitudes are reversed with decreased discharge energy.The in?uence of the electromagnetic force on the wire vibration is not negligibly small under rough cutting conditions,especially with higher discharge frequencies and larger workpiece thicknesses.
5.9.2.Simulation of wire vibration and de?ection
Obara et al.[122],Han et al.[41],and Tomura et al.[163]developed programs for WEDM simulation.The simulation shown in Fig.46is based on the repetition of the following routine:calculation of wire vibration considering the forces applied to the wire,determination of the discharge location considering the gap width between the wire and workpiece,and removal of workpiece at the discharge location.The electromagnetic force and hydro-dynamic force were ignored,and the in?uence of debris particles was not considered.The geometrical simulation error was less than 1.5m m. 5.9.3.Wire breakage
Fig.47shows the temperature distribution along the wire electrode obtained from a heat transfer analysis using FDM [65,90].Zone 3indicates the part of the wire electrode where discharge occurs.Zones 2and 4show the part from the upper and lower feeding points to the upper and lower surfaces of the workpiece,respectively.Although the average temperature is around the boiling point of water,which is used as the working ?uid,the temperature at the point where the preceding discharge occurred is signi?cantly high so that the tensile strength of the wire weakens at this point.Consideration of this has led to the development of adaptive control systems in which the pulse energy is reduced or stopped based on the distribution of discharge locations measured in process [95,120,147].
5.9.4.Optimization of wire electrode composition
Fine wire electrodes with diameter of 30m m or less are becoming popular.To resist the tension force,high tensile strength materials such as tungsten or molybdenum are used.However,since these materials are rare metals,steel wires coated with brass or zinc are being developed [142].Another requirement of the wire electrode is low impedance at the frequency components involved in the discharge current waveform.Thus,the electromagnetic ?eld shown in Fig.45was analysed to investigate the in?uence of the wire electrode and workpiece materials on the discharge current [40].
6.Conclusions and future work
This paper has described the models to simulate the ECM and EDM processes.In the case of ECM,commercial systems based on the potential model have become available.In the case of EDM,the picture of the process is not completely clear probably because the physics of the process is yet to be completely understood.Therefore,in EDM,the complete machining of a workpiece cannot be modelled.Instead individual parts of the process such as estimating the ease of machining or durability of a tool,from a wear viewpoint can be modelled.Much work still needs to be done both in modelling ECM and EDM and this is discussed in the following sections.
Explosion force
Discharge duration
delay time
Reac?on force
Fig.44.Forces acting on wire electrode.
Fig.45.Simulation of electromagnetic ?eld in wire EDM to obtain electromagnetic force acting on wire [161].
Fig.46.Geometrical simulation of wire EDM.
Fig.47.Simulation of temperature distribution along wire electrode (discharge current:90A,discharge frequency:18kHz,wire diameter:0.25mm,workpiece thickness:100mm [90].
S.Hinduja,M.Kunieda /CIRP Annals -Manufacturing Technology 62(2013)775–797
793
6.1.Future work in ECM
To make further progress in modelling the EC milling process, there is a need for geometric pre-processing modules which will shrink/expand the external/internal boundaries of the feature to be machined.Also,a special CAM system is required which can subdivide the tool path into segments,the length of each segment depending on the time step and traverse feed rate.
To save computing time,for each time step,a node should be automatically classi?ed as either active or passive.
Modellers would bene?t from more reliable experimental data especially that for current ef?ciency and polarization voltage. Whilst some data are available for dilute solutions of NaNO3,very little data are available for solutions of H2SO4and HNO3.Another problem with acidic electrolyte solutions is that their electrical conductivity depends on the amount of metal ion concentration.It must be borne in mind that the concentration keeps increasing with usage of these electrolyte solution and very few modellers take this into account.
At the moment,the presence of hydrogen bubbles is accounted for by using empirical equations.The accuracy of an ECM model could be improved by developing more realistic bubbly two-phase models.So far,modelling of two-phase?ow in two-or three-dimensions has not been reported.Modelling of?ow is essential because it has a pronounced effect on the quality of the workpiece, especially if it is a small and deep hole.
6.2.Future work in EDM
There are few EDM simulation programs which have been commercialized.Some machine tools are equipped with empirical simulation tools which can suggest appropriate pulse conditions depending on the workpiece material and working surface area. The machining time can be estimated based on the removal stock [96,169].Considering that it takes many days for machine tool builders to obtain such empirical data each time when a new machine is developed,more effort should be devoted to improving the simulation accuracy.To achieve this,the following problems should be solved.
First of all,it is necessary to obtain correct boundary conditions,especially the energy distribution and arc plasma diameter.Since the simulation of arc plasma is dif?cult,they should be obtained from experiments.Given the correct boundary conditions,powerful multi-physics simulations will realize precise simulation of removal.Since the material removal ef?ciency,the ratio of the ejected to melted volume,can be obtained,the database of pulse conditions could be made off-line. However,since pyrolytic carbon generated on the electrode surface exerts a signi?cant effect on tool wear reduction,as mentioned in Section 4.2.5,quantitative investigation of its protective effect is indispensable.Moreover,to complete the generalized EDM model,integration of the simulation of gap phenomena with the pulse generator control and servo feed control is essential.
https://www.sodocs.net/doc/8f4222461.html,mon future work
For practical use,development of inverse modelling is eagerly anticipated both in ECM and EDM.Preparation of tool electrode is a cornerstone to obtaining high accuracy and productivity.Integra-tion of machine tool controllers and process modelling is another challenging issue.Machine tool builders can design and develop new machine structures and pulse generators using the integrated simulator.The accuracy of adaptive control strategies will be increased by combining both machine and gap simulations. Meanwhile,simulations of hybrid processes between ECM and EDM and synergistic interaction between electrochemical and electrical discharge phenomena in the normal ECM and EDM processes using water-based working?uid are other challenging works.Acknowledgements
The authors would like to express their sincere thanks to J. Atkinson,H.Altena,C.Diver,A.Klink,https://www.sodocs.net/doc/8f4222461.html,uwers,A.Malshe,V. Pattavanitch and M.Zeis for their assistance in the preparation of this paper.We would also like to thank G.Levy,J.P.Kruth,G.Levy, J.A.McGeough,A.De Silva,K.P.Rajurkar,X.D.Yang and other members of the STC-E group for their input.
References
1.Adey RA(1984)in Brebbia CA,(Ed.)Electrostatics Problems,Boundary Element
Techniques in Computer-Aided Engineering,Martinus Nijhoff Publishers,Dor-decht.
2.Alacron PF(1978)Boundary Elements In Potential And Elasticity Theory.
Computers and Structures10:351–356.
3.Alkire R,Bergh T,Sani RL(1978)Predicting Electrode Shape Change with Use
of Finite Element Methods.Journal of the Electrochemical Society125(12): 1981–1988.
4.Altena H(2000)Precision ECM by Process Characteristic Modelling,Glasgow
Caledonian University.(PhD Thesis).
5.Araie I,Sano S,Kunieda M(2008)In?uence of Electrode Surface Pro?le on
Discharge Delay Time in Electrical Discharge Machining.IJEM13:21–28.
6.Bhattacharyya S,Ghosh A,Mallik AK(1997)Cathode Shape Prediction in
Electrochemical Machining Using a Simulated Cut-and-Try Procedure.Journal of Materials Processing66:146–152.
7.Bleys P,Kruth JP,Lauwers B,Zryd A,Delpretti R,Tricarico C(2002)Real-
Time Tool Wear Compensation in Milling EDM.Annals of the CIRP51(1): 157–160.
8.Bommeli B,Frei C,Ratajski A(1979)On the In?uence of Mechanical Perturba-
tion on the Breakdown of a Liquid Dielectric.Journal of Electrostatics7: 123–144.
9.Bortels L,Deconinck B,Van Den Bossche B(1996)The Multi-dimensional
Upwinding Method as a New Simulation Tool for the Analysis of Multi-ion Electrolytes Controlled by Diffusion,Convection and Migration.Part I:Steady State Analysis of a Parallel Plane Flow Channel.Journal of Electroanalytical Chemistry404:15–26.
10.Brebbia CA(1978)The Boundary Element Method for Engineers,Pentech Press,
London.
11.Brookes DS(1984)Computer-Aided Shape Prediction of Electro-chemically
Machined Work Using the Method of Finite Elements,The University of Man-chester.(PhD Thesis).
12.Carslaw HS,Jaeger JC(1959)Conduction of Heat in Solids,Oxford Science
Publications,New York,USA.
13.Cetin S,Okada A,Uno Y(2004)Effect of Debris Accumulation on Machining
Speed in EDM.International Journal of Electrical Machining9:9–14.
14.Chang CS,Hourng LW(2001)Two-dimensional Two-phase Numerical Model
for Tool Design in Electrochemical machining.Journal of Applied Electrochem-istry31:45–154.
15.Christiansen S,Rasmussen J(1976)Numerical Solutions for Two-dimensional
Annular Electrochemical Machining Problems.Journal of the Institute of Mathe-matics and its Applications18:295–307.
16.Clark WG,McGeough JA(1977)Temperature distribution along the gap in
electrochemical machining.Journal of Applied Electrochemistry7:277–286. 17.Collett DE,Hewson-Browne RC,Windle DW(1970)A Complex Variable
Approach to Electrochemical Machining Problems.Journal of Engineering Mathematics4:29–37.
18.Crookall JR(1979)A Theory of Planar Electrode-face Wear in EDM.Annals of
the CIRP28(1):125–129.
19.Dabrowski L,Kozak J(1981)Electrochemical Machining of Complex Shape
Parts and Holes(Mathematical Modelling and Practical Experiments).Proc.
22nd Int.MTDR Conference,Manchester.
20.Danson DJ,Warne MA(1983)Current Density Voltage Calculations using
Boundary Element Techniques.The Int.Corrosion Forum,NACE,Anaheim, California.
21.Das S,Lotz M,Klocke F(2003)EDM Simulation:Finite Element-based Calcula-
tion of Deformation,Microstructure and Residual Stresses.Journal of Materials Processing Technology142:434–451.
22.Das S,Mitra AK(1992)Use of Boundary Element Method for the Determina-
tion of Tool Shape in Electrochemical Machining.International Journal for Numerical Methods in Engineering35:1045–1054.
23.Datta M,Landolt D(1982)High Rate Transpassive Dissolution of Nickel with
Pulsating Current.Electrochemica Acta27/3:385–390.
24.Dauw D(1986)On the Derivation and Application of a Real-Time Tool Wear
Sensor in ED.Annals of the CIRP35(1):111–116.
25.Davey K,Bounds S(1997)A Preconditioning Strategy for the Solution of Linear
Boundary Element Systems using the GMRES Methods.Applied Numerical Mathematics23(4):443–456.
26.DeBarr AE,Oliver DA(1968)Electro-chemical Machining,Macdonald&Co.Ltd,
London.
27.Deconinck J(1992)Current Distributions and Electrode Shape Changes in Elec-
trochemical Systems,Springer-Verlag,Berlin/Heidelberg/New York.
28.Deconinck J(1994)Mathematical Modelling of Electrode Growth.Journal of
Applied Electrochemistry24:212–218.
29.Deconinck D,et al(2011)Study of the Effects of Heat Removal on the Copying
Accuracy of the Electrochemical Machining Process.Electrochimica Acta 56:5642–5649.
S.Hinduja,M.Kunieda/CIRP Annals-Manufacturing Technology62(2013)775–797 794