ABSTRACT: This paper proposes a general applicable numerical algorithm for the simulation of two dimensional electrode shape changes during electrochemical machining processes. The computational model consists of two coupled problems: electrode shape change rate computation and a moving boundary problem. The innovative aspect is that the workpiece shape is computed over a number of predefined time steps by convection of its surface with the velocity proportional with and in the direction of local shape change rate. An example related to the electrochemical machining of a slot in a stainless steel plate in order to demonstrate the strong features of the applied method will be presented.
Electrochemical machining (ECM) has become over the past decades one of the important methods in non-conventional tooling processes. Because ECM is a non-mechanical machining application, regardless the mechanical properties, the metals can be tolled with high removal rates. Common ECM processes include: drilling, duplicating and sinking operations in the manufacturing of dies, creation of cavity holes and slots, compressor and turbine blades etc [1, 2].
The complex behavior of electrochemical processes subjects ECM to several drawbacks and one of the most critical problems is achieving a high degree of accuracy. Therefore the prediction of the final shape of a workpiece, obtained after ECM, is an important issue. Considerable costs and time are saved when reliable predictions are available. Simulations enable to design proper electrode tools and synthesize adequate electrolytes, hence avoiding a lot of "trial and error" cycles.
Different models can be used to describe the electrochemical phenomena that occur in the electrolyte and at the electrode interface during the ECM process. On the other hand one can apply different solution techniques and solution algorithms to solve the governing equations.
The first electrode shape change simulations are based on analytical techniques. McGeough  solved the ECM tool design using the "sin rule". This method can be used to obtain an approximate shape but it fails when the tools have sharp discontinuities caused by the neglecting of stray current effects and the assumption of parallel flux lines.
More recently in the field of analytical solutions for tool design problems Adler and Clifton  used an analytical direct computation of 2D ECM, representing the workpiece by Fourier's series. In addition, an experimental validation of the model is given.
Due to the ever increasing performance of digital computers and numerical methods the calculation of ECM processes highly has become more and more popular. The first numerical simulations of ECM where based on current density distributions obtained using the Potential Model (PM) describing conservation of charge. Several authors applied the PM to compute the electrode change rate and compared the obtained results with experiments. Pretice and Tobias  used the Finite Difference Method (FDM) to solve a changing electrode profile with primary, secondary and tertiary current density distribution. They present a complete study of electrode shape change as a function of Wagner's number for different angles of a notch electrode.
Jain et. al. [5, 6] demonstrated the application of the Finite Element Method (FEM) for the 2D computation of the metal removal rate. Their model incorporated the effects of simultaneous changes in the electrolyte flow velocity and temperature rise.
Deconinck and Van den Bossche [7, 8] applied the BEM to solve Laplace equation, with non-linear boundary conditions to account for the electrode charge transfer reaction at the electrode. They applied the method to model the deburing of a drilled hole.
A study of the workpiece shape prediction during the ECM process is given in . The authors used the Boundary Element Method (BEM) with linear and quadratic elements in order to predict the workpiece shapes under simplified boundary conditions. The method resulted in small errors between the computed and exact solution.
Recently more complicated models have been developed in order to deal with the complex situations that occur during ECM. Chang  investigated a two dimensional, two phase (liquid/ gas) model in order to predict tool shape for changes.
De Silva and Altena [11, 12] developed models in order to increase the precision in ECM. They include the polarization voltage and the influence of electrolyte concentration and conductivity within small gaps in the current density and efficiency models. In order to achieve better precision they used a pulsed voltage to drive the ECM process.
Only a very limited number of publications deal with three dimensional (3D) current density distribution simulations in electrochemical reactors. Applications that consider current density distributions for more complex 3D geometries can be found in the field of cathodic protection [14, 15]. In the field of 3D ECM Kozak [16, 17] presented the concept and prototype development of a computer aided engineering system for the solution of ECM manufacturing problems such as electrode tool design, selection of machining and parameters optimization. Kozak used the FDM in order to solve the mathematical model but the numerical simulations do not deal with real 3D electrode shape change. 3D ECM modeling reaches maturity once with the papers of Purcar and Bortels [18-20] where the principles of advanced CAD integrated approach for 3D electrochemical machining simulations are established.
In all mentioned cases the authors make use of an Euler scheme for the integration of Faraday's law with respect to time, combined with the string or marker theory . The method derives from the Lagrangean approach of changing boundary problems that supposes the mesh of the boundary is attached to the moving front. The method is not well suited to compute bifurcation fronts (swallowtail see figure 1 left). In addition the method shows stability and oscillation problems, high deformation of the boundary elements and inability in naturally dealing topology changes.
(left) string method flaws swallowtail; (right) LSM results under the same boundary conditions; ( - - - ) initial shape of the electrode, ( __ ) final shape of the electrode.
In order to avoid the above mentioned problems (see figure 1 right) this paper proposes a general and accurate electrode growth algorithm based on the Level Set Method (LSM). The method was introduced by Osher and further developed by Sethian [9, 10] based on a comparative study for the numerical modeling of the interface evolution of flame propagation. It has the strong feature that the boundary can freely evolve in such a way that, for example, it can break and naturally merge together. This technique has become very popular and has been applied in a wide range of problems such as fluid mechanics, combustion, manufacturing of computer chips, image processing, structural mechanics, electronic devices optimization [21-24], and more recently in electrochemical modeling [25, 26]. Purcar and all adapted LSM in  to accurately model the electrode growth for general geometries and the results have been validated with experimental data.
This paper proposes a general applicable numerical algorithm for the simulation of two dimensional ECM processes. The computational model consists of two coupled problems:
Electrode shape change rate computation based on the PM and
Moving boundary problem based on the LSM.
Basic principles of the 2D electrode shape change algorithm based on LSM
The continuous process time is divided into a sequence of discreet time steps. At each time step two problems need to be solved: the electrochemical model, which eventually provides the local growing speed of the electrode(s), and the electrode shape change problem, which updates the position of the front. The coupling of the two problems is assured via Faraday's law.
The basic principles of the ECM simulations based on LSM that have been developed in this work are illustrated in figure 2. According to this the total process is divided in two time loops. The first one, being the outer time loop, or macroscale time loop (TF), is related to the calculation of the electrode shape change rate, vLSM, based on the PM and Faraday's law at each time step Î”tF, until the total plating time TF is reached. The inner time loop or microscale time loop deals with the time integration of equations of the Level Set equation, with time steps Î”tLSM within the interval Î”tF in order to find the new shape and position of the electrode.
Structure of the electrode shape change algorithm
Electrochemical process model (Potential Model)
Electrochemical machining is the anodic dissolution process that occurs when an imposed electrical current flows between the workpiece (anode) and the worktool (cathode). The electrolyte usually flows with high velocity in the gap between the anode and cathode in order to intensify the mass/charge transfer with the anode's double layer, and to remove the dissolution products gas bubbles and heat generated in the gap. The electrolyte is mostly an aqueous solution of neutral salts like NaCl, NaNO3 and NaClO3. When a potential difference is applied across the electrodes, several reactions can occur at the anode and cathode.
Although flow, thermal and gas evolution effects might be relevant in ECM, as a first approach it is assumed that the electrolyte is sufficiently well stirred or refreshed such that these effects are neglected and only charge transport with ohmic resistivity effects is considered. Hence the PM  described by Laplace's equation (1), for an electrolyte system having a potential drop U in [V] and electric conductivity ï³ in [S/m], holds:
Ohm's law holds but the conductivity ï³ does not need to be constant. Indeed it is possible to consider a local varying conductivity (e.g. function of the temperature T, ï³(T) or gas).
No current flows through insulating walls and air and therefore the normal current density at each point on these surfaces is zero:
where the subscript n refers to the normal direction.
At the anode, in principle the metal Me dissolves:
The polarisation behaviour of this reaction is quite often accurately described using a linear relation with the efficiency depending on the current density. The polarization can be written as the difference between the (imposed) metal potential V and the potential U in the solution adjacent to the electrode and is function of the amplitude of the current density Jn, normal to the electrode surface:
In which A, B are polarisation constants and E0 is the equilibrium potential. In order to obtain the electrochemical reaction (3) at a given net rate, in a given direction, a driving force is required. This driving force is called the overpotential Î· and is written as:
At the cathode, hydrogen gas and hydroxyl ions are produced:
and the polarization behaviour can also be linearised as in equation (3).
However in the most general situation, the current density distribution can be described as a function of the metal potential V, the electrolyte potential adjacent to the electrode U, and the equilibrium potential E0:
This function is often obtained by measurements and approached using a spline interpolation function .
The PM cannot make a distinction between parallel electrode reactions. In case one needs to include this effect, it is necessary to have as addition information the efficiency Î¸ of the metal removal as a function of the current density.
The ECM rate h at each point of the electrode is proportional to the local normal current density Jn. The direction of change is normal to the electrode surface and the ECM rate equation can be written in vectorial form as:
where is the unit normal oriented outwards to the computational domain, M holds for the atomic weight of the metal in [kg/Mol], z is the number of electrons exchanged in the metal deposition reaction, F is Faraday's constant in [C mol-1], dt the time interval within the local rate of mass is produced in [s], and Î¸ the efficiency of the dissolution. As Jn is positive at the anode, dissolution occurs, while at the cathode, Jn is negative and the electrode grows inward the solution. Remark that equation (8) also holds for insulating boundaries, because there Jn is zero.
When one of the electrodes is moving, a velocity term is to be added to equation (8) yielding:
where is the electrode speed or electrode feed rate. In many practical cases is the constant feed rate of the cathode tool used in the ECM process.
Electrode shape change model (Level Set Method)
The explicit representation of the interface that is used in the classical FEM forbids deep boundary or topological changes, such as creation of holes. This limitation is the main reason of the low performance generally associated to the electrode shape change problems. However, the LSM developed in , which consists in representing the interface with an implicit function, overcomes these kinds of deep changes. The LSM is a numerical technique first developed for tracking moving interfaces. It is based on the idea of implicitly representing the interface as a LS curve of a higher dimension function . The moving interface Î“, considered between the two sub-domains Î©1 and Î©2 with different material properties is conventionally represented by a zero LS function (see figure 2), and the sign convention is given by:
There are many possibilities to express but in most of the LS applications it is described as:
where "d" is the so-called signed distance (SD) function, the distance from the interface to every point of the analyzed domain Î© = Î©1Î©2. This ensures that the material interface Î“ is represented by the zero LS function.
Sign convention and distance function.
The LS values at the nodes show whether an element is cut or not by the material interface. The following criterion is used in order to establish which of the FE nodes are fully enriched:
For plating applications the normal current density distributions along the workpiece determine the electrode growth rate as described by Faraday's law. The LSM uses the local growth rate (9) as the vector speed  in the following convection equation:
The Physical meaning of this equation is that the electrode boundaries are convected proportional to the local growth rate.
Numerical solution method
Equations (1) and (13) generally describe geometries with arbitrary boundary conditions which cannot be solved exactly. The solutions can be obtained using different discretization techniques, which approximate the continuous solution as a set of quantities at discrete locations, the so-called grid or mesh. In this work the finite element method (FEM) is used to both the electrochemical process model and the electrode shape change model.
To reduce the continuum problem described by equation (1) to a discrete system of algebraic equations, the standard FEM is used. The definition domain Î© of equation (1) is discretized into elements (triangles). The FEM formulation for the Laplace equation (1) on the domain can be written as :
Ni and Nk are respectively the weight and shape functions, i, k =1â€¦3 (the vertices of the triangles), Î“el is the boundary of the element and Î©el that part of the domain Î© (in 2D a surface) occupied by an element. Equation (11) can be rearranged for all nodes k of Î© in the following system of equations:
with the system matrix, the vector of unknown potentials in the nodes and the first term of equation (11). The introduction of proper boundary conditions makes the system of equations unique.
When equation (7) is used as boundary condition, equation (15) can be directly solved. Otherwise the righthandside of equation (15) becomes nonlinear at the electrode boundaries and in that case an iterative, Newton-Raphson, process is used to obtain the solution of equation (15).
Level Set Method
Signed Distance Computation and Velocity Extension
Common choice for the SD function is the solution of the Eikonal equation with boundary condition Î¦ = 0 (the zero LS function which coincides with Î“(0)). This method is in literature also referred to as the Fast Marching Method (FMM) . The FMM first initializes the points that are adjacent to the interface with the minimum distance between these points and the interface. Using these adjacent points as boundary condition for the Eikonal equation, the final solution is found by upwind differentiation on the whole computational domain Î©. The above mentioned procedure is used in our implementation.
Since the LSM only needs to be effective for the Î¦ = 0 isocontour, equation (13) is not necessary to be integrate on the entire computational domain. The integration algorithm perfectly fits if it is applied on a narrow band around the interface Î“(0). The computation of the narrow band is performed during the SD computation step. Using FMM, all triangles in certain distance from the moving material interface are marked. The narrow band width is usually computed as the product between the maximum velocity vmax at the moving front and the total LS process time TLSM. Here the narrow band is chosen in such way to not be larger than 5 layers of triangles.
Velocity distribution from the moving interface Î“ to the LS computation region is called the "velocity extension". Malladi et. al. , suggested extrapolation of the velocity from the front on a Cartesian mesh. At each grid point, the value of the velocity at the closest point on the interface is simply used at that grid point. This method works fine for Cartesian grids but for triangular meshes proves to be very time consuming. Sethian used the FMM for the velocity extension . Hence, by constructing the SD function Î¦ = Â± d at each grid point using the FMM, the velocity is simultaneously updated according to equation. The method applied in this work is similar to the one in Sethian, and uses the narrow band extrapolation.
In order to find the new position of the material interface, equation (13) is numerically integrated in space and time with the initial conditions (11) and with a velocity defined in every grid point of the narrow band (the total number of nodes in the narrow band is Nnarrow). The majority of the LSM applications are based on Cartesian grids . In our work, a triangular unstructured mesh, the same as for the XFEM analysis is used for the numerical integration of the LS equation.
In the first step, equation (13) is discretized by using the residual distribution formulation with the standard Galerkin FE shape functions and the multidimensional upwind method . For the residuals distribution to the triangles vertices, the Low Diffusion A (LDA) scheme has been used . The time integration is approached using a second order accurate Petrov-Galerkin (PG) formulation. The set of equations derived from (14) may be written as:
where is the time depending vector of Î¦ in each node [Nnarrow Ã- 1] ; is the vector of the unknown Î¦ at each time step [Nnarrow Ã- 1]; is the global mass matrix [Nnarrow Ã- Nnarrow] and is the global convection matrix [Nnarrow Ã- Nnarrow].
One of the LS demands for such time domain integration scheme is that the scheme must follow some underlying properties as consistency and monotonicity, which are substantial in the integration of the LS equation. These criteria are fulfilled by combining LDA scheme for space integration and PG formulation with mass matrix stabilisation for the time integration .
Time Step Definition
The integration of equation (13) requires that the total process time is divided into a number of time steps . The time step definition is further related to the re-initialisation of the LS equation to the SD function. The concept of re-initialisation is a powerful numerical tool in the LS theory and it has been extensively discussed in .
The time step is chosen such that the front crosses no more than one grid cell of size h each time step:
max||v|| Î”tâ‰¤ h (17)
where the maximum is chosen from all values of v at all possible mesh points in the narrow band.
Numerical Example - Slot electrochemical machining in a stainless steel plate
The 2D electrode shape change problem is created using a cross-section in a plane perpendicular to the middle of the electrodes as presented in figure 4. The workpiece is a stainless steel plate of 30411 mm3. The worktool is 819.90 mm3 and has the bottom part parallel with the workpiece, hence yielding an inter-electrode gap of 100 Î¼m. The electrolyte domain is defined as a box with insulating walls of dimensions 514031.1 mm3.
2D computational domain (dimensions in mm).
The studied configurations.
Functions of the worktool characteristics, four different situations are analysed (figure 5):
Static worktool entirely active (SWA);
Static worktool partially active (SWP); only the bottom part, parallel to the workpiece is active;
Moving worktool entirely active (MWA);
Moving worktool partially active (MWP); only the bottom part, parallel to workpiece is active.
Physico-chemical input parameters
All the necessary process parameters have been extracted from literature [Error: Reference source not found-20]. The physico-chemical data for the ECM process are characteristic for a NaNO3 bath at 27 °C and contain efficiency and polarization data correlations for a variety of NaNO3 concentrations and pulse times. The values corresponding to a NaNO3 concentration of 70 g/l for a time pulse of 1 ms were used in the simulations. The electrical conductivity of the bath at 27 °C is 8.0 S/m.
The anodic polarisation behaviour is described using a linear current density-overpotential relation with fitted parameters as given below (Ja in A/m2):
Ja=2.91 105 (V-U)-9.16 105. (18)
The efficiency of the electrolyte has been retained from  and is obtained from parallel flow cell measurements. This yields a numerical data set as represented in figure 5, directly insertable in the electrochemical simulation model via equation 7.
Polarisation curve and efficiency for the anodic reaction: current density, efficiency.
The efficiency of the metal dissolution (in comparison to the oxygen evolution side reaction), ranges from almost zero, for small current densities, to about 80 %, for current densities above 1.5 106 A/m2. The cathodic polarization behaviour (mainly hydrogen evolution) is assumed to be nearly primary:
Jc=1.0 106 (V-U)-1.0 106. (19)
Since only hydrogen gas is evolved at the cathode, the shape of the electrode remains unaltered during electrolysis and therefore will not be considered within the next section.
Electrode shape change simulations
Using the (FEM), the governing 2D Laplace equation with non-linear boundary conditions that characterizes the electrochemical model is solved over the whole computational domain described in figure 4 at every time step in order to find the normal current density distribution Jn at the workpiece surface. Then, the electrode shape change rate is computed on the common boundary between the electrolyte and the workpiece from equation (8) and the efficiency curves from figure 6. The electrode shape change velocity is further extended at every node of the mesh of the workpiece zone.
In order to capture the ECM profile with high accuracy, the mesh is composed of a very fine unstructured mesh on the workpiece zone (c.a. 10 000 elements).
The microscale time step âˆ†tLSM used for the integration of the level set equation (13) is computed according to equation (17) and is function of the maximum electrode shape change velocity, normal at the electrode surface at the beginning of the simulations, and dimension of the average size of the discretization element at the electrode (c.a. h=1.0e-4 [m]). An overview of the calculated time steps for each of the cases SWA, SWP, MWA and MWP is given in table 1.
Time steps for the integration of the LS equation.
Static worktool entirely/partially active
For a static worktool, as is the case for SWA and SWP, the inter-electrode gap will increase indefinitely as machining proceeds. Both the SWA and SWP configurations need a very high process time, respectively 389.24 and 583.85 s, in order to perforate the plate.
Profile evolution of the ECM slot for the set-ups SWA and SWP.
Figure 7 shows the ECM profile evolution for the SWA and SWP set-ups. It can be seen that the SWA produces a perforation than is four times larger than the worktool width, while that generated by the SWP is only 3 times as wide.
The computations are continued up to the final time of TF = 875.78 s, which yields 2 additional time steps for the SWA and 1 for the SWP. In both cases a rounding of the edges of the perforation can be noticed.
Moving worktool entirely/partially active
However, increasing gap sizes are usually not desirable in ECM, particularly not when an accurate shape must be produced on the workpiece . To avoid the gap size widening a commonly used technique is to work with moving cathode tools. An additional benefit is that in this case the total process time is much shorter when compared to static ECM.
Front and side gap generated at time tF = 19 s in the slot ECM with non-insulated MWA (left) and insulated MWP (right) electrode tool; FG front gap, SGB side gap at the beginning of the tool electrode and SGE side gap at the end of the tool electrode.
The growing of the front gap (FG), as presented in figure, is compensated by the feeding rate of the electrode tool. However, since there is no lateral feed, the growth of the side gap SGE (SGE is defined as the distance along the horizontal coordinate between the cathode tool and the upper tangent point between the straight and curved portion of the working boundary) at the end of the process cannot be compensated (figure 8 left). This leads to wedge contours and overcuts of the working boundary, as shown in figure 8 left and 9 MWA.
Profile evolution of the ECM slot for the MWA and MWP
Using an insulated electrode tool increases the accuracy of the final side gap SGE, which becomes equal with the side gap at the beginning of the electrode tool SGB, as illustrated in figure 7 right and Error: Reference source not found MWP. SGB is defined as the distance along the horizontal coordinate between the corner of the worktool and the working boundary.
Detail of the ECM profile evolution at the perforation zone for the MWA (left) and MWP (right) setups.
For both MWA and MWP configurations, the perforation occurs at tF =20 s. Two more time steps are carried on after the perforation. At the end of the process time the bottom part of the cathode tool is at the same level (0.020) of the bottom face of the plate. A side effect of the perforation is the emergence a "neck". In practice, "necks" are usually unwanted side effects of ECM and a lot of effort is spent in designing tools to avoid them. As for the MWA configuration the lateral wall of the worktool is active, it plays an important role in the rounding of the "neck", see figure 10 left.
On the other hand, in the MWP case, the lateral wall of the worktool is insulated and the active part of the worktool (bottom part) only is not able to round the "neck". The size of the neck is further reduced but it still remains as a protrusion at the exit of the slot, see figure 10 right.
Calculation times for the 2D model, based on the Laplace equation with non-linear boundary conditions and LSM for the electrode shape change, are ca. 10 minutes for each of the static worktool configurations and ca. 15 minutes for each of the moving worktool configurations. The calculations have been carried on a PC Pentium 4, 2.4 GHz, 2 GB RAM.
A new method for the mathematical description of electrode shape change during ECM processes has been presented. The method has the strong feature that it can treat any electrode shape where the speed is governed by the electrode shape change rate. The capacity of the 2D electrode shape change algorithm in dealing with topology changes such as the perforation of a plate has been emphasised. The 2D electrode shape change algorithm based on the LSM naturally handles topology changes and allows continuing the simulations after perforation in order to analyse for example the "neck" smoothing.
Simulations have been carried out for a succession of 4 type of worktools e.g. static worktool entirely active, static worktool partially active, moving worktool entirely active and moving worktool partially active.
To conclude one can say that both the numerical simulations prove that the developed software can easily be used to analyse and improve an ECM process in terms of:
Minimisation of the process variation and optimisation of the time machining quality;
Predictability of the process, for example, by determining rapidly the effects of variation in a certain process parameters;
Calculation/simulation of tool-electrode shape, thus avoiding the costly trial and error tool design for any new applications.
This work was supported within the research program POSDRU/89/1.5/S/57083, project title "Development and Evaluation of an Advanced Mathematical and Numerical Analysis Techniques for Modeling Electrode Shape Changes in Electrochemical Process.