Electro-thermal-based finite element simulation and experimental validation of material removal in static gap single-spark die-sinking electro-discharge machining process

Abstract

Existing single spark models are subjected to too simplistic assumptions such as uniform or point heat source, constant plasma radius, invariable materials properties and constant surface temperature during discharge making them far from reality. In this study, more realistic assumptions including Gaussian type distribution of spark heat flux, temperature dependent materials properties, latent heat of melting and expanding plasma channel with pulse current and time have been made to establish a comprehensive modeling platform. The ABAQUS FEM software has been used to simulate the mechanism of crater formation due to a single discharge. The non-uniform thermal flux was programmed through the DFLUX subroutine. The simulation results show that the temperature of work piece decreases as the discharge time increases while the volume of melted and evaporated material increases. A specially designed single spark experimental set-up was developed in laboratory to carry out a few single spark tests for verification purposes. The obtained craters morphologies were examined by optical microscopy and scanning profilometer. It has been shown that the present approach outperforms other previously developed thermal models with respect to cavity outline and size possessing the maximum confirmation errors of 18.1% and 14.1% in predicting crater radius and depth, respectively. Parametric analysis reveals that the melting boundary moves onward by increasing discharge current, whereas it moves back prolonging discharge time. Finally, a closer proximity to experimental material removal rates than those predicted by analytical approach has been recognized which confirms its more precise generalization capabilities towards the real state EDM process.

Keywords

Electro-discharge machining finite element method material removal rate simulation single discharge temperature distribution

Introduction

Electro-discharge machining (EDM) is indisputably an electro-thermal process that removes material from the work surface through the action of repetitive discrete electrical sparks produced by a direct current (DC) pulse generator between the two electrodes (tool and workpiece) separated by a small distance (gap) while being immersed in a dielectric liquid medium. The thermal energy is then utilized to generate high-temperature plasma channel that erodes minute amounts of the work material through local melting and vaporization. The removed material particulates (debris) are subsequently flushed off from the surface of workpiece via cool dielectric flowing into machining region. Since the EDM does not involve any physical contact between tool and workpiece, mechanical factors such as hardness or strength are no barrier to its applications. It thus provides an economic advantage for making sampling tools, wire drawing and extrusion dies, forging dies, hardened steels, intricate mold cavities and so on. Delicate jobs, being too fragile to withstand cutting forces during conventional machining processes, can also be effectively machined by the EDM.^1–3

Despite the widespread applications of EDM in today’s industry, such as die and mold, aerospace and automobile, the process itself is inherently complex and stochastic in nature as several disciplines of science and branches of engineering are involved. Phenomena such as heat conduction, energy distribution, melting, vaporization, ionization, formation and collapse of gas bubbles all take place in a very short period of time (of the order of microsecond). Therefore, it is too difficult to develop a universally accepted mathematical model predicting its performance measures while entirely covering the involved physical aspects.⁴ Poor quantitative knowledge of the sources of variability affecting this process hinders the identification of its natural tolerance limits. Until now, there is no acceptably comprehensive physics-based process model of the EDM due to the presence of a large number of interconnected variables and complicated process mechanism.⁴

Among several theories suggested for the mechanism of material removal in EDM, the electro-thermal approach has gained further popularity being more consistent with experimental observations.⁴ It proposes intense local temperature gradients associated with a very high energy density caused by ionization within the discharge column formed between the workpiece and the tool. Hence, except for very short pulse on-times (less than 5 µs) in which there is not enough time for the material to be effectively heated and the induced electrostatic forces have been recognized to be responsible for material removal,⁵ the thermal effects resulting in melting and vaporization have been known as the main mechanisms of material removal in EDM process.⁶

Therefore, since EDM’s inception and especially in the early 1970s, many researchers attempted to develop thermal models describing the discharge phenomena (plasma channel) and process physics to model the mechanisms of anode and cathode erosion.⁷ These early researches have all been based on obtaining analytically exact closed-form solutions through solving heat conduction problem applying different boundary conditions and other simplifying assumptions such as cylindrical spark plasma,⁸ uniform (disk) heat source,⁹ point heat source,⁶ expansion of plasma channel radius with only pulse on-time while ignoring the effect of discharge current¹⁰ and constant thermo-physical properties of materials¹¹ taken over the temperature domain concerned in EDM. However, recent experimental observations made through high-speed video cameras, optical spectroscopic analysis and fiber optics^12,13 have shown that the spark plasma grows with time and is of barrel shape, the cathode arc root being much smaller than the anode arc root due to electron emission while both being smaller than the bulk plasma. It has also been indicated that the distribution of heat flux within the discharge channel is not uniform but changing with time and position over the area where the spark develops. Moreover, the assumption of constant material properties is not in line with the true physical aspects of the process as it involves extremely high temperatures exceeding the work melting point. Accordingly, it is necessary to apply temperature-dependent thermo-physical properties to simulate the process mechanism more accurately. To sum up, all the above referred works can be criticized in some parts considering such simplifying assumptions which severely limit their applicability and hence require suitable modifications.

Salah et al.¹⁴ applied finite difference schema in two-dimensional (2D) case to study temperature distribution in EDM process and obtained material removal rate (MRR) resulted from single discharge. It was indicated that the use of temperature-dependent thermal conductivity is of crucial importance enhancing modeling accuracy and fitting experimental results better than the case when a constant conductivity is assumed.

Shankar et al.¹⁵ conducted a comprehensive attempt based on finite element method (FEM) to analyze the spark profiles and reported that the plasma shape is non-cylindrical having the smallest cross-sectional area at the middle of discharge channel. However, their results are in conflict with those experimental evidences referenced by advanced technological setups^12,13 that revealed a non-cylindrical barrel shape for a developed plasma channel.

Ahn and Chung¹⁶ performed a numerical analysis of unsteady approach on alumina–titanium carbide composite under EDM conditions with variable discharge currents and duty cycles. It was found that increasing both current intensity and duty cycle results in higher amounts of erosion rates.

Marafona and Chousal¹⁷ developed an electric–thermal model considering the effect of Joule heating as the main source of thermal energy increasing discharge channel temperature using the ABAQUS/Standard FEM software. Comparing their results with those experimental values reported by AGIE SIT Corporation, USA, for anode and cathode erosion rates,^6,11 a good agreement was claimed. Nevertheless, the effect of latent heat of melting was ignored and constant material properties were used as well.

More recently, Somashekhar et al.¹⁸ applied FEM based on the electro-thermal theory to approximate the dimensions of a typical crater formed in micro-EDM. Large deviation between experimental and simulated result was noted, attributed to the assumption of a constant spark radius (2 µm) during discharging.

Although studying the consequences of occurring a single-pulse discharge over a work surface may not completely reveal the whole nature of EDM process in its real state multi-spark phenomenon, it is of great importance helping understand the fundamental mechanism of metal erosion by spark discharge. Additionally, it can be used as a useful tool in parametric analysis investigating the effects of different parameters on the process measures such as crater size, geometry and MRR. This provides a versatile mean to predict process behavior under different operating conditions yielding a cost-effective and economical approach eliminating the need of running exhaustive experimentations. In short, based on the review of related literature considering their restrictions and limitations, the main objectives of the current research can be stated as follows:

Shed a fundamental insight on the principle working mechanism of EDM process using more realistic assumptions being in line with the true nature of EDM. These mainly involve the non-uniform Gaussian distribution of thermal plasma, growth of spark radius, variable material properties, latent heat of fusion and the fraction of discharge energy entering the workpiece.

Study the actual crater morphologies obtained through conducting a few single-spark tests using a special single-discharge experimental setup and verify the simulation results. This also provides a unique basis to evaluate and compare the predictive accuracy of the present model with other approaches against the experimental findings.

Methodology

Electro-thermal modeling of single-spark EDM process

The discharge occurrence in EDM can be modeled as heating the work surface by plasma formed during the process. As cited before, previous researchers have mainly used disk heat source with uniformly distributed heat flux⁹ or point heat source model⁶ for thermal modeling of EDM process. However, these assumptions are simplistic and not in line with the real EDM process. In fact, recent experimental observations through spectroscopic techniques done by Kojima et al.¹² and Descoeudres et al.¹³ have revealed that the temperature distribution within the plasma channel is non-uniform but rather very similar to a Gaussian sharing in which the temperature is maximum at the center of plasma channel. Moreover, the real shape of a crater formed by a single spark can best be approximated as Gaussian distributed.¹⁹ Hence, it is imperative to apply a heat flux with Gaussian distribution as it is more reliable and realistic to characterize the heat flux inside the EDM plasma in terms of relevant process parameters controlling the discharge energy.

Figure 1 illustrates the ideal case in which a single spark has been modeled as a heat source with Gaussian-distributed heat flux with appropriate boundary conditions considering the axisymmetric nature of domain of the analysis. The next sections will elaborate more on the details needed to implement the FEM simulation of a single discharge.

Figure 1.

Thermal model of single-spark EDM process: (a) schematic representation of workpiece domain and (b) axisymmetric model with designated conditions.

Modeling assumptions

As the EDM is a highly complex and stochastic process, several suitable assumptions are needed to make its modeling approach mathematically tractable. They are as follows:

As the volume of material removed per each discharge is so small compared to workpiece and tool electrodes’ volumes, both the workpiece and tool are considered as semi-infinite bodies.

The investigated domain for thermal analysis is assumed to be axisymmetric with respect to r–z plane (see Figure 1).

Only one spark happens during each pulse, that is, every pulse duration (on-time + off-time) is equivalent to a single spark.^6,20

Only a constant fraction of spark thermal energy is transferred to the workpiece, the rest is used for tool wear as well as dissipating via convection and radiation through the dielectric medium.⁶ This assumption will be discussed more in section “Energy sharing in gap” and later on at the end of section “Effects of modeling assumptions on the accuracy of results and comparison with other approaches.”

Conduction is regarded as the main mode of heat transfer from the spark channel to the workpiece surface^6,9,14 as the contribution of other mechanisms such as convection and radiation is negligible.

The workpiece is assumed to be homogeneous and isotropic while its thermo-physical properties such thermal conductivity, specific heat and density are considered temperature dependent due to the high temperature gradients created during discharging.

All the materials whose temperature exceeds the melting temperature are assumed to be effectively removed by the dielectric rush upon the collapse of plasma channel. In other words, there is no deposition of recast layer at the bottom of crater cavity formed at the end of each discharge.

Governing equation

It has been established that the main mode of heat transfer from the plasma channel to the work surface in EDM process is by conduction. Therefore, the principal equation, the Fourier transient heat conduction equation with no internal heat generation, which has to be solved, can be stated as follows

\frac{1}{r} \frac{\partial}{\partial r} (k_{t} r \frac{\partial T}{\partial r}) + \frac{\partial}{\partial z} (k_{t} \frac{\partial T}{\partial z}) = ρ C_{p} \frac{\partial T}{\partial t}

(1)

where T is the temperature; r and z are the radial and axial coordinates, respectively; t is the time; ρ is the density; $C_{p}$ is the specific heat and $k_{t}$ is the thermal conductivity of workpiece material.

Input heat

Three main issues contribute to the formation of input heat in EDM, namely, plasma channel radius, heat flux function and energy share within gap during sparking. In order to provide suitable boundary conditions for the incident plasma on the work surface, they should be addressed properly and are explained in the following subsections.

Selection of spark radius growth function

The discharge channel radius is not constant, expands with time during discharging and depends on several factors such as spark energy, tool and workpiece materials as well as dielectric medium as was proved by Erden²¹ and Natsu et al.²² The most comprehensive relation proposed so far which is also used in this research is the one attributed to Ikai and Hashiguchi,²³ known as the semi-empirical equivalent heat source radius written as

R_{pc} (μ m) = 2.04 \times I^{0.43} \times t^{0.44}

(2)

where $I$ is the discharge current in Amperes and $t$ the time in microsecond which varies between 0 and the end of pulse on-time ( $[0, T_{on}]$ ). It should be noted that this equation is valid for spark energy levels up to 670 mJ.²³ So, care has been paid so as not to violate this condition by properly selecting the domain of variation in input parameter settings controlling the discharge energy in both the simulation and experimental steps.

Selection of heat flux function

The probability density function of Gaussian distribution for a random variable ( $p (r)$ ), regarded as the intensity or plasma heat flux ( $q_{w} (r)$ ) imparting to the workpiece, can be written as²⁴

q_{w} (r) = p (r) = \frac{1}{σ \sqrt{2 π}} e^{- \frac{r^{2}}{2 σ^{2}}}

(3)

where $σ$ is the standard deviation and $1 / σ \sqrt{2 π}$ the maximum value of distribution. As this function converges to 0 while the distance ( $r$ ) tends to infinity, hence, it is necessary to assign some large enough values for its argument to obtain a representative criterion for the bottom of crater formed due to a single pulse. Usually, a range of six sigma ( $- 3 σ to + 3 σ$ ) is considered to evaluate the domain in which the response drops to 0.25% of its initial value (99.75% of the values lie between $r = μ - 3 σ$ and $r = μ + 3 σ$ ). Therefore, $R = 3 σ$ or $σ = R / 3$ , upon substitution into equation (3)

p (r) = \frac{3}{R \sqrt{2 π}} e^{- 4.5 \frac{r^{2}}{R^{2}}}

(4)

Assuming $q_{0} = 3 / R \sqrt{2 π}$ as the maximum heat flux at the center of spark ( $r = 0$ ), the above equation can be written as

q_{w} (r) = q_{0} e^{- 4.5 \frac{r^{2}}{R^{2}}}

(5)

Hence, the total amount of thermal power supplied by the spark over the work surface is

\begin{matrix} \oint q_{w} (r) dA = \int_{0}^{R} q_{w} (r) 2 π rdr = \int_{0}^{R} q_{0} e^{- 4.5 \frac{r^{2}}{R^{2}}} 2 π rdr \\ = - \frac{π R^{2}}{4.5} q_{0} e^{- 4.5 \frac{r^{2}}{R^{2}}} | \begin{matrix} R \\ 0 \end{matrix} \\ = \frac{π R^{2} q_{0}}{4.5} (1 - e^{- 4.5}) = 0.2191 π q_{0} R^{2} \end{matrix}

(6)

On the other hand, the effective amount of heat rate entering the workpiece can be identified to be a fraction of the overall electrical power, that is, $F_{c} VI$ , which is converted to thermal energy and absorbed by the cathode workpiece for erosion purposes. So

F_{c} VI = 0.2191 π q_{0} R^{2}

(7)

where $F_{c}$ stands for the cathode (workpiece) fraction of discharge power. Solving for $q_{0}$ and substituting into equation (4) yields

q_{w} (r) = \frac{4.57 F_{c} VI}{π R^{2}} e^{- 4.5 {(\frac{r}{R})}^{2}}

(8)

Since the plasma channel radius (R) expands with time till collapsing at the end of pulse on-time (see equation (2)), it can be inferred that the thermal flux in EDM depends on both discharge time (t) and radial position (r).

Energy sharing in gap

The real effective fraction of discharge energy entering each electrode (tool or workpiece) depends on a number of factors, such as thermal properties, flushing conditions, gap distance, conductivity of dielectric, polarity, discharge current and pulse duration.^25–27 It is extremely difficult if not possible to define and identify a unique physical quantity fully reflecting the material properties and EDM conditions altogether to be used in assigning the energy share of each electrode.

In this regard, the literature presents a large discrepancy in selecting a reference value for the fraction of energy dispersing into the cathode workpiece (F_c ). It has generally been chosen as a constant percentage of the total discharge energy ranging from 14% to 50% by different researchers.²⁸ Yeo et al.²⁹ conducted a comparative study between five analytically based electro-thermal models of EDM suggested by earlier researchers and concluded that the predictive accuracy of MRR models can be enhanced should a proper value be assigned for F_c . Joshi and Pande³⁰ scrutinized further into the effects of various amounts of F_c on erosion rates of steel workpiece and concluded that this fraction does not vary substantially using copper tool in kerosene dielectric. On the other hand, using water as dielectric, Shankar et al.¹⁵ suggested that about 40%−50% of thermal energy can be received by the 5 Cr die steel as the workpiece.

As can be inferred, there is not a universally accepted value for F_c as these reported amounts have all been obtained under specific circumstances with limited applicability. In this research, F_c has been set constant at 18% for simulation purposes. This has been experimentally verified to be a reliable quantity by DiBitonto et al.,⁶ which later was also reconfirmed to be used as a proper value of F_c while EDM of steel using copper tool in kerosene dielectric medium.³⁰ The validity of such a selection will be examined and discussed more at the end of section “Effects of modeling assumptions on the accuracy of results and comparison with other approaches.”

Boundary and initial conditions

According to Figure 1(b), the boundary conditions over the border B ₁ can be written as

- k \frac{\partial T (r, 0, t)}{\partial z} = {\begin{matrix} q_{w}, r \leq R \\ h_{c} (T - T_{0}), r > R \\ 0, off time \end{matrix} on B_{1}

(9)

where q_w is the thermal flux generated by a single discharge expressed by equation (8), h_c the convective heat transfer coefficient due to the heat exchange between the upper face of workpiece and dielectric, T ₀ the temperature of dielectric, k the thermal conductivity of workpiece and R the plasma channel radius (equation (2)). The boundaries B ₂ and B ₃ are supposed to be located far enough from the discharge location and can then be specified as insulated borders while the boundary B ₄ acts as a symmetry line with no heat flux passing through it. Hence, $\partial T / \partial n = 0 : on B_{2}, B_{3}, B_{4}$ , where n is a unit vector normal to the boundary, and $T_{i} = T_{0} = 25^{\circ} C$ is the initial condition for all the domain at t = 0.

Latent heat of fusion due to phase change

In EDM, material removal takes place mainly by melting and vaporization phenomena. However, the contribution of vaporization is negligible compared to melting as a result of tremendously large dynamic pressures within the plasma channel which further restrict material evaporation.^31,32 Therefore, to enhance the simulation accuracy, it is essential to include the latent heat of fusion which has been neglected by many researchers.^6,16–18 To that end, the specific heat (C_P ) is amended to yield the effective specific heat (C_Peff ), using the following relation³³

C_{Peff} = C_{P} + \frac{L_{H}}{Δ T}

(10)

where L_H is the workpiece latent heat of melting and ΔT = T_m − T ₀ the temperature difference between melting point temperature and dielectric temperature (for the AISI 304 steel: L_H = 247 kJ/kg to be released or absorbed over the temperature range of solidus temperature T_s = 1477 °C to liquidus temperature T_m = 1535 °C).³⁴

Finite element formulation and solution technique

Applying the variational method,³⁵ the principle governing equation (equation (1)) can be written in final assembled finite element (FE) form as

[C_{pe}] {{\overset{\cdot}{T}}_{e}} + ([K_{e}^{d}]) {T_{e}} = {q_{e}} + {q_{e}^{c}}

(11)

where $[C_{pe}] = \int_{Ω} ρ C_{p} {N} {N}^{T} d Ω$ is the elemental capacitance or specific heat matrix; $[K_{e}^{d}] = \int_{Ω} {[B]}^{T} [K] [B] d Ω$ the elemental diffusion conductivity matrix; ${q_{e}} = \int_{Γ} {N} q_{w} d Γ$ the elemental heat flux vector; ${q_{e}^{c}} = \int_{Γ_{1}} h (T_{0} {N} - {N} {N}^{T} {T_{e}}) d Γ$ the elemental convection surface heat flow vector; $[B]$ the shape function derivative matrix; $[K]$ the conductivity matrix; $[N]$ the element shape function vector; ${T_{e}}$ the element’s nodal temperature vector ( $T = {N}^{T} {T_{e}}$ ) and $Ω$ , $Γ$ and $Γ_{1}$ represent the domain of analysis, the domain of boundary element where heat flux spreads till the end of discharge and the domain of boundary element where convective heat transfer occurs, respectively (see Figure 1). Finally, in differential form, equation (11) can be written as

[C_{G}] {{\overset{\cdot}{T}}_{G}} + [K_{G}] [T_{G}] = {Q_{G}}

(12)

where $[C_{G}]$ is the global capacitance or specific heat matrix, $[K_{G}]$ the global conductivity matrix, ${Q_{G}}$ the global heat flux vector, ${T_{G}}$ the global temperature vector and ${{\overset{\cdot}{T}}_{G}}$ the time derivative of ${T_{G}}$ .

The FEM has an important feature over the analytical approaches, and that the temperature dependency of work material properties can be incorporated into the calculations. Table 1 shows the variable thermo-physical properties of the AISI 304 stainless steel as the material used in both simulation and experimental works, while Table 2 lists the variations in convective heat transfer coefficient with temperature over the top surface of work surface (see Figure 1).

Table 1.

Temperature-dependent thermo-physical properties of AISI 304 stainless steel used in simulations.³⁴

Temperature T (°C)	Conductivity K (W/m K)	Specific heat C_p (J/kg K)	Density ρ (kg/m³)	Temperature T (°C)	Conductivity K (W/m K)	Specific heat C_p (J/kg K)	Density ρ (kg/m³)
0	51.9	450	7872	600	35.6	773	7669
75	51.3	486	7852	675	32.8	846	7636
100	51.1	494	7845	700	31.8	1139	7625
175	49.5	519	7824	730	30.1	1384	7612
200	49	526	7816	750	28.9	1191	7602
225	48.3	532	7809	775	27.5	950	7590
275	46.8	557	7763	800	26	931	7578
300	46.1	566	7740	1000	27.2	779	7552
325	45.3	574	7717	1500	29.7	400	7268
375	43.6	599	7727	1540	29.7	400	7218
400	42.7	615	7733	1690	29.7	847	7055
475	40.2	662	7720	1840	29.7	847	6757
500	39.4	684	7711	1890	29.7	400	6715
575	36.6	749	7680	2860	29.7	400	5902

Table 2.

Temperature dependency of convection heat transfer coefficient (h) for kerosene dielectric medium.³⁶

Temperature T (°C)	Convection heattransfer coefficient h (W/m² K)
20	688.9
200	688.9
600	317.7
720	836.2

Mesh, FEM parameter settings and model continuum geometry

The workpiece is modeled using a 2D continuum with a dimension of 600 µm × 600 µm. In order to capture the extreme temperature gradients arising due to a single-pulse discharge, very small-sized elements should be assigned at the vicinity of discharge location. Hence, an area of dimension 200 µm × 200 µm with a much denser mesh quality is created within the entire domain of analysis, whereas beyond that the mesh density gradually decreases conveying a non-uniformly distributed FE mesh with much more elements mapped toward the heat-affected zone (Figure 2).

Figure 2.

2D mesh of a quarter of workpiece used in simulations.

To code a non-uniform Gaussian distribution of heat flux with regard to time and distance imposed over the top surface of the workpiece, the user-defined subroutine DFLUX has been programmed in FORTRAN language and called during the simulation using ABAQUS FEM software, implicit (standard) module. It works by placing the local origin at the center of plasma channel (top left corner of workpiece) and calculating the value of heat flux at each time increment considering the radial distance of each node surrounding this center.

The boundary (B ₂) should be far enough from the axis of symmetry (spark center) to be assumed as an adiabatic surface. Byron Bird et al.³⁷ introduced the parameter “thermal penetration depth ( $δ_{D}$ )” as an extreme distance beyond which the region can be considered insulated without any heat exchange in conduction heat transfer analysis. The thermal penetration depth can be written as

δ_{D} = 4 \sqrt{α t}

(13)

where $α$ is the thermal diffusivity ( $K / ρ C_{p}$ ) and $t$ the discharge time. For a pulse on-time of 600 µs as maximum discharge duration, this value would be about 200 µm. Hence, the borders B ₂ and B ₃ being 600 µm far away from the spark center (Figure 1) may confidently be supposed as adiabatic boundaries.

Table 3 summarizes the most important parameters used for the FEM thermal analysis of a typical single-spark EDM process. To select the right element size, a number of simulations were done using elements of different sizes. 0.5, 1, 2 and 3 µm square-type elements were tried and the nodal temperatures were obtained and recorded. Figure 3 shows the comparison of temperature distribution profiles along a typical radial path of 10 µm length on top of the workpiece surface. As is clear, subtle discrepancies can be observed among these cases, especially when the element size is taken to be 1 or 2 µm. As an instance, the peak temperatures at the plasma center are 25362.57 °C, 25362.62 °C, 25362.64 °C and 25358.99 °C for the four aforementioned element sizes, respectively. Similar results have also been obtained along other paths. It is to be mentioned that all the simulations were done on an Intel Pentium IV of 2 Quad Core CPU with 3.23 GB RAM and the times taken to fulfill the tasks were about 83, 20, 6 and 2 min, respectively. Therefore, in view of getting a reliably stable solution within a suitable time frame needed to run the simulation procedure, the optimum element size can be considered to be 2 µm providing a higher computational efficiency. In this case, the considered area will have a total number of 20,201 elements (Figure 2).

Table 3.

Parameters used for a typical single-discharge FEM simulation.

Type of analysis	Two-dimensional transient thermal
Element used	DCAX4 (four-node linear axisymmetric heat transfer quadrilateral)
Solution method	Iterative Newton–Raphson
Time step	Fixed
FEM domain size	600 × 600 µm²
Meshing type	Structured meshing
Discharge current (A)	20
Spark on-time (µs)	100
Discharge voltage (V)	25

FEM: finite element method.

Figure 3.

Comparison of simulation results with different element sizes (I = 20 A, T_dis = 100 µs and V_dis = 25 V).

Experimental details

The literature reveals very limited experimental studies on the shape and size of a real crater formed as a result of single discharge. As there is no commercially available apparatus to produce single sparks, an in-house single-spark experimental setup which was specifically designed and built in the laboratory was used to verify the FEM simulations of crater formation under the same conditions. It produces a unique discharge each time by creating a single pulse, the characteristics of which could be selected on the machine manually. Figure 4(a) and (b) illustrates a sample of tool/workpiece used in single-spark experiments and the experimental rig with a close view of discharge cell and different enumerated components, respectively. Side flushing was applied using kerosene as the dielectric liquid while the work specimen was submerged within the tank. A minimum gap distance needs to be attuned so as to be able providing suitable conditions leading to the occurrence of a single-spark over the work surface. A manual adjusting wheel driving the head and positioning the tip of cone-shaped tool electrode over the work surface has been used to control and maintain a constant distance between the tool and the workpiece as the spark gap along the Z-axis. Once the tool comes into contact with the work surface, an alarm buzzer alerts the operator that a zero gap width has been reached. The gap span can then be adjusted and controlled by revolving the manual wheel in a reverse direction with the help of an installed dial indicator (micrometer) having 0.01-mm accuracy.

Figure 4.

(a) Sample of tool/workpiece. (b) Single-spark experimental setup with a close view of discharging cell.

Before running tests and to better observe and measure the dimensions of cavities formed due to single sparks, the workpiece surfaces were first ground and then polished by a very fine emery sheet till a mirror-like finish was obtained. The discharge waveforms of EDM, that is, discharge pulse voltage and discharge pulse current, reflect the amount of discharge energy released during sparking and it is of vital importance to accurately measure and record these values to verify the results of FEM simulations. A Tektronix two-channel digital storage oscilloscope, model TDS 1002, 60 MHz sampling rate, was employed and suitably connected to the experimental setup (Figure 4(b)) to capture the discharge waveforms and quantify the relevant amounts of discharge current and gap voltage. Figure 5 shows a sample of the discharge waveforms of the first single-spark trial according to the parameters: V_dis = 20 V, I_dis = 40 A and T_dis = 20 µs, obtained by the oscilloscope during sparking. Table 4 summarizes the conditions that could normally be set on the machine to carry out single-spark experiments. Several replications were made for each parametric combination to ensure achieving reliable data.

Figure 5.

Typical discharge waveforms captured by oscilloscope and characterized for spark energy: (a) voltage pulse and (b) current pulse (source voltage V_s = 250 V, discharge voltage V = 20 V, discharge current I_dis = 40 A, pulse on-time T_on = 22 µs, discharge delay time t_d = 2 µs, and discharge duration T= 20 µ; where these values have been obtained based on the scales set for each unit on the relevant coordinate axes; 50 V and 10 µs for the vertical and horizontal units, respectively).

Table 4.

Experimental single-spark EDM conditions.

Tool material	Copper
Work material	AISI 304 stainless steel
Dielectric	Kerosene
Gap distance (mm)	0.01
Source voltage (V)	250
Discharge voltage (V)	20
Discharge current (A)	40, 50, 60
Pulse on-time (µs)	22, 125, 200

Optical metallographic microscope, METALLUX 3, manufactured by LEITZ Company, Germany, equipped with Image Analyzer software has initially been used to observe crater morphologies at 50× magnification and approximately measure the corresponding radius of each crater cavity formed. Figure 6 illustrates the microscopic image of the respective crater formed. Preliminary evaluation through optical microscope estimates the average crater radius as 48.2 µm. As can be seen, a somehow circular-type crater with a portion of molten metal condensed around its edge is evident. The is due to the abrupt collapse of discharge power at the very end of pulse on-time as well as the cooling effect of cold dielectric rushing into the discharge spot, which, in turn, results in the distribution and protrusion of molten material around the crater rim.³⁸

Figure 6.

Microscopic image of the crater formed (I_dis = 40 A, T_dis = 20 µs and V = 20 V).

In addition, each crater has been examined by Form Talysurf Series 2, a three-dimensional (3D) topography measurement instrument, from Taylor Hobson Company, to better analyze and reveal the crater surface morphology. Figure 7(a) shows the 3D image of the crater cavity produced as per the discharge conditions illustrated in Figure 5. The shape of the cavity has been examined at different sections and Figure 7(b) depicts a typical cross-sectional view chosen for further analysis. As it is clear, the shape of the cavity is neither hemispherical nor disk-like as was assumed by previous researchers,^6,39–41 but rather bowl shape in appearance. Crater diameter and depth have been measured as the major geometrical features. The diameter is 95.2 µm (i.e. radius of 47.6 µm) measured as the pick-to-pick distance and maximum depth of 27.6 µm.

Figure 7.

Surface morphology of a typical single-spark depression: (a) 3D view and (b) crater shape and dimensions (I_dis = 40 A, T_dis = 20 µs and V_dis = 20 V).

As an example, the corresponding FE simulation of a single spark with respect to crater cavity size and shape will be illustrated later for comparison with the current typical experimental finding.

Results and discussions

Evolution of temperature and crater formation

Solution of a heat conduction problem usually indicates the time dependence of temperature at a fixed location in space. An alternative solution, particularly useful for moving boundary problems, is to calculate how specified temperatures move through the heat conducting medium at specified times, using an “isotherm migration method.”⁴² Figure 8(a)–(f) shows the temperature distribution profiles and propagated melting fronts for a representative single-discharge setting as per Table 3, obtained at three distinct conditions, the initial state (t = 5 µs), middle of spark on-time (t = 0.5T_dis = 50 µs) and end of discharge duration (t = T_dis = 100 µs). As can be seen, the isothermal profiles are like quartile concentric circles from the discharge center.

Figure 8.

Temperature distribution contours and propagated molten border at three different discharge steps: (a, b) t = 5 µs, (c, d) t = 0.5 T_dis = 50 µs and (e, f) t = T_dis = 100 µs for I = 20 A and V_dis = 25 V.

From Figure 8(a), it can be observed that the work surface experiences very high temperatures at the early stages of discharge, up to about 8.6 × 10⁴ °C as the maximum temperature which occurs at the center of plasma channel due to the nature of Gaussian flux distribution. This temperature is absolutely high enough to melt and even vaporize the workpiece. Meanwhile, the molten boundary, obtained by eliminating all the elements whose temperatures exceeded 1535 °C (the melting point of AISI 304), spreads only to a small volume of material, 19.5 µm/16.5 µm as the ratio of radius/depth of crater as shown in Figure 8(b). Temperature profiles and the formed molten cavity right at the center of discharge duration have been depicted in Figure 8(c) and (d), respectively. It is interesting to note that the peak temperature now gets a smaller value, decreased to about 3.4 × 10⁴ °C, more than 50% lower than the previous case. The reason can be attributed to the growth of plasma channel with time which results in lower density of discharge energy over the involved area.⁴³ However, the temperature is still large enough to melt the workpiece as the volume of melted puddle extends to a relatively larger volume, radius 45.5 µm × depth 39 µm. The decreasing trend of temperature distribution continues till the end of discharge duration when the highest temperature now gets the lowest amount compared to the two aforementioned cases (about 2.5 × 10⁴ °C, Figure 8(e)). Nonetheless, the induced heat has got adequate time penetrating into the bulk material as the yielded molten cavity reached to its largest value having a size of radius 58 µm × depth 48.5 µm under the simulated conditions (Figure 8(f)). In general, the utmost temperature varies between 2.5 × 10⁴ °C and 8.6 × 10⁴ °C, which is on the analogous order of 6 × 10⁴ K as theoretically predicted by Eubank et al.⁸

Figure 9 illustrates the temperature history of the point locating at the center of plasma channel over the time interval of discharge duration. It can easily be observed that the temperature increases drastically within the first microseconds after the discharge initiation due to the instantaneous application of spark flux, and that there exists an extremely sharp temperature gradient during the early stages of discharge after which a continuous decreasing trend is followed until the end of discharge time.

Figure 9.

Temperature history of the point at the center of plasma channel (I = 20 A, T_dis = 100 µs and V_dis = 25 V).

Figure 10 depicts the temperature distribution profiles alongside both the radial and subsurface directions obtained at two different times, the middle and end of discharge. It can be inferred that the temperature profiles along the axial way possess rather higher gradients than those at radial direction for both simulated conditions. Furthermore, as expected, the feature of Gaussian trend is also reflected by the temperature profiles. As the discharge time passes more, the volume of melted and evaporated material gets higher values while the peak temperature decreases progressively.

Figure 10.

Temperature distribution profiles along radial and subsurface directions at the middle and end of discharge duration (I = 20 A, T_dis = 100 µs and V_dis = 25 V).

The experimental validation of temperature variations and magnitudes with respect to time and position is extremely difficult if not possible due to a very short time interval and size over which a single discharge takes place. However, the final resultant effect of a single spark impinging on the workpiece surface, that is, the formation of crater cavity, can be verified in terms of its size and outline. The results obtained before through running single-spark experiments are used here to verify the simulation results.

Model validation

The FE simulation of crater formation has been done at the same machining conditions reported in section “Experimental details” and the results have been shown in Figure 11(a) and (b) as the molten pool and final resultant cavity, respectively. Only half of the workpiece was depicted as a cross-sectional view due to the axisymmetric assumption of the analysis domain we made before. Elements kill technique has been applied to obtain the border of cavity. It works by removing all the elements whose temperatures exceed the melting temperature of the work material (1535 °C) at the end of discharge making it possible to visualize a hollow area as the crater boundary. Hence, the corresponding values of simulated dimensions are 39 and 31.5 µm for the crater radius and depth, respectively. The discrepancies between the FE-generated and experimental values (Figure 7) are due mainly to the deposition of a portion of melted material over the crater surface known as recast layer. In practice, not all the molten metal can be evacuated at the end of the pulse, which subsequently gets re-solidified. Therefore, the real crater is to some extent shallower than the ideally obtained one. Moreover, the measured crater radius is bigger than the simulated one which is definitely due to the splash of molten material around the crater rim as a result of sudden drop of plasma channel pressure at the end of discharge time and quenching effects of molten material.

Figure 11.

FE simulation of crater formation: (a) molten cavity and (b) crater cavity formed at the end of discharge time (I_dis = 40 A, T_dis = 20 µs and V_dis = 20 V, Experiment 1, Table 5).

Table 5 summarizes the results of verification experiments along with the simulation values for all the three different discharge conditions considered in this research. The last two columns show the relative percentage errors in predicting crater radius and depth. Reasonable modeling accuracy can be observed as the amounts of all errors lie within acceptable range in view of engineering discipline; the highest quantity of which is 18.1% corresponding to Experiment 1 in predicting crater radius. Generally, rational consistency can be noted between the experimental and simulation findings in terms of crater dimensions (radius and depth), hence validating the outcomes of FEM simulations.

Table 5.

Comparison between simulation and experimental results.

Experiment number	Discharge parameters		Crater dimensions (µm)				Relative prediction error (%)
			FE simulation		Experiment		Error in radius	Error in depth
	I (A)	T_on (µs)	Radius	Depth	Radius	Depth
1	40	22	39	31.5	47.6	27.6	18.1	14.1
2	50	125	83.5	64.5	94.2	58.4	11.4	10.5
3	60	200	105.5	79	116.3	70.8	9.3	11.6

V_dis = 20 V. Relative prediction error (%) = (|Experimental value − Simulation value|) / Experimental value × 100.

Effects of modeling assumptions on the accuracy of results and comparison with other approaches

Here, it would be instructive to compare the effect of constant and variable thermo-physical properties of workpiece material and also the inclusion and neglect of latent heat of melting on the accuracy and credibility of simulation results with respect to the experimental findings. For this purpose, while ignoring the role of latent heat of fusion (247 kJ/kg), the average values of thermal conductivity (K_ave : 39.2 W/m K), specific heat (C_Pave : 700 J/kg K) and density (ρ_ave : 7673 kg/m³) over the temperature range from solid to molten state of the AISI 304 steel (see Table 1) have been considered and applied for the FE simulations. Moreover, it has been decided to run the numerical analyses and obtain crater shape and dimensions (radius and depth) using the circular disk heat source model with expanding plasma channel radius^9,11 as well as the point heat source model,⁶ as the other two heat flux functions mostly used in previous researches,⁷ and compare them with the results of this study.

Table 6 illustrates the different values of craters’ radii and depths obtained according to the discharge parameters set for the first verification experiment (Experiment 1, Table 5) along with the amounts of relative verification errors using the aforementioned approaches. First, it is clear that the consideration of variable material properties together with melting heat brings about significant improvement on the accuracy of simulated results. Using a constant set of material properties while excluding the latent heat of fusion will increase the simulation errors from 18.1% to 23.3% and 14.1% to 19.6% in predicting crater radius and depth, respectively. Second, it can be seen that applying either disk heat source model or point heat source to characterize the spark heat flux will also deteriorate the confirmation results (crater shape and size) compared to the Gaussian heat flux proposed in this study. Figure 12((a) and (b)) shows the schematic illustration of a cavity formed when disk and point heat sources were employed. The crater cavity formed due to the disk heat source presents a shallow flat bottom bowl appearance having a radius of almost three times larger than its depth while a semi-spherical cavity was obtained when applying the point heat source model. Neither of which can appropriately approximate the true figure of their experimental counterpart as shown in Figure 7. Hence, the selection of Gaussian distribution of heat flux can be justified here as it can better estimate the real form of crater (see Figures 7 and 11). Moreover, as can be seen from Table 6, when selecting the circular-type heat flux, the error values become 230.9% and 88.4% whereas being 28.2% and 121% for the case of point heat source in predicting crater radius and depth, respectively. These error values are by far larger than those yielded by the present model; therefore, the suitability and logics of our modeling assumptions can be assured.

Table 6.

Comparison of crater geometry obtained from different modeling approaches.

Simulation approach		Current model	Constant properties without melting heat	Disk heat source^9,11	Point heat source⁶
Crater dimensions (µm)	Radius	39	36.5	157.5	61
	Depth	31.5	33	52	61
Relative prediction error (%)	Error in radius	18.1	23.3	230.9	28.2
	Error in depth	14.1	19.6	88.4	121

Discharge conditions have been set as per Experiment 1, Table 5.

Figure 12.

Crater cavity formation: (a) circular disk heat source model and (b) point heat source model (I = 40 A, T_dis = 20 µs and V_dis = 20 V).

It should be noted that these discrepancies in verification errors arise due to the simplistic forms of heat flux functions proposed for the disk and point heat source models which are far away from its more realistically complicated form derived and employed in our simulations, that is, the Gaussian distributed (equation (8)). For the disk heat source model^9,11

q = \frac{F_{c} VI}{π r_{c}^{2}}

(14)

and for the point heat source model⁶

q = \frac{F_{c} VI}{2 π r^{2}}

(15)

where q is the heat flux, V is the discharge voltage, I is the discharge current, F_c is the portion of discharge power entered to the cathode (workpiece), r_c is the radius of plasma channel at the cathode surface and r is the radial distance from the origin. Furthermore, the relationship of plasma channel growth specifying the size and shape of discharge channel takes the form as r_c = 0.788t ^0.75, depending only on the pulse duration, while in this study a more comprehensive formula has been used (equation (2)) which takes into account the effects of both pulse time and discharge current, hence providing a more efficient platform leading to more accurate results.

The last issue that needs farther investigation is the fraction of discharge energy entering the workpiece. A comparative analysis is needed to unveil the effect of different values of F_c on the predictive accuracy of our FE model. Table 7 lists the different values of crater dimensions obtained using various F_c s ranging from 0.05 to 0.5. As a general rule, the higher the F_c , the bigger would be the crater dimensions due to the higher amount of discharge power transferred to the cathode workpiece. Moreover, the values of radius errors present a decreasing trend by increasing the F_c from 0.05 to 0.4, whereas the depth errors show an ever-increasing style after the second simulation. In order to get a fair conception about the suitable numerical value assigned for this fraction, the arithmetic mean of the two verification errors can be considered as the decisive measure. Here, as can be seen, if the F_c is chosen a value between 0.1 and 0.2, the average error will get its lowest amount (almost 16%) compared to other cases. This is also in agreement with DiBitonto et al.⁶ who recommended that the energy distribution factor should be about 18% for cathode to provide good correlation between theoretical and experimental results. As mentioned before, the same amount has been used in our model building and its validity confirmed here.

Table 7.

Effect of different values of cathode fraction of discharge power on the accuracy of FE predictions (Experiment 1, Table 5).

F_c	Simulated crater dimensions		Relative verification error (%)
	Radius (µm)	Depth (µm)	Error in radius	Error in depth	Average error
0.05	29	19.5	39.1	29.4	34.3
0.1	35	26	26.5	5.8	16.2
0.2	41	33	13.9	19.6	16.8
0.3	44.5	37	6.5	34.1	20.3
0.4	47	40	1.3	44.9	23.1
0.5	49	42.5	2.9	54	28.5

Parametric analysis

Having verified the FEM simulations, it is now imperative to perform a parametric analysis to investigate the effects of most dominant parameters modulating the discharge energy, that is, discharge current and discharge duration (spark on-time), on crater size (radius and depth). The results would be of paramount importance predicting the productivity (MRR) of EDM process as this characteristic is directly proportional to the size of crater volume. It also helps understand the quality of how a parameter affects the MRR while working under different machining conditions.

Effect of discharge current on crater dimensions

Figure 13(a) illustrates the effect of discharge current on crater radius and depth while keeping the other parameters constant (T_dis = 100 µs and V_dis = 40 V). As is clear, both the radius and depth get larger values varying from 51.5 to 89 µm and 42.5 to 79 µm, respectively. This is due to the higher levels of discharge energies released during sparking when discharge current increases from 10 to 50 A. However, it is remarkable to note that the ratio of the developed crater depth to radius (d/r) increases from 0.83 to 0.89. It implies that the crater depth grows at a higher rate than that of crater radius indicating that the melting border advances into the workpiece as a result of high discharge current to avoid re-solidification of molten material. In other words, the material ejection efficiency increases progressively by increasing discharge current with more material removed instead of deposited over the work surface. Hence, thinner recast layer can be expected while increasing discharge current.⁴⁴

Figure 13.

Effects of discharge parameters on crater dimensions: (a) discharge current and (b) discharge duration.

Effect of discharge duration on crater dimensions

The effect of discharge duration on crater dimensions (radius and depth) at a steady level of discharge current (30 A) and voltage (40 V) has been shown in Figure 13(b). Similar to the effect of discharge current, the crater size increases by providing larger discharge energies resulted from longer discharge times. Nonetheless, unlike the discharge current, the quotient between crater depth and radius (d/r) decreases from 0.88 at T_dis = 10 µs to 0.84 at T_dis = 200 µs. It conveys that the crater depth increases at a lower rate than the crater radius does, that is, the melting front recedes while discharge duration increases which means more melted material gets re-solidified and spread over the work surface instead of removed effectively. This phenomenon can best be ascribed to the drop of density of discharge energy at the spot of sparking due to the expansion of plasma canal radius at prolonged pulse durations, thereby lowering the plasma flushing efficiency and hence MRR.⁴³

Based on the outcomes of the above analyses, it can generally be concluded that the sparks having higher discharge currents along with shorter discharge durations are more favorable in terms of the amount of the material removed from the workpiece. Higher discharge currents accompanied with shorter pulse durations provide more effective sparks possessing greater electrical discharge density, and hence impulsive force leading to more material ejected instead of re-solidified over the crater surface. This result, solely obtained by the simulation and parametric analysis of the size of single-spark crater cavity, is indeed in agreement with the results obtained by other researchers who conducted experimental works on the real state multi-spark EDM process.⁴⁵

Model extension to real state multi-discharge EDM

The ideally theoretical material removal or erosion rate due to the cumulative effects of single-spark discharges, MRR_Sim (mm³/min), considered in a series of consecutive pulses during the continuous operation of EDM process, can be computed as

MR R_{Sim} (m m^{3} / \min) = \frac{60 \times 10^{- 3} \times C_{v}}{T_{on} + T_{off}}

(16)

where $C_{v}$ is the crater volume (µm³); T_on and T_off are the pulse on-time and the pulse off-time, both in microseconds, respectively.

To examine the generalization capabilities of our single-spark model to the real state EDM process in which series of successive discharges take place during machining, a comparative study has been made here. The experimental and theoretical values of MRRs have been adopted from DiBitonto et al.,⁶ which is regarded as one of the pioneering research works on EDM containing benchmark data quite often referred in EDM literature. In addition, the numerically simulated MRRs (MRR_Sim ) have been calculated as per equation (16) through the development of a program written in FORTRAN to carry out the required calculations estimating the crater cavity based on the outputs of FEM software under the same discharge conditions.

Table 8 lists the amounts of MRRs obtained using the three different aforementioned approaches in ascending order of discharge energies, whereas Figure 14 shows these MRR values against discharge energies.

Table 8.

Comparison of MRR between experimental, theoretical and the present FEM model.

Experiment number	Machining conditions			Discharge energy W (mJ)	MRR (mm³/min)
	Current I (A)	Pulse on-time T_on (µs)	Pulse off-time T_off (µs)		Experimental⁶	Analytical solution⁶	Current research
1	2.34	5.6	1	0.33	0.3	13.82	8
2	2.83	7.5	1.3	0.53	1.6	17.26	9.08
3	3.67	13	2.4	1.19	3.1	21.78	11.99
4	5.3	18	2.4	2.39	8.4	35.58	23.02
5	8.5	24	2.4	5.1	23.2	63.79	41.22
6	10	32	2.4	8	32	77.18	53.84
7	12.8	42	3.2	13.44	50.5	100.33	65.05
8	20	56	3.2	28	89.7	164.65	110.26
9	25	100	4.2	62.5	125	207.2	139.82
10	36	180	4.2	162	226	304.56	237.24
11	44	240	5.6	264	246	373.09	253.01
12	58	420	7.5	609	346	494.03	355.11
13	68	560	10	952	559	579.7	400.83

MRR: material removal rate.

V_dis = 25 V. Workpiece: steel; tool: copper; dielectric: kerosene; discharge energy (W) = V_dis × I × T_on .

Figure 14.

Comparison of MRR for cathode erosion among experimental, analytical and FEM models’ results.

As can be observed, the current numerical model exhibits the same behavior as that offered by either experimental or theoretical one following the same trend of variations. It is interesting to note that the current numerical model presents a much closer margin than the analytical approach with respect to the experimental data except for the last point where a large deviation can be seen. This is because of the fact that in this position, the discharge energy has exceeded 670 mJ, and hence, the relation introduced for plasma channel expansion (equation (2)) is no longer valid. Therefore, it can generally be confirmed that the current FE model is more suitable and accurate than the theoretical model in predicting the MRR under different levels of discharge energies. It should be noted that the theoretical results⁶ have been obtained based on the analytical solution of equation (1) assuming a point heat source model, which results in hemispherical crater cavity. Besides, all the material properties (thermal conductivity, specific heat and density) have been supposed to be constant over the whole temperature range involved, and that the inclusion of latent heat of fusion has also been neglected, while in this research their roles have all been considered to make the numerical model as comprehensive and reliable as possible.

On the other hand, it has been shown that the volume of the material removed as a result of single discharge experiment is less than that simulated by FEM due to the formation of recast layer but still higher than the one attributed to its equivalent counterpart as an isolated single pulse within a series of on-going discharge pulses occurring under the real conditions of EDM process.²⁸ The discrepancy between gap environment in single-spark and multiple discharge cases, and boundary conditions which give rise to different temperature distributions and hence different MRRs can be mentioned as the key reasons for this disparity in removal amounts which needs further research. Moreover, various kinds of pulses such as normal, short, open and arc can take place during EDM leading to no or different amount of material removal, whereas their effects cannot be incorporated into the numerical simulation. These are the main causes responsible for the errors between the MRRs obtained by our simulations and the real experimental values.

Conclusion

In this research, attempts have been made to develop a more comprehensive and detailed physics-based electro-thermal simulation, experimentally validate and parametrically analyze the effects of discharge parameters on the characteristics of single-spark erosion which could, in turn, manifest the true process mechanism of EDM. To those ends, several more realistic assumptions such as the Gaussian distribution of heat flux within plasma channel, spark radius equation based on discharge current and duration, the latent heat of fusion, temperature-dependent material properties, convection heat transfer between workpiece and dielectric and the percentage of discharge energy transferred to the workpiece have all been applied to make the simulation strategy as complete and accurate as possible.

In brief, the following can be mentioned as the main findings of this research:

Due to the temporal growth of plasma channel radius and the subsequent decrease in electrical discharge density, the peak temperature of workpiece, at the center of plasma channel, decreases with the passage of discharge time.

The volume of melted and evaporated material increases by going more into the discharge time as more heat penetrates into the bulk material.

The developed FEM model presents much better predictive accuracy compared with previous models based on disk or point heat sources. It has a maximum verification error of 18.1% and 14.1%, whereas the disk heat source possesses 230.9% and 88.4%, and the point heat source 28.2% and 121% in estimating crater radius and depth, respectively.

The inclusion of variable materials properties and latent heat of melting has a significant effect on the accuracy of simulation results. It reduces the verification errors from 23.3% to 18.1% and 19.6% to 14.1% in predicting crater radius and depth, respectively.

Parametric analyses indicate that the molten depth advances as the discharge current increases while it recedes as the discharge duration increases. The first implies that steadily deeper crater cavities are created by progressive increase in discharge current while the latter implies the formation of more re-solidified materials on the crater surface.

When extending the results of single-spark model to predict the MRR of real state multi-spark EDM, a closer proximity with experimental results than those of analytical solution can be observed. This is undoubtedly due to the inclusion of more realistic assumptions being in line with the real-life EDM considered in this study.

Footnotes

Declaration of conflicting interests

The authors declare that there is no conflict of interest.

Funding

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

References

Jain

. Advanced machining processes (12th reprint). New Delhi, India: Allied Publishers, 2010.

Shrivastava

Dubey

. Electrical discharge machining–based hybrid machining processes: a review. Proc IMechE, Part B: J Engineering Manufacture 2014; 228(6): 799–825.

Rajurkar

Sundaram

Malshe

. Review of electrochemical and electrodischarge machining. Procedia CIRP 2013; 6: 13–26.

Schumacher

. After 60 years of EDM the discharge process remains still disputed. J Mater Process Tech 2004; 149: 376–381.

Singh

Ghosh

. A thermo-electric model of material removal during electric discharge machining. Int J Mach Tool Manu 1999; 39: 669–682.

DiBitonto

Eubank

Patel

. Theoretical models of the electrical discharge machining process. I. A simple cathode erosion model. J Appl Phys 1989; 66(9): 4095–4103.

Abbas

Solomon

Bahari

MdF

. A review on current research trends in electrical discharge machining (EDM). Int J Mach Tool Manu 2007; 47: 1214–1228.

Eubank

Patel

Barrufet

. Theoretical models of the electrical discharge machining process. III. The variable mass, cylindrical plasma model. J Appl Phys 1993; 73(11): 7900–7909.

Yeo

Kurnia

Tan

. Electro-thermal modeling of anode and cathode in micro-EDM. J Phys D: Appl Phys 2007; 40: 2513–2521.

10.

Panda

Bhoi

. Electro-discharge machining—a qualitative approach. Mater Manuf Process 2006; 21: 853–862.

11.

Patel

Barrufet

Eubank

. Theoretical models of the electrical discharge machining process. II. The anode erosion model. J Appl Phys 1989; 66(9): 4104–4111.

12.

Kojima

Natsu

Kunieda

. Spectroscopic measurement of arc plasma diameter in EDM. CIRP Ann: Manuf Techn 2008; 57: 203–207.

13.

Descoeudres

Hollenstein

Demellayer

. Optical emission spectroscopy of electrical discharge machining plasma. J Phys D: Appl Phys 2004; 37: 875–882.

14.

Salah

Ghanem

Atig

. Numerical study of some aspects of electric discharge machining process. Int J Mach Tool Manu 2006; 46(7–8): 908–911.

15.

Shankar

Jain

Sundararajan

. Analysis of spark profiles during EDM process. Mach Sci Technol 1997; 1(2): 195–217.

16.

Ahn

Chung

. Numerical analysis of the electro-discharge machining process for alumina-titanium carbide composite-II. Unsteady state approach. Korean J Chem Eng 2002; 19(4): 694–702.

17.

Marafona

Chousal

JAG

. A finite element model of EDM based on the Joule effect. Int J Mach Tool Manu 2006; 46(6): 595–602.

18.

Somashekhar

Mathew

Ramachandran

. Electrothermal theory approach for numerical approximation of the µ-EDM process. Int J Adv Manuf Tech 2012; 61(9–12): 1241–1246.

19.

Nowicki

Dmowska

Podolak-Lejtas

. A new method of investigating crater and flash made by individual discharge using scanning profilometers. Wear 2011; 270(3–4): 121–126.

20.

DiBitonto

. Measurement of EDM light emission: preliminary test results. AGIE corporation internal report, Locarno, 27 March 1986.

21.

Erden

. Effect of materials on the mechanism of electric discharge machining (EDM). J Eng Mater: T ASME 1983; 105: 132–138.

22.

Natsu

Shimoyamada

Kunieda

. Study on expansion process of EDM arc plasma. JSME Int J C 2006; 49(2): 600–605.

23.

Ikai

Hashiguchi

. Heat input for crater formation in EDM. In: Proceedings of the 11th international symposium for electro machining—ISEM XI, EPFL, Lausanne, 17–21 April 1995, pp.163–170. Lausanne, Switzerland: Presses polytechniques et universitaires romandes.

24.

Papoulis

Pillai

SU.

Probability, random variables and stochastic processes. 4th ed. New York: McGraw Hill, 2002.

25.

Perez

Rojas

Walder

. Theoretical modeling of energy balance in electroerosion. J Mater Process Tech 2004; 149: 198–203.

26.

Motoki

Hashiguchi

. Energy distribution at the gap in electric discharge machining. CIRP Ann: Manuf Techn 1967; XIV: 485–489.

27.

Konig

Wertheim

Zvirin

. Material removal and energy distribution in electrical discharge machining. CIRP Ann: Manuf Techn 1975; 24(1): 95–100.

28.

Hinduja

Kunieda

. Modeling of ECM and EDM processes. CIRP Ann: Manuf Techn 2013; 62(2): 775–797.

29.

Yeo

Kurnia

Tan

. Critical assessment and numerical comparison of electro-thermal models in EDM. J Mater Process Tech 2008; 203: 241–251.

30.

Joshi

Pande

. Thermo-physical modeling of die-sinking EDM process. J Manuf Process 2010; 12: 45–56.

31.

Benzerga

Lhiaubet

Meyer

. Improvement of a proposed model for heat conduction in an electrode submitted to an electric discharge: application to the indirect determination of the anodic and cathodic voltage drops. J Appl Phys 1985; 58: 604–606.

32.

Lhiaubet

Meyer

. Method of indirect determination of the anodic and cathodic voltage drops in short high-current discharges in a dielectric liquid. J Appl Phys 1981; 52: 3929–3934.

33.

Dixit

Joshi

Davim

. Incorporation of material behavior in modeling of metal forming and machining processes: a review. Mater Design 2011; 32: 3655–3670.

34.

Sabbaghzadeh

Azizi

Torkamany

. Numerical and experimental investigation of seam welding with a pulsed laser. Optics Laser Tech 2008; 40(2): 289–296.

35.

Fish

Belytschko

A first course in finite elements. New York: John Wiley & Sons, 2007.

36.

Rebelo

Kornmeier

Batista

. Residual stresses after EDM-FEM study and measurement results. Mater Sci Forum 2002; 404–407: 159–164.

37.

Byron Bird

Stewart

Lightfoot

. Transport phenomena. New York: John Wiley & Sons, 2007.

38.

Panda

Bhoi

. Analysis of spark eroded crater formed under growing plasma channel in electro-discharge machining. Mach Sci Technol 2005; 9: 239–261.

39.

Van Dijck

Dutre

. Heat conduction model for the calculation of the volume of molten metal in electric discharges. J Phys D: Appl Phys 1974; 7(6): 899–910.

40.

Beck

. Transient temperatures in a semi-infinite cylinder heated by a disk heat source. Int J Heat Mass Tran 1981; 24(10): 1631–1640.

41.

Jilani

Pandey

. An analysis of surface erosion in electrical discharge machining. Wear 1983; 84: 275–284.

42.

Reddy

Gartling

DK.

The finite element method in heat transfer and fluid dynamics, 2nd ed. London: CRC Press, 2010.

43.

Assarzadeh

Ghoreishi

. Statistical modeling and optimization of the EDM parameters on WC-6%Co composite through a hybrid response surface methodology-desirability function approach. Int J Eng Sci Tech 2013; 5(6): 1279–1302.

44.

Shabgard

Ahmadi

Seyedzavvar

. Mathematical and numerical modeling of the effect of input-parameters on the flushing efficiency of plasma channel in EDM process. Int J Mach Tool Manu 2013; 65: 79–87.

45.

Assarzadeh

Ghoreishi

. Neural-network-based modeling and optimization of the electro-discharge machining process. Int J Adv Manuf Tech 2008; 39: 488–500.