# Iterative surface warping to shape craters in micro-EDM simulation

Anthony Surleraux, Jean-Philippe Pernot, Ahmed Elkaseer, Samuel Bigot

► To cite this version:

Anthony Surleraux, Jean-Philippe Pernot, Ahmed Elkaseer, Samuel Bigot. Iterative surface warping to shape craters in micro-EDM simulation. Engineering with Computers, 2016, 32 (3), pp.517-531. 10.1007/s00366-016-0439-0 . hal-05042455

HAL Id: hal-05042455

https://hal.science/hal-05042455v1

Submitted on 22 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Iterative surface warping to shape craters in micro-EDM simulation

Anthony Surleraux<fnref n="1" /> · Jean-Philippe Pernot<fnref n="2" /> · Ahmed Elkaseer<fnref n="3" /> · Samuel Bigot<fnref n="1" />

**Abstract** This paper introduces a new method for simulating the micro-EDM process in order to predict both the tool's wear and the workpiece's roughness. The tool and workpiece are defined by NURBS patches whose shapes result from an iterative crater-by-crater deformation technique driven by physical parameters. Through hundreds of thousands of local surface warping, the method is able to compute the global as well as the local shapes of the tool and workpiece. At each step, the warping vector and function are computed so as to be able to generate a spherical crater whose volume is also controlled. While acting very locally to simulate the real physical phenomenon, not only the method can evaluate the tool's wear from the overall final shape at a low resolution level, but it can also estimate the workpiece's roughness from the high resolution level. The simulation method is validated through a comparison with experimental data. Different simulations are presented with an increase in computation accuracy in order to study its influence on the results and their deviation from expected values.

**Keywords** Micro electrical discharge machining · A priori wear estimation · Surface roughness · Local deformation · NURBS surface warping · Multi-resolution

<figure><img src="image_2.png" /></figure>

## 1 Introduction

Electrical discharge machining (EDM) is a manufacturing process that consists in removing parts of a material with electrical discharges and is characterized by its ability to machine any conductive material regardless of its hardness [1].

Although various forms of EDM exist, all of them share the same concept as described in Fig. 1. Two electrodes (the tool and the workpiece) are separated by a dielectric fluid. Both electrodes are submitted to an electrical current and as the gap between the electrodes diminishes, the intensity between them increases until it reaches what is called the dielectric breakdown voltage. At this point, the dielectric cannot act as an insulator anymore and allows current to flow from one electrode to another leading to the apparition of a plasma channel. The plasma's temperature ranges from 8000 to 12,000 °C and in some cases can reach up to 20,000 °C [2]. This leads to evaporation and melting of both the tool and the workpiece. When the current is stopped, the dielectric fluid rushes back where the plasma stood and evacuates resulting debris.

Micro-EDM (μEDM) shares the same underlying concepts of EDM. It simply tackles with dimensions in the order of the micron. However, while the process setup is modified in order to reduce its effect, tool wear often becomes the main factor of imprecision. Figure 2 describes the tool wear phenomenon. As material is removed from the workpiece, craters also appear on the tool, changing its shape. This is what is referred to as tool wear. As the process goes on, change in the shape of the tool can be drastic and, in the scope of μEDM, has major repercussions on the final profile of the workpiece. While this phenomenon also occurs in classical EDM, its influence on the final shape is significantly smaller. When using conventional machining

<footnote>Samuel Bigot
bigots@cardiff.ac.uk</footnote>

<footnote>Jean-Philippe Pernot
jean-philippe.pernot@ensam.eu</footnote>

<footnote n="1">¹ Manufacturing Engineering Centre, Cardiff University,
Cardiff, UK</footnote>

<footnote n="2">² Arts et Métiers ParisTech, LSIS UMR CNRS 7296,
Aix-en-Provence, France</footnote>

<footnote n="3">³ Faculty of Engineering, Port Said University, Port Said,
Egypt</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Principle of EDM</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Tool wear influence in EDM. a Tool and workpiece before machining. b Tool and workpiece after machining. The black workpiece is the desired final profile while the orange one is the actual profile that is obtained</figcaption></figure>

strategies, the electrodes' shapes quickly deviates from the original ones. Thus, as of now, μEDM milling is the preferred machining strategy for the fact that proven methods exist to mitigate the influence of tool wear on the final result [3], while for similar applications die-sinking μEDM may require a dozen or more tools before obtaining the desired geometrical tolerances. Being able to predict the tool wear more accurately would enable us to design more efficient machining strategies [4], in particular for die-sinking EDM where using extra volumes on specific parts of a tool electrode could compensate partially for the wear and drastically reduce the number of electrodes required. To achieve this, a better theoretical understanding of the wear phenomena is required so as to be able to simulate a priori, i.e., before machining, the expected final shape as well as the surface's roughness.

Until now, most modeling works that have been conducted focus on the development of single crater (or single spark) models in which physical equations are used under certain assumptions to determine the shape of the resulting craters of a single spark. A few have proposed simulation methods of μEDM usually using a grid of points to model

the involved geometries [5, 6]. Those existing methods result in inaccurate surfaces with discontinuity predicted as steps and sharp corners, which might be the cause of significant cumulative errors in surface approximation when simulating a long sequence of discharges. A more detailed analysis is proposed in Sect. 2.1.

Against this background, our aim is to develop a new approach for the accurate modeling of the local deformations altering the complex 3D shaped micro-electrodes. In order to do so, and instead of using a cloud of points to model the geometries, the use of parametric surfaces is made in order to preserve continuity. The idea is to mimic the iterative crater-by-crater appearing phenomenon through an iterative surface warping approach. At each step, the warping vector and function are computed so as to be able to generate a spherical crater whose volume is also controlled. Thus, through hundreds of thousands of local surface warping, the method is able to compute the global as well as the local shapes of the tool and workpiece. More precisely, the proposed approach has the following advantages:

* it uses NURBS to define the geometries to be able to model complex shapes at different resolutions: shape of the craters locally, and shape of the wear globally;

* it mimics the physical μEDM process while performing local surface warping to generate micro-craters;

* it is driven by physical parameters so as to generate an accurate approximation of the wear and roughness;

* it can be inserted in a shape optimization loop so as to find the optimal shape of the initial tool that would produce a targeted workpiece.

Considering the first advantage, it is clear that the proposed approach is completely different from the existing surface deformation techniques found in the literature (see Sect. 2.2). Through a huge amount of local deformations, the method can evaluate the overall tool wear as well as the workpiece's roughness. This underlines the multi-resolution character of the proposed approach even if multi-resolution surfaces are not used. Considering the last above mentioned advantage, our final aim is to be able to set up a complete shape optimization framework that would help finding the shape of the initial tool to be used to obtain a targeted workpiece after machining. Of course, the proposed crater-by-crater simulation method is part of this framework. It is called several times during the optimization.

The paper is decomposed as follows. First, Sect. 2 reviews the state-of-the-art in both μEDM modeling and surface deformation techniques. Then, the crater-by-crater deformation technique is presented and detailed in Sect. 3. In Sect. 4, the proposed approach is then validated through several experimentations where numerically simulated
magnitudes are compared to experimentally measured ones. The final section concludes this paper and sketches future perspectives.

## 2 Literature review

Since our approach is based on an iterative surface warping driven by physical parameters derived from the real μEDM process, this section will review the state-of-the-art in both μEDM modeling and surface deformation.

### 2.1 μEDM modeling

#### 2.1.1 Introduction

Electrical discharge machining can be considered as a heat transmission problem for which the temperature distribution in the workpiece shall be sought in order to determine the shape of the crater. Three main components are to be studied: energy distribution, size and shape of the discharge channel and material ejection. A portion of the energy is evacuated in the dielectric and by radiation, another portion is evacuated in both electrodes by conduction and a negligible part is lost by convection in the workpiece. The exact repartition of the thermal energy relies on different factors that include the boiling temperature of the electrodes.

Finally material evacuation has to be considered. Unfortunately all of the molten metal is not evacuated by the dielectric fluid and considering so will lead to a theoretical material removal rate (MRR) that will be greater than experimental data. This introduces the concept of plasma flushing efficiency which is defined by the fraction of the molten metal that is actually removed. Work by Descoeudre et al. [7] has established that the plasma flushing efficiency is dependent on the thermal expansion coefficient of the electrode, the amount of molten material, the plasma channel radius, the thermal properties of material and the flushing conditions.

#### 2.1.2 Thermo-electrical modelling

These EDM erosion models are based on the research of solutions for the heat conduction problem described by the following differential equation:

$$ \frac{\partial T}{\partial t} - \alpha \Delta T = 0 \qquad (1) $$

where T is the temperature (K) and Δ is the Laplace operator.

$$ \text{with } \alpha = \frac{K_t}{\rho C_p} \qquad (2) $$

where K_t is the thermal conductivity of the material (in W m⁻¹ K⁻¹), ρ is the density (kg m⁻³) and C_p is the specific heat capacity (J kg⁻¹ K⁻¹). Some models take the melting heat m into account in the formulation of the thermal diffusivity:

$$ \alpha' = \frac{K_t}{\rho(C_p + \frac{m}{T_m})} \qquad (3) $$

where T_m is the melting temperature of the material.

The boundary conditions relate to the geometry being considered. Several approaches of the problem exist. Some have used a semi-infinite cylinder, a disk heat input and assumed that the energy was equally shared among the cathode and anode [8, 9]. This was followed by a model using a two-dimensional heat flow model bounded by an adiabatic cylinder [10]. Dibitonto et al. [11] used different sources for each electrode: a disk heat source for the anode and a point heat source for the cathode erosion. Another element of the boundary conditions is the plasma heat flux that can be defined as the following in the case of an uniform temperature distribution:

$$ q = \begin{cases} F_c \frac{UI}{\pi r_c^2} & \text{disk source} \\ q_0 & \text{point source} \end{cases} \qquad (4) $$

where F_c is the fraction of energy being transferred, U is the gap voltage applied to the electrodes, I is the current, q_0 is a constant and r_c is the plasma channel radius.

It has been observed that there is variation in the plasma radius. There is an expansion followed by stabilization. A common function describing the plasma radius, R\_plasma(t) as a function of time is given as follows:

$$ R_{\text{plasma}}(t) = K t^n \qquad (5) $$

where K and n are coefficients depending on the machining parameters and material properties. Those values are determined empirically.

Concerning the plasma heat flux, researchers have established that the temperature distribution in the plasma channel is not uniform and have adopted a Gaussian distribution as following [7]:

$$ q(r, t) = \frac{F_c UI_{\text{ON}}}{2\pi R_{\text{plasma}}(t)} e^{-4.5r^2/R_{\text{plasma}}^2(t)} \qquad (6) $$

where t\_ON is the duration of the spark.

In order to provide a better understanding of the various elements that have been described, a model from Beck is presented here in more details [12]. As visible on the parameterisation described on Fig. 3, the heat source used is that of a disk. The geometry is a semi-infinite cylinder insulated everywhere except for the region in contact with the heat flux. Solving Eq. 1 yields the following solution:
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Parameterisation of beck's model</figcaption></figure>

$$T(r, z, t) = T_i + \frac{2q \cdot r_c}{K_t} \left( \frac{r_c B(z, t)}{r_0^2} + \sum_{i=1}^{\infty} \frac{C_i(z, t) J_0(\lambda_i r) J_1(\lambda_i r_c)}{2[\lambda_i r_c J_0(\lambda_i r_c)]^2} \right) \quad (7)$$

where $T_i$ is the initial temperature, $J_0$ and $J_1$ are the Besel functions of the first kind of the first and second order respectively and $\lambda_n$ are the roots of $J_1(\lambda_n r_0)$ and:

$$B(z, t) = \sqrt{\alpha t} \cdot \operatorname{ierfc}\left(\frac{z}{2\sqrt{\alpha t}}\right) \quad (8)$$

$$C_i(z, t) = e^{-z\lambda_i} \left( 1 + \operatorname{erf}\left[\lambda_i \sqrt{\alpha t} - \frac{z}{2\sqrt{\alpha t}}\right] \right) - e^{z\lambda_i} \operatorname{erfc}\left[\lambda_i \sqrt{\alpha t} + \frac{z}{2\sqrt{\alpha t}}\right] \quad (9)$$

Using the material properties for aluminium and process parameters observed experimentally, it is possible to draw the temperature distribution using Eq. 7 (see Fig. 4).

From this, it is possible to define the shape of a crater by defining an isotherm on the figure which is usually chosen as the melting temperature affected by a coefficient in order to take into account material that is melted but not flushed away.

### 2.1.3 Electro-mechanical modeling

An electro-mechanical model for EDM is convenient in the case of short spark duration as the material isn't sufficiently exposed in order to heat. Singh et al. [13] have proposed such a model exploiting the stress induced by electrostatic forces. This model is based on the assumption that the

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Temperature distribution from Beck's model</figcaption></figure>

plasma region is electrically neutral as the positive ion density is nearly equal to the electron density. The potential difference can't be accommodated in this neutral region creating a thin sheath near the cathode where there is a charge imbalance and where the potential gradient is formed. The presence of such a sheath creates a strong electrical field at the cathode that induces a negative charge. This surface charge is pulled outwards by the field resulting in stress. The research concludes with the observation that at low duration pulses, the crater depth is independent of the spark duration.

### 2.1.4 Geometric modeling

Geometric modeling is opposed to the previous models in the fact that it is based on simulation and not formulation. The concept lies in knowing the geometry of the tool and the workpiece and that one can derive the gap separating these two elements. This enables us to determine the probability that a discharge will occur between two points. Once the most probable discharge location is determined,
removal on both the electrode and the workpiece can be assessed. Most models deal with two dimensional axisymmetric geometries [14]. A 2D geometry profile is extracted from the 3D geometry. The profile is then expressed as a discrete function, i.e., a set of points. The algorithm then proceeds to determine the discharge location and deleting the relevant nodes.

### 2.1.5 Conclusion

Most of the theoretical models mainly tackle with the modeling of a single discharge in order to derive the resulting crater shape and dimensions. None of these methods provide a satisfying solution for the simulation of an entire process mainly due to the fact that the calculations involved are very long for a single crater. Considering that hundreds of thousands of discharges can occur in a single process, the computation time factor is quite relevant. A recent geometrical effort involving the use of a Z-map has been developed specifically for the case of µEDM milling with good results [15]. However the inherent limitations of a Z-map have led to the idea of developing a purely geometrical method mainly focused on providing an accurate representation of the resulting geometries as well as decent computation speed.

From this analysis, it clearly appears that there is a need for developing a 3D surface deformation technique which should not have to solve complex PDE at each step. Hence, as illustrated on Fig. 4, such an approach has to be able to model one-by-one the apparition of sphere-like craters so as to follow an isotherm whose shape are driven by pre-processed physical parameters of the µEDM process.

## 2.2 Surface deformation techniques

As mentioned previously, simulation efforts of the µEDM process have focused on the use of a cloud of points (usually set in a regular grid). In order to eliminate the issue of interpolation approximation, the proposed approach will make use of a parametric surface as the deformable geometry. As a result, the different deformation techniques that are available to that specific geometric representation will be discussed. Of course, deformation techniques have been proposed in many different fields and on many representations [16, 17]. Additionally, with the purpose of avoiding expensive computations as much as possible, only purely geometric deformation techniques will be studied. Approaches that make use of physical equations or properties won't be discussed [18, 19].

The different approaches will be classified using a few criteria relevant to the problem at hand.

* Algorithmic speed considering the high number of discharges that occur during a µEDM manufacturing process (the exact number depends on the parameters being used but can easily be in the millions), the chosen method must be fast in order to be repeated numerously.

* Local/global control this criterion illustrates the fact that it is needed to be able to deform a very specific location of our geometry without affecting the rest of it. The global deformation will result from a set of local deformations.

* Model preservation in order to avoid added complexity in an iterative process, it is desirable that no change in the initial model be made. For example if the initial number of patches describing the geometry is one, a method preserving the model will keep it that way.

* Shape control this criterion is linked to the variety of deformation that are available.

Criteria linked to the user are not considered here. Interactivity is not a concern since the considered process is to be fully automatic. Criteria characterizing the fact that the method is applicable to multiple trimmed patches is not to be used since in our case, a unique patch will be used. Also, the shape control is not a critical criterion since the main interest lies in the insertion of simple spherical shapes several hundreds of thousands of times per simulation.

The definition of a NURBS surface involves three elements: the knot vectors, the control points and their weights. Changing any of these parameters will have an effect on the resulting surface. However, working directly with these can be tedious and designers usually tend to work directly with the surface itself and indirectly change the position of control points or the value of weights. The modification of knots is hardly used, even indirectly, as it is highly less intuitive than weight modification or control point repositioning. Anyhow, three categories can be identified:

* Control point repositioning the most common way of modifying NURBS curves or surfaces. The user generally specifies one or multiple curve or surface points to be displaced. Since the problem is often under-constrained, an energy function to be minimized is added to this optimization problem [16, 20, 21].

* Weight or knot modification this method is used in the specific case where the displacement of a point towards or further away from a control point is required. While local control can be increased similarly to control point repositioning, the nature of this method prevents it to be used for shape-specific deformations. Also the modification of knots is not intuitive and the deformations are
<table><caption>**Table 1** Performance of some methods in regards to various criteria</caption><thead><tr><th></th><th>Algorithmic speed</th><th>Local control</th><th>Model preservation</th><th>Shape control</th></tr></thead><tbody><tr><td>Control point repositioning</td><td>⊕</td><td>⊕⊕</td><td>⊕⊕</td><td>⊕</td></tr><tr><td>Weight, knot sequence modification</td><td>⊕⊕</td><td>⊕⊕</td><td>⊕⊕</td><td>⊕</td></tr><tr><td>Surface warping</td><td>⊕</td><td>⊕⊕</td><td>⊕⊕</td><td>⊕⊕</td></tr></tbody></table>

indirectly performed. The use of non-linear constraints slows down the deformation process [22–26].

* Surface warping this deformation strategy relies on the notion of control point repositioning. However, instead of specifying a set of constraints linked to points’ positions or derivatives’ values, the control points are moved in respect to the following general formula:

$$ \mathbf{P}'_{i,j} = \mathbf{P}_{i,j} + f(i,j) \cdot \mathbf{N}(i,j) \quad (10) $$

where $f(i, j)$ is the warp function and $N(i, j)$ is the warp direction. A certain number of common strategies exist such as surface flattening or bending [27–29].

A comparison of the previously described methods is presented in Table 1.

While control point repositioning with the use of a certain number of constraints seems to be a very good candidate, the shape control is lacking. It would be possible to specify an exact shape for the deformation but at the cost of using a high number of constraints which adds to the algorithmic complexity. Weight and/or knot sequence modification being too restrictive, surface warping remains the viable candidate.

As a conclusion, local surface warping will be used to insert one-by-one a large number of spherical craters in a heavily refined (1000 × 1000 control points) untrimmed NURBS patch of degree 3 in each direction. Using such a multi-resolution approach, it will be possible to get the tool’s wear as well as the workpiece’s roughness.

# 3 The μEDM simulation framework

## 3.1 Simulation process overview

In the proposed approach, both the geometry of the tool and the geometry of the workpiece are defined by means of NURBS patches. To allow the insertion of thousands of craters, the surfaces of the tool $S_t$ and workpiece $S_w$ are heavily refined using the Boehm’s knot insertion algorithm [30]. As a result the surfaces’ control points will be a lot closer to it, hence local control is significantly increased. At each step of the insertion process the location of each crater (one on each electrode) is determined while identifying the shortest distance between the tool and the workpiece since it is considered that

the electrical spark will happen on the less resistive path, i.e., the shortest one. Minimum distance computations are done using an optimization method known as particle swarm optimization (PSO) which is a simple numerical optimizer that does not require the use of the gradient of the objective function [12]. A crater is then inserted in each of those locations by moving the surrounding control points. As explained in the literature review, the shape of a crater can be assimilated to a part of sphere. If the computed minimum distance $d$ exceeds the value of the minimum distance required for a spark to appear (known in EDM as the gap distance $M_g$) then the tool is moved down along the $z$ axis with an increment of $\Delta z$. Otherwise, if the computed minimum distance is smaller than $M_g$ then the PSO algorithm returns four values ($u_t, v_t$) and ($u_w, v_w$) corresponding to the parametric coordinates of the craters’ centres respectively on the tool (subscript $t$) and workpiece (subscript $w$). The algorithm then moves the control points located in the surrounding of the craters’ centres so that two craters of volumes $\Delta V_t$ and $\Delta V_w$ are inserted into the tool and workpiece. The deformation technique is similar to surface warping considered as a geometric deformation technique [30]. The process ends when the desired depth $D_g$ is met. The overall algorithm can be described in pseudo-code form as in Fig. 5. The different steps of the crater insertion process are further detailed in the next subsections.

## 3.2 Volume enclosed by a NURBS patch

In the proposed algorithm, it is mandatory to control the volume $\Delta V_e$ removed on the electrode $e$ (with $e \in \{t, w\}$) at each iteration of the simulation process. To do so, the idea is to compute the difference between the initial volume $V_{e0}$ and the volume at a given iteration $V_e$. Generally speaking, the Green–Ostrogradsky’s theorem states that for a vector field $\mathbf{F}$ and a region $V$ enclosed by a surface $S$:

$$ \int_V \mathrm{div} \, \mathbf{F} \, dV = \oint_S \mathbf{F} \cdot \mathrm{d}\mathbf{S} \quad (11) $$

In the proposed approach, the surface $S_e$ to be deformed is defined by a single NURBS patch of degree $p_e$ in $u$ and $q_e$ in $v$ (Fig. 6a) and defined by its parametric Eq. [30]:

$$ \begin{aligned} \mathbf{S}_e(u, v) &= \frac{\sum_{i=0}^{n_e} \sum_{j=0}^{m_e} w_{ij} N_{ip}(u) N_{jq}(v) \mathbf{P}_{ij}}{\sum_{i=0}^{n_e} \sum_{j=0}^{m_e} w_{ij} N_{ip}(u) N_{jq}(v)} \\ &= {}^t[x(u, v), y(u, v), z(u, v)] \end{aligned} \quad (12) $$
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Flowchart of the simulation's algorithmic process</figcaption></figure>

where $\mathbf{P}_{ij}$ represents one of the $(1+n_e) \times (1+m_e)$ control points to be moved, and $w_{ij}$ the weights associated to each control points to be moved. As explained in Sect. 2.2, the weights will be taken equal to 1 in the proposed approach

and won't be considered as unknowns of the deformation process.

The initial surfaces are, at the moment, defined by hand while paying attention to their parameterization. Future
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Enclosed volume computation</figcaption></figure>

<table><caption>**Table 2** Comparison of speed and accuracy for the sampling-based method</caption><thead><tr><th>Nb. points</th><th>Time (ms)</th><th>V<sub>w</sub> (μm<sup>3</sup>)</th><th>D<sub>b</sub> (%)</th><th>D<sub>f</sub> (%)</th></tr></thead><tbody><tr><td>25</td><td>24</td><td>283329</td><td>15.66</td><td>100</td></tr><tr><td>50</td><td>97</td><td>360872</td><td>7.42</td><td>24.74</td></tr><tr><td>100</td><td>379</td><td>344769</td><td>2.63</td><td>6.33</td></tr><tr><td>200</td><td>1508</td><td>334685</td><td>0.38</td><td>1.59</td></tr><tr><td>400</td><td>3397</td><td>336873</td><td>0.27</td><td>0.70</td></tr><tr><td>800</td><td>24652</td><td>335950</td><td>0.00</td><td>0.09</td></tr></tbody></table>

works might include interpolating the initial geometries from a different format (meshes, cloud of points...).

The surface $S_e$ not being a closed surface, the idea is to choose the field $\mathbf{F} = {}^t[0, 0, z]$ so that:

$$ \oint_S \mathbf{F} \cdot d\mathbf{S} = \int_{S_e} \mathbf{F} \cdot d\mathbf{S}_e = \int_{U_e} \mathbf{F} \cdot \left( \frac{\partial \mathbf{S}_e}{\partial u} \times \frac{\partial \mathbf{S}_e}{\partial v} \right) du dv \quad (13) $$

where $U_e$ is the patch's parametric domain. The patch is not trimmed in our simulation tool. Thus the enclosed volume $V_e$ becomes (Fig. 6b):

$$ V_e = \int_{U_e} z \cdot \left( \frac{\partial x}{\partial u} \frac{\partial y}{\partial v} - \frac{\partial y}{\partial u} \frac{\partial x}{\partial v} \right) du dv \quad (14) $$

It can then be approximated while discretizing the parametric space $U_e$ with $(N_e^u \times N_e^v)$ sample points. Thus, for a given sample point $\mathbf{S}_{k,h} = \mathbf{S}_e(u_k, v_h)$, and a

<table><caption>**Table 3** Comparison between sampling-based (800 points in each direction) and control points-based (1000 points in $u$ and $v$) methods</caption><thead><tr><th>Method</th><th>Time (ms)</th><th>V<sub>w</sub>(μm<sup>3</sup>)</th><th>D<sub>b</sub> (%)</th><th>D<sub>f</sub> (%)</th></tr></thead><tbody><tr><td>Sampling</td><td>3397</td><td>335950</td><td>100</td><td>1.79</td></tr><tr><td>Control points</td><td>61</td><td>335137</td><td>0.24</td><td>100</td></tr></tbody></table>

neighborhood $[u_k, u_{k+1}] \times [v_h, v_{h+1}]$, the surface $S_e$ can be approximated locally by a quadrangle which can be decomposed in two triangles $T_{k,h}^1$ and $T_{k,h}^2$ and the enclosed volume becomes (Fig. 6c):

$$ V_e = \sum_{k=1}^{N_e^u-1} \sum_{h=1}^{N_e^v-1} \left[ V_e(T_{k,h}^1) + V_e(T_{k,h}^2) \right] \quad (15) $$

where the volume $V_e(T_{k,h}^1)$ enclosed by $T_{k,h}^1$ can be computed using Eq. 14 and a simple parameterization of the triangle (Fig. 6c):

$$ V_e = \int_{U_e} \frac{1}{3} (s_{k,h}^z + s_{k+1,h}^z + s_{k,h+1}^z) \frac{(e_k^x e_h^y - e_k^y e_h^x)}{||e_k \times e_h||} du dv \quad (16) $$

With $e_k = s_{k+1,h} - s_{k,h}$ and $e_h = s_{k,h+1} - s_{k,h}$

$$ V(T_{k,h}^1) = \frac{1}{6} (e_k^x e_h^y - e_k^y e_h^x) (s_{k,h}^z + s_{k+1,h}^z + s_{k,h+1}^z) \quad (17) $$

Here, $x, y$ and $z$ stands for the coordinates of the vectors represented on figure (Fig. 6c). A similar formula applies for the triangle $T_{k,h}^2$.

In order to be able to tune the discretization parameters $N_w^u$ and $N_w^v$, several sampling have been performed on a final surface $S_w$ and the results are reported in Table 2. In this table, $D_b$ is the ratio to the best estimation, and $D_f$ the ratio to the fastest estimation. As expected, doubling the number of sample points in each direction leads to four times more sample points over the entire surface, hence the total computation time is also greater by a factor of four. Overall, it is noticed that acceptable results (those that would be within 5 % of the true value of the volume) would require computational time of several hundreds of milliseconds.

Even if the simulation tool doesn't need to sample the whole surface but just the portion being deformed, it is still interesting to explore other strategies. An alternative would be to consider that the control points are sufficiently close to the real shape to be used in an accurate calculation of the surface's enclosed volume. In this case, formulas similar to 15 and 17 can be obtained directly from the control points and not anymore from the sampling points (Fig. 6a). A comparison between the sampling-based and control points-based methods is presented in Table 3. From those results, it is clear that the control points-based method is
<figure><img src="image_9.png" /><figcaption>Fig. 7 Definition of the warping vector $\omega_w^{[k]}$ for a crater to appear on the workpiece</figcaption></figure>

much faster than the sampling approach without however reducing a lot the quality of the approximation. Thus, this method has been adopted since the computation of the enclosed volume will be performed hundreds of thousands of times during the μEDM simulation.

### 3.3 Volume to be removed for each crater

Each electrical spark transfers a certain amount of energy to the tool, the workpiece and the dielectric fluid. Here, it is considered that the amount of energy brought to each element is the same at each spark. As such, it is desirable to remove the same volumes $\Delta V_t$ and $\Delta V_w$ when simulating the insertion of all the craters. These volumes are experimentally obtained by measuring the mean radius $\bar{R}_e$ and mean depth $\bar{D}_e$ of actual craters (see Sect. 4.1). Then, considering that the craters are domes, the volumes to be removed $V_e^r$ are computed using the following formula:

$$V_e^r = \pi \frac{\bar{D}_e}{6} (3\bar{R}_e^2 + \bar{D}_e^2), \quad e \in \{t, w\} \qquad (18)$$

From these volumes and domes, the support spheres can be identified, i.e., the spheres of radii $R_e$ equal to the dome's radius. As explained, these two radii remain constant for the two surfaces for the crater-by-crater simulation.

### 3.4 Craters insertion

Following the flowchart of Fig. 5, for a given depth of the tool, if the minimum distance is smaller than the gap distance, the craters insertion process starts. Craters are inserted one by one on each surface $S_e(e \in \{t, w\})$ and centered on $S_e(u_e, v_e)$. For sake of clarity, the superscript [k] has not been put on the parametric coordinates $u_e$ and $v_e$ even if these values refer to the kth craters (one on each

<figure><img src="image_10.png" /><figcaption>Fig. 8 Definition of the support sphere and removed volume $\Delta V_w^{[k]}$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9 Plane $\Pi_w^{[k]}$ definition and a control point's projection</figcaption></figure>

surface). First, to identify the displacement directions, the two warping unit vectors are computed as follows:

$$\omega_e^{[k]} = sg(e) \frac{\mathbf{S}_t(u_t, v_t) - \mathbf{S}_w(u_w, v_w)}{\|\mathbf{S}_t(u_t, v_t) - \mathbf{S}_w(u_w, v_w)\|} \qquad (19)$$

$$\text{with } sg(e) = \begin{cases} 1 & \text{for } e = t \\ -1 & \text{for } e = w \end{cases}$$

Figure 7 represents a two-dimensional version of the process after having found the minimum distance. For sake of clarity, solely the part of the workpiece is represented even if the same strategy applies to the tool. The next step consists in identifying which control points need to be moved in the surrounding of the two points $S_e(u_e, v_e)$.

Once the radii of the two spheres identified (Sect. 3.3), the location of the spheres' centres has to be computed (one for the tool and one for the workpiece). As illustrated on Fig. 8, the centre of the sphere lies on the spark line. Its
<figure><img src="image_12.png" /><figcaption>**Table 4** Machining parameters</figcaption></figure>

<table><caption>**Table 5** Experimental results</caption><thead><tr><th>Experiment</th><th>1</th><th>2</th></tr></thead><tbody><tr><td>Energy level (index)</td><td>300</td><td>13</td></tr><tr><td>Voltage (V)</td><td>60</td><td>60</td></tr><tr><td>Current (index)</td><td>20</td><td>20</td></tr><tr><td>Time on (ms)</td><td>5</td><td>5</td></tr></tbody></table>

<table><thead><tr><th>Experiment</th><th>1</th><th>2</th></tr></thead><tbody><tr><td>Hole depth (μm)</td><td>50.8</td><td>50.2</td></tr><tr><td>Tool vertical wear (μm)</td><td>12.5</td><td>11.3</td></tr><tr><td>Roughness Ra (μm)</td><td>1.27</td><td>0.82</td></tr><tr><td>Workpiece crater diameter (μm)</td><td>15</td><td>3</td></tr><tr><td>Workpiece crater depth (μm)</td><td>3</td><td>1</td></tr></tbody></table>

exact position depends on the volume $V_e^r$ that needs to be removed. In order to find the location, an iterative dichotomy method (also known as binary search or bisection method) is used. At each step, the surface $S_e$ is deformed and the intersecting volume (the hashed part of Fig. 8) between the sphere and the surface is computed while using the previously introduced formula $(\Delta V_e^{[k]} = V_e^{[k]} - V_{e0}^{[k]})$. If the volume obtained is smaller than $V_e^r$ the sphere is moved towards the surface, and if it is bigger it is moved away from it. The process carries on until the obtained volume $\Delta V_e^{[k]}$ falls within a specific tolerance $T_v$. Once the $C_e^{[k]}$ adequate positions are found, it is possible to determine the $N_e^{[k]}$ control points of the two surfaces that need to be moved. This is done by computing for each control point the distance that separates them from the centre of the sphere. If the distance is smaller than the radius of the sphere, the control point is added to the list of points to be displaced. At the end, two lists of control points are obtained.

In order to displace the control points to mimic the shape of a sphere, a reference is needed. Let $\Pi_e^{[k]}$ be the plane that includes the centre of the sphere $C_e^{[k]}$ and that has $\omega_e^{[k]}$ as normal vector. Then, the new position of all the control points $\mathbf{P}_{e,j}^{[k]}$ to be moved ($e \in \{t, w\}$ and $j \in \{1...N_e^{[k]}\}$) are computed as follows (Fig. 9):

$$ \tilde{\mathbf{P}}_{e,j}^{[k]} = \mathbf{P}_{e,j}^{[k]} + f_{e,j}^{[k]}(r_{e,j}^{[k]}).\omega_e^{[k]} $$ 
with $f_{e,j}^{[k]}(r_{e,j}^{[k]}) = \sqrt{R_e^2 - r_{e,j}^{[k]}^2} - [\mathbf{P}_{e,j}^{[k]} - \pi_{e,j}^{[k]}].\omega_e^{[k]}$ (20) 
and $r_{e,j}^{[k]} = \|\pi_{e,j}^{[k]} - \mathbf{C}_e^{[k]}\|$, 
where $\pi_{e,j}^{[k]}$ is the projection of $\mathbf{P}_{e,j}^{[k]}$ on the plane $\Pi_i^{[k]}$. In this formulation, a unique index $j$ has been used to identify the $j$th control point $\mathbf{P}_{e,j}^{[k]}$ of $S_e$ to be moved to generate the $k$th crater.

This process is repeated iteratively until no more craters can be inserted for the actual depth. Then, the tool is moved down along the z axis with an increment of $\Delta_z$ and the craters insertion process starts again (Fig. 5).

## 4 Experimental validation

Three experiments were conducted in order to assess the accuracy of the simulation tool. The two first experiments were designed to give a rough idea of the simulation's performances through indirect measurements of the surfaces' quality (vertical wear and roughness) while the third used direct surfaces comparison with the use of the Hausdorff metric.

### 4.1 Measurements

For the first two experiments, experimental results related to roughness and crater volumes were measured with the help of a white light interferometer. The experimental vertical tool wear was measured on the μ EDM machine using the contact detection feature. On the other side, measurements on the NURBS patches involved in the simulation were realised by numerically implementing the formula for average roughness [31]. Several measurements in different locations were made then averaged. The last experiment's resulting surfaces were measured with the help of a micro-tomographer [32] that made it possible to obtain meshed three-dimensional models of the tool and workpiece after the experiment.
## 4.2 Vertical wear and surface roughness evaluation

As a validation of the new simulation process, two experiments were conducted on a Sarix SX-200 µEDM machine equipped with a wire dressing unit. A tungsten carbide rod, with a nominal diameter of 290 µm was used as tool electrode. Ultra-Fine Grained aluminium (Al1070) with an average grain size of 0.6 µm was chosen as workpiece material to minimise the material's inhomogeneity while aiming at improving µEDM predictability, as suggested in [33].

Using the wire dressing unit, the tip of the tool electrode was machined flat for the first experiment, while for the second experiment a curved shape was introduced (Fig. 10a). The electrodes were then used to erode the workpiece down to a 50 µm depth objective. Machining parameters are shown in Table 4. The final shape of the tool after machining is depicted on Fig. 10b for the experiment 2.

With the help of a white light interferometer, experimental crater dimensions were determined and tabulated in Table 5. The knowledge of those dimensions was then used to infer a value for the workpiece craters’ volume using Eq. 18. The volume of the tool craters’ was determined with the use of the wear ratio derived from the hole depth and the tool vertical wear. As a result, in the case of experiment 1, the target volume $\Delta V_e$ expected to be removed per crater was 279.20 $\mu m^3$ for the workpiece and 71.69 $\mu m^3$ for the tool. In the case of experiment 2, the target volume $\Delta V_e$ expected to be removed per crater was 4.05 $\mu m^3$ for the workpiece and 1.04 $\mu m^3$ for the tool.

For the simulation of the two experiments, both the tool and workpiece have been modeled by an untrimmed NURBS patch of degree 3 in each direction and defined by a network of 1000 × 1000 control points.

In a first time, three simulations of Experiment 1 were performed using different tolerances $T_v$ to assess the influence of the computational precision on roughness results (Table 6). Measurements on the NURBS patches involved in the simulation were realised by numerically implementing the formula [31] for the average roughness $R_a$. Several measurements in different locations were made and averaged. From those results, it is clear that the tolerance level has a significant influence on the final simulated roughness.

In a second time, one simulation was performed for Experiment 2 to assess agreements in terms of geometry deformation due to tool wear, as well as achieved roughness. Based on the previous results, this simulation used a tolerance level $T_v = 1\%$. Figure 11 presents an overview of the evolution of the workpiece during the simulation rendered as a wireframe of the NURBS' surface. The total runtime of the process was around 20 h on a personal

computer (i5 3570 K at 4.5 Ghz) and the number of craters inserted was 375494. It means that the proposed crater-by-crater insertion algorithm (Sect. 3.4) has required 375494 local surface warping to generate the final shape of the workpiece (Fig. 11g). Actually, in our context, the need for simulating accurately the local surface deformation (to get the tool wear and workpiece roughness) increases a lot the computation time. Of course, it cannot be compared to other applications targeting a global deformation inducing a single surface warping unable to capture shape deformations at the scale of the roughness. One can finally notice that, at the end of the simulation, not all the craters are visible on the outer surface since most of them have been iteratively removed during the crater-by-crater iterative deformation process. A detail of the craters is provided in Fig. 12.

The simulation results are in accordance with those of the second experiment for the values of vertical tool wear and roughness (Table 7) as the experimental and simulated surfaces deviate by 0.89 % for the tool vertical wear and 6.09 % for the roughness. Additionally, the tool's profile can be compared with the simulated one (Fig. 13). It appears that the simulation does not lead to the dissymmetry present in the experimentation. This can be due to several factors that haven't been included yet in the simulation, notably the influence of the dielectric flow and an eventual alignment error while machining the tool by wire EDM and/or during the die-sinking process.

During the simulation process, a certain number of difficulties had to be managed with the displacement of control points. An evident case to avoid is displacing control points leading to self-intersecting geometries. This is often the case when dealing with the sides. As a result, the geometries must be taken into account when defining the warping vector. This will be optimised in the future.

Another issue linked to the surfaces' parameterization occurred. A simple example highlighting this issue is to consider a scenario beginning with a flat surface. As it is being deformed, the control points' displacements will lead to some areas of the surface having a smaller density of control points than others. This phenomenon is visible on Fig. 11g on the outer flanks of the workpiece. A solution would be to regularly re-parameterize the surface to avoid this. It would however add to the computation time. This will be further studied in the future.

## 4.3 Shape difference assessment

Following the previous experiments, an additional one was devised with the objective of directly comparing three-dimensional meshes of the experimental results with meshed NURBS surfaces from the simulation.
<figure><img src="image_13.png" /><figcaption>**Fig. 11** Evolution of the workpiece's simulated shape: **a** initial configuration, **b** 1417 craters, **c** 24238, **d** 64543, **e** 251425, **f** 285401 and **g** 374294</figcaption></figure>

The tool used is of a spherical shape obtained through wire-dressing of a tungsten carbide electrode of nominal diameter 300 µm. The tool was then measured through the use of micro-tomography [32] (with a resolution of 1 µm) and the result exported as a mesh visible in Fig. 14. The nominal diameter of the sphere is 250 µm while the measured diameter through interpolation with a sphere results in 246 µm. The experimental machining parameters are given in Table 8. The feature machined on the workpiece was also measured in the same manner and allowed the meas-

<table><caption>**Table 6** Effect of volume removal precision on roughness (Experiment 1)</caption><thead><tr><th>Tolerance T<sub>v</sub> (%)</th><th>10</th><th>5</th><th>1</th></tr></thead><tbody><tr><td>Volume removed (µm<sup>3</sup>)</td><td>172058</td><td>169978</td><td>167586</td></tr><tr><td>Average volume per crater (µm<sup>3</sup>)</td><td>286.76</td><td>283.30</td><td>279.31</td></tr><tr><td>Roughness Ra (µm)</td><td>1.87</td><td>1.88</td><td>1.43</td></tr></tbody></table>

urement of the craters directly from the meshed model. The dimensions of those are given in Table 9.

<figure><img src="image_14.png" /><figcaption>**Fig. 12** Detail of the craters on the workpiece. The scale is given in the bottom-left</figcaption></figure>

<table><caption>**Table 7** Comparison experimental vs. simulation results (Experiment 2)</caption><thead><tr><th>Characteristic</th><th>Experi. (µm)</th><th>Simu. (µm)</th><th>Deviation (%)</th></tr></thead><tbody><tr><td>Tool vertical wear</td><td>11.3</td><td>11.2</td><td>0.89</td></tr><tr><td>Roughness Ra</td><td>0.82</td><td>0.87</td><td>6.09</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Fig. 13** Tool profiles from experimentation (a) and simulation (b). Profiles differences with simulation extra volume in black (c)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14** Meshed model of the measuring of the initial tool</figcaption></figure>

A sphere was subsequently modelled in the simulation software as a single NURBS patch and refined by knot insertions. The resulting patch is defined with a network of a 1000 × 1000 control points. Using the same tolerance parameter $T_v$ as for the experiments described in Sect. 4.2 and the experimental crater dimensions measured on the 3D models, a simulation was started and ended after 34 h and 18201 discharges. The resulting NURBS surfaces visible in Fig. 15 were then converted as three-dimensional
<figure><img src="image_17.png" /><figcaption>Fig. 15 NURBS surfaces of the tool (left) and the workpiece (right)</figcaption></figure>

meshes (a million points were sampled uniformly on each patch) and processed with the use of Meshlab (http://mesh-lab.sourceforge.net/) and GOM inspect (http://www.gom.com) in order to cut and align those meshes with the experimental ones using the Gaussian best fit method.

In order to compare those surfaces, the Hausdorff metric was used [34] as defined in Eq. 21 where X and Y are two sets of points.

$$d_H(X, Y) = \max\{\sup_{x \in X} \inf_{y \in Y} d(x, y), \sup_{y \in Y} \inf_{x \in X} d(x, y)\} \quad (21)$$

The results are tabulated in Table 10 and can be observed in Fig. 16.

Those results can be interpreted in comparison to the bounding box diagonal length of each model in order to appreciate the scale of the Hausdorff distance when compared to the dimensions of the models. The bounding box diagonal of the tool was 366.84 µm while the workpiece's was 366.64 µm. Relatively to those dimensions, the two meshes are within 0.4 % on average for the tool and 0.8 % for the workpiece. In light of those results, the simulation method can be considered validated as it provides with a sensibly accurate surface profile of the tool and workpiece after machining. Additionally, the maximum distance is due to debris that have resolidified during the

process (observable in Fig. 16 as the bright blue spot on the workpiece for example) and that the simulation does not account for.

## 5 Conclusions

In order to overcome issues linked with the wear phenomenon in µEDM, it is important to be able to predict said wear. A new method of simulation involving NURBS surfaces was presented. It relies on an iterative crater-by-crater insertion algorithm. Through hundreds of thousands of local surface warping, the method is able to compute the global as well as the local shapes of the tool and workpiece. At each step, the warping vector and function are computed so as to be able to generate a spherical crater whose volume is also controlled. This is a completely different approach from the existing surface deformation techniques found in the literature. While acting very locally to simulate the real µEDM manufacturing process, not only is it possible to evaluate the tool wear from the overall final shape at a low resolution level, but also to estimate the workpiece roughness from the high resolution level.

The overall simulation framework being now set up, several optimization of this promising approach are foreseen. Although some discrepancies were found between experimental data and simulation results, those values remain acceptably close (deviation of 0.89 % for the tool vertical wear and 6.09 % for the roughness and a low Hausdorff metric for the third experiment). As discussed, crater dimensions measurements were of paramount importance in the final results and are extremely important if a good accuracy was to be obtained. Another approach could be to use theoretical values for crater dimensions in lieu of the experimental ones or, alternatively, derive those values from a database built

<figure><img src="image_18.png" /><figcaption>Fig. 16 Hausdorff metric mapped on the meshes of the experimental tool (left) and workpiece (right). The histogram is displayed on the left of each model</figcaption></figure>
<table><caption>**Table 8** Machining parameters</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Energy level (index)</td><td>200</td></tr><tr><td>Voltage (V)</td><td>90</td></tr><tr><td>Current (index)</td><td>20</td></tr><tr><td>Time on (ms)</td><td>5</td></tr><tr><td>Objective depth (μm)</td><td>100</td></tr></tbody></table>

<table><caption>**Table 9** Experimental craters dimensions</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Workpiece crater diameter (μm)</td><td>13.30</td></tr><tr><td>Workpiece crater depth (μm)</td><td>4.42</td></tr><tr><td>Tool crater diameter (μm)</td><td>12.40</td></tr><tr><td>Tool crater depth (μm)</td><td>4.39</td></tr></tbody></table>

<table><caption>**Table 10** Hausdorff metric results</caption><thead><tr><td></td><td>Tool</td><td>Workpiece</td></tr></thead><tbody><tr><td>d<sub>H</sub> min (μm)</td><td>0.000000</td><td>0.000107</td></tr><tr><td>d<sub>H</sub> max (μm)</td><td>8.629291</td><td>14.886533</td></tr><tr><td>d<sub>H</sub> mean (μm)</td><td>1.449477</td><td>3.073571</td></tr><tr><td>d<sub>H</sub> RMS (μm)</td><td>2.521132</td><td>3.626015</td></tr></tbody></table>

from electrical measurements of the process [35]. Efforts should also be put into preserving the surfaces' parameterization as well as reducing the computational times, even if this last criterion is not crucial in the considered context. For example, the way the volume computation is performed could be optimized so as to only consider local variations of the volume, and not global variations, on a NURBS surface defined by 1000 × 1000 control points. The minimum distance computation could also be improved while considering the multi-resolution character of the simulation and avoid looking for a global minimum distance on the entire patch.

Finally, the new simulation tool will be integrated in a shape optimization loop so as to be able to predict the optimal shape of the initial tool for a given shape of the final workpiece.

**Acknowledgments** This work was supported by the Engineering and Physical Sciences Research Council [EP/J004901/1]. The authors are grateful to Prof. Maurizio Santini for the micro-tomography analysis.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

## References

1. Ho KH, Newman ST (2003) State of the art electrical discharge machining (EDM). Int J Mach Tools Manuf 43(13):1287-1300

2. Dhanik S, Joshi SS (2005) Modeling of a single resistance capacitance pulse discharge in micro-electro discharge machining. J Manuf Sci Eng Trans ASME 127:759-767

3. Bissacco G, Hansen HN, Tristo G (2011) Valentinic J (2011) Feasibility of wear compensation in micro edm milling based on discharge counting and discharge population characterization. CIRP Ann 60:231-234

4. Rajurkar KP, Yu ZY (2000) 3d micro-edm using cad/cam. CIRP Ann Manuf Technol 49(1):127-130

5. Tao J, Ni J, Shih AJ (2012) Modeling of the anode crater formation in electrical discharge. Mach J Manuf Sci Eng 134(1):011002-1

6. Bigot S, Surleraux A, Bissacco G, Valenticic J (2012) A new modelling framework for die-sinking micro edm. In: Proceedings of the 9th international conference on multi-material manufacture, pp 51-55

7. Descoeudres A, Hollenstein C, Walder G, Perez R (2005) Time-resolved imaging and spatially-resolved spectroscopy of electrical discharge machining plasma. J Appl Phys 38:4066-4073

8. Snoeys R, VanDijck FS (1971) Investigation of electro discharge machining operations by means of thermo-mathematical model. CIRP Ann 20:35-37

9. Snoeys R, VanDijck FS, Peters J (1972) Plasma channel diameter growth affects stock removal in edm. CIRP Ann 21:39-40

10. VanDijck FS, Dutre WL (1974) Heat conduction model for the calculation of the volume of molten metal in electric discharges. J Phys D 7:899-910

11. DiBitonto DD, Eubank PT, Patel MR, Barrufet MA (1989) Theoretical models of the electrical discharge machining process. I. A simple cathode erosion model. J Appl Phys 66:4095-4103

12. Beck JV (1981) Transient temperatures in a semi-infinite cylinder heated by a disk heat source. Int J Heat Mass Transf 24(10):1631-1640

13. Singh A, Ghosh A (1999) A thermo-electric model of material removing during electric discharge machining. Int J Mach Tools Manuf 39:669-682

14. Tricarico C, Delpretti R, Dauw DF (1988) Geometrical simulation of the edm die-sinking process. CIRP Ann Manuf Technol 37(1):191-196

15. Heo S, Jeong YH, Min B-K, Lee SJ (2009) Virtual edm simulator: Three-dimensional geometric simulation of micro-edm milling processes. Int J Mach Tools Manuf 49(12-13):1029-1034

16. Pernot J-P (2004) Fully free form deformation features for aesthetic and engineering designs. PhD thesis, Universita di Genova and Institut National Polytechnique de Grenoble

17. Zhang D, Lu G (2004) Review of shape representation and description techniques. Pattern Recognit 37(1):1-19

18. Nealen A, Mueller M, Keiser R, Boxerman E, Carlson M (2005) Physically based deformable models in computer graphics. Comput Graph Forum 25(4):809-836

19. Botsch M, Sorkine O (2004) On linear variational surface deformation methods, zhang, guojun lu, review of shape representation and description techniques. IEEE Trans Vis Comput Graph Pattern Recognit 14(1):213-230

20. Pourazady M, Xu X (2006) Direct manipulations of nurbs surfaces subjected to geometric constraints. Comput Graph 30(4):598-609

21. Xu G (2013) K. chuen Hui, W. bing Ge, and G. zhao Wang, Direct manipulation of free-form deformation using curve-pairs. Comput Aided Des 45(3):605-614

22. Juhasz I (1999) Weight-based shape modification of nurbs curves. Comput Aided Geom Des 16:377-383
23. Hui SM, Li YF, Ju T, Zhu X (2001) Modifying the shape of nurbs surfaces with geometric constraints. Comput Aided Des 33:903–912

24. Xie H, Qin H (2001) Automatic knot determination of nurbs for interactive geometric design. In: Shape modelling international conference, pp 267–276

25. Zheng W, Bo P, Liu Y, Wang W (2012) Fast b-spline curve fitting by l-bfgs. Comput Aided Geom Des 29(7):448–462

26. Hua S-M, Lib Y-F, Jua T, Zhua X (2001) Modifying the shape of nurbs surfaces with geometric constraints. Comput Aided Des 33(12):903–912

27. Zhang W, Leu MC (2009) A spatial warping method for free-form modeling based on a level-set method. Comput Aided Des 41(11):765–771

28. Wang CCL (2008) Wirewarping: A fast surface flattening approach with length-preserved feature curves. Comput Aided Des 40(3):381–395

29. Higashi M, Aoki N, Kaneko T (2002) Application of haptic navigation to modify free-form surfaces through specified points and curves. In: ACM symposium on solid modelling and applications (SM’02)

30. Piegl L, Tiller W (1996) The NURBS Book. Springer, Berlin

31. Degarmo EP, Black J (2003) Materials and processes in manufacturing. Wiley, New Jersey

32. Santini M, Guilizzoni M, Fest-Santini S (2013) X-ray computer microtomography for drop shape analysis and contact angle measurement. J Colloid Interface Sci 409:204–210

33. Bigot S, Bissacco G, Valentincic J (2011) Die-sinking micro edm for complex 3d structuring—research directions. In: International conference on multi-material micro manufacture

34. Cignoni P, Scopigno CRR (1998) Metro: measuring error on simplified surfaces. Comput Graph Forum 17(2):167–174

35. Ferri C, Ivanov A, Petrelli A (2008) Electrical measurements in uedm. J Micromech Microeng 18(8):13–25