Modeling of Material Removal Rate and Surface Roughness Generated during Electro-Discharge Machining

: This study reports on the numerical model development for the prediction of the material removal rate and surface roughness generated during electrical discharge machining (EDM). A simpliﬁed 2D numerical heat conduction equation along with additional assumptions, such as heat e ﬀ ect from previously generated crater on a subsequent crater and instantaneous evaporation of the workpiece, are considered. For the material removal rate, an axisymmetric rectangular domain was utilized, while for the surface roughness, a rectangular domain where every discharge resides at the end of previous crater was considered. Simulated results obtained by solving the heat equation based on a ﬁnite element scheme suggested that results are more realistic by considering instantaneous evaporation of the material from the workpiece and the e ﬀ ect of residual heat generated from each spark. Good agreement between our model and previously published data validated the newly proposed models and demonstrate that instantaneous evaporation, as well as residual heat, provide more realistic predictions of the EDM process.


Introduction
Electrical discharge machining (EDM) is a popular nonconventional machining process, which is capable of machining materials regardless of their hardness and strength. EDM is an electrothermal process, which erodes material as long as they possess some degree of electrical conductivity. In this process, due to an established electrical field between the workpiece and the tool, a plasma channel is generated in the gap. Ions and electrons in the plasma channel collide with the surfaces of the electrodes, which eventually results in localized heating of the surfaces. This highly intense heat flux results in a very high temperature at the surfaces, and consequently, melting and evaporation of material occurs. At the end of every discharge, part of the molten material is removed from the surface and is flushed away by the dielectric fluid, leaving behind a small crater. The remaining molten materials eventually resolidify and form a recast layer. Depending on the frequency, this process can happen several to hundreds of times in a second, thus creating multiple overlapping craters. Although the EDM process is a combination of different phenomena, such as dielectric break down, plasma formation, heat conduction, fluid dynamics, electrical phenomenon, chemical phenomenon, etc., thermal effects with a small error can be considered as the dominating phenomenon in the EDM process [1][2][3][4][5].
Many researchers have studied the EDM process through three different approaches; namely analytical, experimental, and numerical approaches. Van Dijck and Dutre [6] developed a two-dimensional by Salonitis et al. [29]. Their model predicts that both MRR and surface roughness increases with the increase of discharge current, voltage, and pulse duration.
In most of the modelling works of the EDM process, authors, for sake of simplicity, have tended to ignore many important aspects and parameters of the process. The evaporation of the material is one of the parameters usually ignored. In our work, we have focused, among the other aspects, on the evaporation of the material during discharge, and the remaining heat from previous discharges has been considered too. These parameters have been carefully evaluated and presented. Most important studies in the EDM process assume a single-discharge, while in an actual EDM process, multiple discharges happen that affect the EDM performance. In this study, new numerical models are developed for predicting the material removal rate and surface roughness of surfaces undergoing EDM as a function of process parameters separately. For MRR, a thermal model based on realistic assumptions, such as Gaussian heat distribution and instantaneous evaporation (IE), is developed. On the other hand, for surface roughness, the model is based on a simplified 2D numerical heat conduction equation where the heat effects of each discharge on subsequent discharges are taken into consideration.

Model Description
In this study, COMSOL Multiphysics 5.3 (COMSOL Inc., Burlington, MA, USA), a commercial finite element software package, was utilized for the modeling of heat transfer in the EDM process. For simplification, a 2-dimensional continuum was considered for the simulation. Two different set of simulations were performed, one for MRR and the other for surface roughness. For MRR prediction, a single discharge on the top surface of a cylindrical coordinate was implemented. However, for evaluation of the surface roughness, first, a single discharge at the edge of the rectangular computational domain was implemented. After creation of the first crater, the next discharge was positioned at the edge of the first crater. Afterward, the next discharge was applied at the edge of the second crater, and so on. This procedure for the first and second discharge is shown in Figure 1. In both sets of simulations, the melting temperature isotherm curve evaluated from the last time step of each simulation was considered as a boundary of the crater. Since the evaporated material was in gas form, it was expected to escape instantly, and evaporated material from the workpiece was constantly removed as the workpiece was heated, while the melted part was removed at the end of the simulation. A mesh with quadrilateral elements were generated using the COMSOL software. In order to achieve better accuracy for the evaluation of surface roughness, the top region of the computational domain, where craters were generated, smaller elements with a size of 0.25 µm were utilized. The total number of elements in the simulation was about 420,000.
In both sets of our numerical simulation, the following assumptions were made: 1. Gaussian heat flux with expanding plasma channel was considered.
The primary mode of heat transfer was conduction. Heat transfer via radiation and convection was ignored.
Effect of latent heat of fusion was taken into account. 6.
In every cycle, only one discharge was assumed. 7.
Required time for breakdown of dielectric fluid was considered equal to 500 ns as suggested by Almacinha et al. [22] The only distinction in assumptions of both set of simulations was that the calculation of the surface roughness considered the residual heat transfer from previous discharges in the formation of the next crater.

Governing Equations
The Fourier heat transfer equation is utilized to model the process [30]: where T is temperature (K), α is thermal diffusivity (m 2 /s), ∇ 2 is Laplace operator and t is time (s). Thermal diffusivity is defined as [30]: where k, ρ, and cp are the thermal conductivity (W m −1 K −1 ), density (kg/m 3 ), and heat capacity (J kg −1 ·K −1 )of the material, respectively. For the first model, the Laplace operator in cylindrical coordinates was considered: while for the second model, the simulation was performed in Cartesian coordinates:

Heat Source
Different techniques to model the heat source have been proposed in the literature. DiBitonto [7] proposed a point heat source, while many researchers [2,6,9,[31][32][33] have considered a uniform disk-like heat source. Results of the temperature measurement of the plasma channel have shown that the temperature distribution in the plasma channel is not uniform [18,34]. Many researchers [3,8,18,35,36] have considered a Gaussian heat distribution as a more realistic heat source model. A Gaussian heat distribution is based on the probability distribution function of a Gaussian distribution of variable r [37]: where σ is the standard deviation of the distribution. Since 99.7% of values reside in three standard deviation, σ = Rp/3, therefore: In an EDM process, p(r) is heat flux distribution (Qw). If the maximum heat flux is Qmax, therefore:

Governing Equations
The Fourier heat transfer equation is utilized to model the process [30]: where T is temperature (K), α is thermal diffusivity (m 2 /s), ∇ 2 is Laplace operator and t is time (s). Thermal diffusivity is defined as [30]: where k, ρ, and c p are the thermal conductivity (W m −1 K −1 ), density (kg/m 3 ), and heat capacity (J kg −1 K −1 )of the material, respectively. For the first model, the Laplace operator in cylindrical coordinates was considered: while for the second model, the simulation was performed in Cartesian coordinates:

Heat Source
Different techniques to model the heat source have been proposed in the literature. DiBitonto [7] proposed a point heat source, while many researchers [2,6,9,[31][32][33] have considered a uniform disk-like heat source. Results of the temperature measurement of the plasma channel have shown that the temperature distribution in the plasma channel is not uniform [18,34]. Many researchers [3,8,18,35,36] have considered a Gaussian heat distribution as a more realistic heat source model. A Gaussian heat distribution is based on the probability distribution function of a Gaussian distribution of variable r [37]: where σ is the standard deviation of the distribution. Since 99.7% of values reside in three standard deviation, σ = R p /3, therefore: In an EDM process, p(r) is heat flux distribution (Q w ). If the maximum heat flux is Q max , therefore: Machines 2019, 7, 47 5 of 17 The total energy applied over the surface can be calculated by integrating Equation (7): On the other hand, this value is equal to the power applied to the workpiece: It can be written as follows: By substitution of Q max from Equation (10) into Equation (8), the heat flux distribution is given as follows: when the discharge is applied at distance d from the origin, the Gaussian heat distribution is given by: where F c , V, I, R p are the energy fraction to the workpiece, discharge voltage, discharge current, and plasma radius, respectively. In Equation (12), d is the position where the spark is implemented. The most important factor for simulating the EDM process is perhaps the fraction of energy that is transferred to the workpiece. There is a huge discrepancy in suggested values for the fraction of energy from 14% to 50%. In this research, a value of 0.183, which is suggested by DiBitonto et al. [3] and Yeo et al. [28] was used. Later, a numerical study by Izquierdo [18] reconfirmed this value to be a proper value.

Plasma Radius
Both constant and expanding plasma heating radii have been used by different researchers [18][19][20]24,38,39]. In this research, the following semi-empirical relation for "equivalent heat input radius," proposed by Ikai and Hashiguchi [40], was used: (13) where I is the discharge current in amperes, t is discharge time, and R p is the plasma radius in micrometers. For solving Equation (1), appropriate boundary conditions were required. Boundary conditions were as follows (see Figure 1): For boundary A-B: For boundary B-C and C-D and D-A: where n is the normal vector.
Since the workpiece was in direct contact with a dielectric, the initial temperature in the computational domain was considered as the dielectric temperature: Machines 2019, 7, 47 6 of 17 Latent heat of melting consumed a significant amount of the energy supplied to the electrode. Therefore, the latent heat of melting was taken into account by means of the effective heat capacity c p,eff .
where λ is the latent heat of melting (J/kg) and ∆T is the temperature difference between the melting point of the material and dielectric temperature. Figure 2a,b illustrates the change in radius and the depth of the crater obtained from a single discharge simulation with a spark duration for different values of current. For low values of discharge time, the radius increased but then approached a constant value. Also, a higher current led to a larger radius. This trend can be explained by considering the dependence of the radius of the plasma channel on both the time and discharge current (Equation (12)). For short discharge times, the depth of the crater showed similar trends. However, during longer discharge times, depending on the discharge current, the depth of crater either approached a constant or decreased. This can be explained by taking into account the trends of heat flux and heat input radius: as the discharge time increased, the heat density on the surface decreased while the plasma radius continued to increase. In addition, when instantaneous evaporation was in effect, while the radius of the crater showed no significant change (Figure 2a), the depth of the crater showed higher values for a short discharge time. Due to a higher heat flux for a short discharge time, the evaporation of material caused deeper craters. As seen in Figure 2b, as the discharge time increased, the depth of crater with instantaneous evaporation approached the depth of the crater without instantaneous evaporation. Similarly, for higher values of the current, more evaporation was expected, and hence the model with IE showed a deeper crater. It is worth noting that there was no significant difference between the radius of the crater predicted by the model with and without instantaneous evaporation. Figure 2c-f shows increasing crater radii for increasing values of current.

Effects of Process Parameters on Material Removal Rate
For the evaluation of material removal rate, Equation (1) with a 2D axis-symmetry domain was solved. Steel with a thermal conductivity of 65 W/m K, density of 7545 kg/m 3 , heat capacity of 575 J/(kg K), latent heat of fusion of 247 kJ/kg, and melting point of 1808 K was considered in the simulation [23]. The shape of the crater can be evaluated by removing elements with temperature higher than the melting point. When the shape of the crater is known, volume of the crater can be evaluated using Equation (18): where D i is volume of a horizontal disk in the crater. The material removal rate MRR was calculated using the following equation: where V c is crater volume, t on is pulse-on time and t off is the pulse-off time. A program using MATLAB 2018 (MathWorks Inc., Natick, MA, USA), was developed to calculate MRR based on Equations (18) and (19).
Machines 2019, 7, x www.mdpi.com/journal/machines = + (19) where Vc is crater volume, ton is pulse-on time and toff is the pulse-off time. A program using MATLAB 2018 (MathWorks Inc., Natick, MA, USA), was developed to calculate MRR based on Equations (18) and (19).   The dependency of MRR on the pulse duration for different currents is shown in Figure 3. For all the discharge currents, MRR was found to increase with pulse duration to a maximum value and then decrease again. The reason for this trend might be the generation of lower heat densities at a higher pulse duration, and hence, MRR decreased. In addition, higher values of discharge current were also related to the higher MRR. Figure 4 illustrates the effect of the discharge current on MRR for different discharge voltages. Obviously, any increase in current or voltage led to a higher heat flux, which in return caused melting and the removal of a higher amount of material (Equation (3)). The dependency of MRR on the pulse duration for different currents is shown in Figure 3. For all the discharge currents, MRR was found to increase with pulse duration to a maximum value and then decrease again. The reason for this trend might be the generation of lower heat densities at a higher pulse duration, and hence, MRR decreased. In addition, higher values of discharge current were also related to the higher MRR. Figure 4 illustrates the effect of the discharge current on MRR for different discharge voltages. Obviously, any increase in current or voltage led to a higher heat flux, which in return caused melting and the removal of a higher amount of material (Equation (3)).

Effect of Process Parameters on Surface Roughness
As explained earlier, for the surface roughness evaluation, several discharges on the surface along the x-axis were considered. Figure 5 shows the temperature distribution after the fourth discharge. This figure shows a higher bulk temperature of the workpiece. In addition, a slight increase in the crater depth after each discharge can be seen. Figure 6 shows the depth of the first four craters.

Effect of Process Parameters on Surface Roughness
As explained earlier, for the surface roughness evaluation, several discharges on the surface along the x-axis were considered. Figure 5 shows the temperature distribution after the fourth discharge. This figure shows a higher bulk temperature of the workpiece. In addition, a slight increase in the crater depth after each discharge can be seen. Figure 6 shows the depth of the first four craters. In all cases, after four discharges, the crater depth converged to a constant value. Henceforth, the crater was used to calculate the surface roughness.   According to Figure 7, the surface roughness increased with the increase in both discharge time and discharge current. This can be explained with the increase in discharge time, where more energy was transferred to the workpiece, hence resulting in a bigger crater. In addition, due to the lower According to Figure 7, the surface roughness increased with the increase in both discharge time and discharge current. This can be explained with the increase in discharge time, where more energy was transferred to the workpiece, hence resulting in a bigger crater. In addition, due to the lower energy densities at a higher discharge time, the dimensions and volume of the crater grew slowly. As a result, the rate of increase in the surface roughness decreased. This behavior is in agreement with experimental work by Kiyak and Cakir [41] and Keskin et al. [42].
Machines 2019, 7, x 12 of 18 Machines 2019, 7, x www.mdpi.com/journal/machines energy densities at a higher discharge time, the dimensions and volume of the crater grew slowly. As a result, the rate of increase in the surface roughness decreased. This behavior is in agreement with experimental work by Kiyak and Cakir [41] and Keskin et al. [42].  Figure 8 presents the effects of current on surface roughness for different discharge voltages. As explained earlier, an increase in either discharge voltage or current caused more energy transfer to the workpiece. Therefore, larger craters and consequently higher surface roughness occurred. A similar trend is reported by Kiyak and Cakir [41].   Figure 8 presents the effects of current on surface roughness for different discharge voltages. As explained earlier, an increase in either discharge voltage or current caused more energy transfer to the workpiece. Therefore, larger craters and consequently higher surface roughness occurred. A similar trend is reported by Kiyak and Cakir [41].   Figure 8 presents the effects of current on surface roughness for different discharge voltages. As explained earlier, an increase in either discharge voltage or current caused more energy transfer to the workpiece. Therefore, larger craters and consequently higher surface roughness occurred. A similar trend is reported by Kiyak and Cakir [41].

Model Validation
In order to validate the model, both MRR and surface roughness were compared against the previously published results. Figure 9 shows the material removal rate obtained from our model with and without considering instantaneous evaporation, as well as numerical analysis by Joshi [21] and experimental results reported by DiBitonto [7]. Table 1 provides detailed information about the experimental input parameters used for comparison. As it can be seen, there is good agreement between our model and experimental results, thereby validating our numerical scheme. Compared to the standard FEM model, the FEM model with instantaneous evaporation had 1.5% better agreement with the experimental values. Figure 10 shows a comparison of the model with numerical and experimental results from Almacihna et al. [22]. For discharge energies of 420, 660, and 1340 mJ, Almacinha provided two sets of data; lower values of MRR correspond to experiments without flushing while higher values were experiments when flushing was activated. For the first two points, the model shows very good agreement with the experiments when flushing was performed. This was due to the fact that in the model, no recast layer was considered and all molten material was removed, a condition that is more likely to happen when fluid flushing is active in the EDM process. The discrepancy at high discharge energies is possibly due to the invalidity of the heat input radius for energies higher than 670 mJ.
In Figure 11, the surface roughness from our model, and theoretical analysis of Salonitis [29] and experimental results are shown. Table 2 provides the detailed information about the experimental input parameters used for comparison. In Salonitis's work, the surface roughness was calculated from the average heights of two overlapping craters. Parameters in their model were evaluated based on a crater generated by a single discharge. According to Figure 11, the surface roughness predicted by our model was closer to the experimental data than for Salonitis' model. A possible explanation for the discrepancy between Salonitis and our model is that Salonitis' model considered many simplifying assumptions such as uniform heat source. The deviation between our model and Salonitis model and experimental values were 6.4% and 7.5%, respectively. Another comparison of surface roughness obtained from the model with experimental values [22] are shown in Table 3. The model predicts the increasing trend of surface roughness with an increase in discharge energy. There is no exact relation between mean surface roughness (Ra) and ten-point mean roughness (Rz). According to Whitehouse [43] and DIN 47 [44], ten-point roughness (Rz) is between 3 to 10 times larger than mean surface roughness (Ra). Via consideration of the mentioned relation between Ra and Rz, numerical results (Ra) obtained from the model lie between these two bounds (i.e., 1/3 and 1/10 times of Rz), therefore validating our model.          Figure 11. Comparison of surface roughness obtained from a numerical model and experimental values from Salonitis [29].

Conclusions
Numerical models based on the heat conduction equation for the evaluation of crater geometry, MRR, and surface roughness of the EDM process were proposed. It can be concluded that more realistic results can be obtained by considering the effects of the residual heat generated from each spark and instantaneous evaporation. Additionally, the dependency of both MRR and surface roughness on process parameters were determined. Good agreement between our model and experimental data have been illustrated. Comparisons of models with and without instantaneous

Conclusions
Numerical models based on the heat conduction equation for the evaluation of crater geometry, MRR, and surface roughness of the EDM process were proposed. It can be concluded that more realistic results can be obtained by considering the effects of the residual heat generated from each spark and instantaneous evaporation. Additionally, the dependency of both MRR and surface roughness on process parameters were determined. Good agreement between our model and experimental data have been illustrated. Comparisons of models with and without instantaneous evaporation demonstrates that for a shorter discharge time and higher energies, evaporation plays a significant part in the removal of material. The model with IE had 1.5% better prediction performance regarding MRR. Also, the residual heat in each spark could not be ignored in the evaluation of the surface roughness. Therefore, authors recommend evaluating MRR and surface roughness via a more comprehensive 3D thermal model. In addition, if the effect of dielectric flow on surface roughness can be incorporated, better correlation with empirical data can be expected. In recent years, ceramics have gained more attention in industries. Despite this fact, theoretical modeling of EDM for ceramic materials is still scarce. Future work will also be focused on the numerical investigation of EDM for non-conductive ceramic materials.
Author Contributions: A.P. conceived the idea of this research; A.R. conducted the modeling work and wrote part of the paper under the supervision of C.M. and A.P.; D.T. assisted in editing the paper.
Funding: This research study was funded by Nazarbayev University under the project "Multi-scale Investigation of the Machining Behavior of Non-Conductive Ceramics Using Electro-Discharge Machining" (grant No. 090118FD5324). The APC was funded by 090118FD5324.

Acknowledgments:
The authors would like to express their deepest gratitude to Pauline Mcloone for checking the English of this paper.

Conflicts of Interest:
The authors declare no conflict of interest. Laplace operator Q max maximum heat flux σ standard deviation of distribution