Theoretical and experimental analysis of plasma radius expansion model in EDM: a comprehensive study

Electrical Discharge Machining (EDM) is an established non-conventional process, which is particularly efficient for the processing of hard-to-cut materials, in order to obtain high dimensional accuracy and surface integrity. However, in order to determine the appropriate parameters for machining novel materials, it is necessary to investigate the EDM process in depth, both by experiments and numerical models, taking into consideration the fundamental physical phenomena occurring during this process and be able to predict the surface morphology and microstructural alterations under various conditions. One of the challenging issues of EDM simulation models that still remain open is the representation of the evolution of plasma channel radius, for which various approaches have been proposed such as a linear, power law, or a more complex piecewise relation, in respect to time. Thus, in this work, the effect of different relations for the plasma channel radius evolution on energy absorption coefficient, plasma flushing efficiency (PFE), and crater morphology is compared under various conditions with a numerical model, which is also compared to experimental data. The results indicate that the energy absorption coefficient is dependent on the plasma column radius function, as slower growth of plasma channel leads to lower absorption coefficient and more efficient material removal, whereas a lower variation and different trends under different conditions were observed regarding PFE values, in respect to the power law exponent. Finally, the crater dimensions were shown to be consistently more narrow and deeper with higher exponents; thus, based on actual observations of indicative craters, it was revealed that the appropriate values for the exponent of the power law plasma radius function are below 0.25.


h melt
Artificial heat transfer coefficient that is utilized to calculate the amount of energy that results in the melting of the material on the eroding front (W/(m 2 K)) I P Pulse-on current (A) k Material's thermal conductivity (W/mK) LH Material's latent heat of melting (J/kg) n Exponent coefficient that is correlated with the plasma channel growth rate (-) Q Volumetric heat source (W/m 3 ) q conv Convection boundary heat flux (W/m 2 ) q melt Artificial, specifically defined heat flux that controls the temperature on the molten front (W/m 2 ) q pl (r,t) Plasma channel boundary heat flux (W/m 2 ) q rad Radiation boundary heat flux (W/m 2 ) r The distance from the center of the plasma channel (m) R 0 pl Maximum plasma radius (m) R pl (t) The plasma radius over time (m) Material's melting point (K) υ melt The velocity of the crater's wall due to the material erosion (m/s) V Close circuit voltage (V) ε Emissivity coefficient of steel (-) ρ Material's density (kg/m 3 ) σ(t) The standard deviation of the Gaussian distribution (m)

Introduction
EDM is one of the most popular non-conventional machining processes, which involves the use of a large number of spark discharges in order to remove material from the workpiece surface. The discharges occur as a result of a sufficient voltage difference between an electrode acting as cathode and the workpiece acting as anode, which gives rise to a plasma channel of high energy density [1]. Subsequently, high temperatures occur in the gap between anode and cathode, melting or ablating the workpiece material [2,3]. The EDM process is schematically presented. Thus, this process is particularly effective during processing of hard-to-cut materials, such as steels [4,5], titanium [6][7][8], or nickel alloys [9], as it does not involve mechanical contact between a tool and the workpiece, and it can process every electrically conductive material, regardless of its mechanical strength [10]. Various parameters can affect the efficiency of EDM, such as the pulse-on current, pulse-on time, machining voltage, duty factor, and thermophysical properties of the workpiece material. Due to the large number of parameters, researchers often conduct a considerable amount of experiments for each material and perform optimization studies [11]. EDM process has been extensively studied during the last decades, as it has been proven important for the construction of various parts for the automotive, aerospace, and biomedical industries [12]. Although this process can be successfully carried out in industrial practice, its application for advanced or novel materials is not straightforward and requires optimization based on both the efficiency of the process, indicated by material removal rate (MRR) and tool wear ratio (TWR), and on the surface integrity, indicated by surface roughness or modification of the microstructure [11]. However, due to obvious difficulties, both in directly monitoring the process and determining the procedure of the creation of heat affected zone (HAZ) and white layer (WL), reliable simulation models are required in order to provide a deeper understanding of the complex phenomena occurring in this process. In fact, modeling of EDM process is not a trivial task, as the occurrence of complex thermophysical phenomena should be taken into consideration. Although simulation models for EDM have been developed for several decades, there are issues which remain unresolved still nowadays, and thus, EDM remains not entirely understood or explained by a universal model. Furthermore, it is inevitable to adopt reasonable simplifications for the models in order to perform simulations in a feasible timeframe, given the limitations of the available computational power, but the accuracy of the simulations should be ensured by optimizing the computational model parameters based on experimental observations (Fig. 1).
As the thermal effects are dominant in EDM, most simulation models focus on modeling the effect of heat input into the workpiece by using a heat source with a Gaussian distribution [13,14]. This heat source represents a single spark, and heat transfer is mainly considered to take place by conduction [15,16]. Then, the calculated MRR resulting from the action of the spark is compared to the experimental [13,17]. Other important issues with the models are the determination of energy distribution between electrode and workpiece, which cannot be considered constant, but dependent on process parameters [18][19][20][21], and the determination of the PFE, which is relevant to the percentage of molten material removed from the workpiece by every Fig. 1 Schematic presentation of the EDM process spark. For PFE, the simplification that it is equal to 100% is definitely unacceptable [13,21], as it contradicts the formation of WL, which is evident by many experimental works.
Especially for the plasma column, although some authors neglected the variability of its spatial dimension, most authors take into account the plasma channel radius by semi-empirical formulas or even adopt time-dependent relations, as it was shown that the variation of the plasma channel radius can directly affect the simulation results, especially for shorter discharge times [22]. Semi-empirical formulas usually correlate the radius with discharge power and pulse-on time [2,23,24] or pulse-on time and discharge current, but formulas using gap voltage exist as well [25], although other parameters, such as properties of the electrode and dielectric, polarity, and workpiece material, can also affect it [20]. These formulas are developed based on the experimental results after fitting procedure, while direct measurements of the plasma column dimension, e.g., using high-speed cameras and spectrometers, are rare [26,27]. These formulas are considerably popular, can be used for a specific range of conditions, and can provide fairly good approximation of the dimensions of the plasma column, something that is fundamental for the application of heat source in the model. However, phenomena such as the rapid expansion of the plasma column during the first microseconds after its creation can only be taken into consideration by time-dependent formulas.
Regarding the modeling of plasma column expansion, apart from the obvious difficulties in observing this phenomenon directly, there is also a lack of established models for plasma dynamics describing the different states of plasma column during EDM, i.e., creation, expansion, and collapse [28]. Only a few authors have adopted comprehensive approaches relevant to the time evolution of plasma column. For example, Mujumdar et al. [28] attempted to create a model for plasma column taking into account chemical phenomena, apart from electrical, and thermal ones. Using this model, they directly predicted heat flux to the workpiece and electrode under various conditions. Eubank et al. [29] developed a model for the calculation of plasma radius, temperature, and pressure evolution using a fluid mechanics equation, an energy balance equation, a radiation equation, and an equation of state, whereas thermophysical properties were calculated by taking into consideration some fundamental reactions occurring during the process. Pandey and Jilani [30] adopted an iterative procedure for the determination of plasma channel radius based on the assumption that cathode spot temperature is constant during the pulse-on time and equal to the boiling point temperature of the electrode material. Results were obtained for various electrode materials, and this approach was considered sufficiently accurate compared to experimental data. Zhang et al. [31] determined a formula for the spark column expansion, based on a differential equation describing spark volume increase and also took into account variable heat transfer time, depending on the distance of points from the center of the spark. Chu et al. [32] presented a comprehensive methodology for the calculation of plasma column radius, temperature, and pressure including various components, such as a model for breakdown in the dielectric medium, which takes into consideration the nucleation of bubbles and production of discharge and models for the initial stages of plasma column formation and expansion stage until the end of each discharge. Dhanik and Joshi [33] also developed a model for plasma column, including the nucleation and growth of bubbles at the first stage and appropriate fluid dynamics and heat transfer models for the heating stage.
Shabgard et al. [34] used a differential equation for the calculation of plasma column radius, which took into consideration the contribution of electric field created by charged particles to the acceleration of plasma channel, the contribution of magnetic field created by the movement of the charged particles, the contributions of internal and external pressure of plasma channel, and the surface tension in the interface between plasma channel and dielectric. Especially for magnetic field-assisted EDM process, Shabgard et al. [35] also used a slightly modified model of their original one, including the effect of the additional magnetic field in the differential equation of plasma expansion. Gholipoor et al. [36] compared the results of the model presented by Shabgard et al. [34] and showed that the calculated plasma radius values are slightly higher to those of semi-empirical models.
However, in the majority of relevant works, due to computational power limits, a simpler approach was followed using a single or a piecewise function to model the time dependence of plasma column radius. In specific, the following types of functions have been already proposed in the relevant literature: linear, power law, combined power law with semi-empirical and piecewise. For the linear function, it can be noted that only a few works have used it, such as the work of Singh and Ghosh [37], who analyzed plasma formation by a fluid mechanics equation, in order to establish a thermoelectric model for the calculation of the electrostatic force on the electrode and the stress distribution on the workpiece during a discharge. In their work, plasma radius was assumed to vary linearly with time for small discharge durations.
On the other hand, the power law function is by far the most popular in the relevant literature. In these works, the assumption made is relevant to the continuous expansion of plasma column radius during the discharge time [30]. Izquierdo et al. [38] argued that modeling of plasma radius with a power law function including a constant term is in line with the common assumption that the expansion of plasma channel is abrupt for a few microseconds and then stabilizes. Similar approaches were adopted by Shao and Rajurkar [16] and Guo et al. [17]. Kliuev et al. [39] used a similar function with the addition of a second constant term which could be determined by an optimization procedure. Schneider et al. [40] also employed a function with a second constant term, which was dependent on maximum spark radius, initial radius, and discharge duration. In every case where power law was adopted, the exponent has values below 1.0. In many of these works [16,17], the exponent of the power law function is 0.75; however, it is worth noting that an exponent with a value around 0.2 has also been suggested in the works of Revaz [41] and Perez [42] based on various experimental observations, in contrast to the works which used larger exponents. More specifically, Revaz et al. [41] commented that the exponent value of 0.75 corresponds only at the very beginning of the discharge and determined that the most suitable value for the exponent was close to 0.2.
Assarzadeh and Ghoreishi [43,44] noted that the discharge channel radius is dependent on various factors and utilized a combined semi-empirical/power law relation for the growth of discharge channel radius, including an additional term dependent on discharge current. In this relation, time was varied between 0 and T on value, after which the plasma channel is assumed to collapse, and this relation was considered accurate up to a specific spark energy level. Similar approaches were used by Guo et al. [45] and Xie et al. [46]. Vishwakarma et al. [47] also introduced a relation which included other terms rather than time, such as the discharge length and discharge power, as well as empirical constants. The latter approach was also mentioned by Kumar et al. [20] and Kansal et al. [48].
Rajhi et al. [49] conducted a comprehensive comparison regarding the use of constant equivalent plasma column radius functions, time-dependent radius function, and time-dependent radius function with a semi-empirical term. The models which employed a constant radius, even if different semi-empirical functions were used, showed similar results both regarding the highest temperature and the time when it was reached, as well as similar temperature variation in respect to time. For the time-dependent models, although the time to reach the maximum temperature was similar, significant differences were obtained for the temperature profiles. Regarding crater morphology, it was shown that both crater diameter and depth varied for different radius models, and especially for crater diameter, considerably different results could be obtained. One model showed shallow crater morphology, whereas the others predicted a more hemispherically shaped one. Finally, MRR prediction was also found to be dependent on the radius model, with some models overestimating MRR and others underestimating its values, whereas the accuracy of some of the models was also limited only to cases with low or high values of discharge current.
Liu and Guo [50,51] developed a plasma radius growth model with a piecewise function employing two empirical constants. The first expression, which was a power law function in respect to time, was considered valid up to a certain critical time, and then the plasma radius was considered constant. The estimation of empirical constants and critical time was based on a considerable amount of experimental data [51]. This model reflects some experimental observations, which indicated that plasma is not continuously expanding during the entire pulse duration [51]. Li et al. [52] used also a piecewise function, but the function employed for the expansion stage was linear in respect to time.
Natsu et al. [27] performed observations with highspeed camera and confirmed that the expansion of plasma is completed after a few microseconds after the breakdown of the dielectric. Thus, they proposed a piecewise model, termed as first stage expansion model (FSEM), and compared it to a power law model which was inferior to the FSEM regarding the crater dimensions. Izquierdo et al. [53] also compared power law functions to piecewise ones (FSEM), including a linear part for an extremely brief time period and a constant part. The simulation results were compared to experimental ones regarding MRR and surface roughness indicators. At first, it was shown that a FSEM model with an optimal value for the final plasma radius can provide an excellent estimation for MRR and acceptable error for surface roughness. Thus, they considered it necessary to include the first stage of rapid expansion in their model. However, the use of an improper final value for plasma radius leads to a much higher error level for FSAE models, compared to power law ones.
From the analysis of works relevant to the modeling of plasma column radius in the relevant literature, it can be observed that there is a lack of a comprehensive comparison of the results produced between different plasma radius expansion functions for various different conditions based on experimental validation of both crater dimensions and microstructure observations. Moreover, the variation of PFE values in respect to different plasma radius expansion functions has also not been considered yet in the literature. Thus, in the present work, a numerical model for the prediction of crater morphology, HAZ, and WL formation during EDM is employed in order to comprehensively investigate the effect of the plasma column radius expansion on the results of the EDM simulations. Three different types of functions are compared, namely, constant, linear, and power law functions, and their effect on the geometry of produced craters, energy absorption coefficient, and PFE is analyzed, based also on actual microscope observations of craters produced during EDM of 60CrMoV18-5 steel under various conditions.

Materials and methods
For the current study, a FEM model was employed to simulate the material ablation during a single spark. Taking into the account the fuzzy and inherent chaotic nature of the process in microscale, the adoption of a "typicalaverage spark" approach is reasonable, justified, and up to a point necessary. This choice is in line with the literature where Klocke et al. [6] deduced that the material removal from a continuous discharge process cannot be studied in a representative way by single-discharge experiments. The current modeling methodology has already been validated regarding its adequacy to simulate the material removal mechanism during EDM [54]; thus, in the current study, and keeping in mind that the main aim is to define how the boundary conditions regarding the temporal evolution of the plasma radius affect the simulation results, the FEM model will be discussed only in brief, emphasizing on the main topic of the study.
The general model assumptions and necessary simplifications are: • All the sparks are assumed to be identical • Conduction is considered as the only mechanism of heat transfer into the material volume [6] • A Gaussian boundary heat source is utilized to model the heat flux from the plasma channel to the workpiece [21,55] • The percentage of the spark energy that it is absorbed by the workpiece depends on the machining parameters • The PFE also depends on the machining parameters • Deformed geometry feature is employed to model and simulate the material erosion mechanism • A homogenous and isotropic material with temperature-dependent thermophysical properties is considered

Governing equations and boundary conditions of the FEM model
The Fourier's law that describes the heat transfer due to conduction is expressed as with T the temperature in K, ρ the material density in [kg/m 3 ], Cp the material specific capacity in J/kg K, k the material thermal conductivity in W/mK, and Q a volumetric heat source (or a heat sink) in W/m 3 . For the current model, four boundary heat sources are considered, namely, the heat flux from the plasma channel, the heat losses due to convection and radiation, and finally an artificial, specifically defined heat flux that controls the temperature on the molten front. Thus, the boundary heat flux term of Eq. 1 can be expressed as The heat flux due to the plasma channel is modeled by a Gaussian boundary heat source that it is mathematically described as with V the close circuit voltage in V, I P the pulse-on current in A, f w the percentage of the spark energy that it is absorbed by the workpiece, r the distance from the center of the plasma channel in m, and σ(t) the standard deviation of the Gaussian distribution that can be expressed as with R pl (t) the plasma radius over time in m. As it is already mentioned in introduction and based on the relevant literature review, the estimation of the plasma radius consists of an area of active research. It has been deduced that among the different semi-empirical equations for the definition of the plasma radius depending on the machining parameters, a suitable relationship that can also be combined with deformed geometry feature is [2,54] with R 0 pl the maximum radius that the plasma channel obtains during the discharge in m and T on the pulse-on time in s. Thus, although through Eq. 5 the maximum plasma radius can be estimated, the temporal evolution of the plasma channel radius is not clear, and hence, it is also not clear how the temporal growth of the plasma channel affects the obtained results from the simulation. Hence, in order for this topic to be fully and in-depth investigated, in the current study, four different approaches for the plasma radius growth are tested. More specifically, the first and most common approach is that the plasma channel has a constant radius for the whole discharge duration, and thus The three other approaches assume an exponential growth of the plasma radius over time but with different growth rates, and hence, three different exponents are adopted and tested. The time-dependent plasma radius can be expressed as with n the different exponent that is correlated with the plasma channel growth rate which has the following values: 0.05, 0.15, and 0.25. These values were chosen in accordance with the relevant literature which indicates that an appropriate value for the exponent of the power law is around 0.2 [35,36]. However, given that this subject was not thoroughly investigated, a wider range of exponents was adopted.
The plasma radius as function of time for the four different plasma radius models is depicted in Fig. 2. It is deemed that these four models can adequately describe the different approaches that have been suggested in the relevant literature regarding the plasma channel growth during the discharge time.
The heat losses due to convection and radiation are respectively calculated as with q conv and q rad the heat fluxes due to convection and radiation respectively in W/m 2 , h diel the heat transfer coefficient between the workpiece and the dielectric fluid with value 10 5 (W/(m 2 K)) [56], and ε = 0.75 the typical emissivity coefficient of steel. Finally, T amb and the T diel are the ambient temperature and the dielectric temperature respectively that are both considered equal to 293.15 K [K].
The material melting rate and the material removal rate (i.e., the erosion rate) are directly correlated through the energy that is required in order for the eroding front to switch over from solid into liquid. In the current modeling amb method, two distinct process phases are regarded: the 1st of the material removal and the 2nd of the WL formation. During the 1st phase, it is considered that the material on the eroding front remains just over the melting temperature, whilst the whole molten material is removed. In the 2nd phase, no material removal occurs, the material overheats (i.e., reaches temperatures above its melting point), and thus, the WL is gradually formed. The 1st phase of the material erosion can be modeled and simulated by employing an artificial, specifically defined, boundary heat flux: with q melt the heat flux that results in the material melting on the eroding front in W/m 2 and T ph.ch. the material's melting point in K. The h melt (W/m 2 K) represents the artificial heat transfer coefficient that is utilized to calculate the amount of energy that results in the melting of the material on the eroding front. The q melt is used in Eq. 11 to define the material erosion rate, i.e., the eroding velocity of the crater's wall: with υ melt the velocity of the crater's wall due to the material erosion in m/s and LH the material's latent heat of melting in J/kg.
During the 2nd phase, zero erosion is considered, and thus The time point for which the h melt is zeroed (and hence the υ melt is also zeroed) is when the simulation crater's volume equals to the experimentally defined mean crater's volume. The fine tuning of the energy absorption coefficient f w , which is depending on the machining parameters (see Eq. 3), is based on the experimental results about the average white layer thickness (AWLT). More specifically, the simulation results must not only be in line with the experimental regarding the crater's volume, but also predict and calculate a similar AWLT. Thus, the f w is defined in such a way that at the end of the pulse-on time, the formed crater has the foreseen volume and the overheated area has the width of the AWLT. The AWLT is experimentally measured by dividing the area of the WL with its respective length. The AWLT of the simulation results is calculated in a similar way. In Fig. 3, the aforementioned measurements for experimental and simulation results are depicted.
Finally, the PFE, and based on its definition, is calculated by assuming that during the discharge the entire molten material is removed by the workpiece. Hence, to estimate the PFE, the material erosion rate (i.e., the υ melt ) is not zeroed during the pulse duration, and thus, an optimal maximum crater volume is defined. The PFE results from the quotient of the real crater volume that has been previously calculated to the optimal one.
The simulations were carried out for four different machining parameters' combinations, while the model validation and fine tuning was done based on a previously conducted extensive study regarding the machining of 60CrMoV18-5 (CALMAX) steel with EDM [11]. More specifically, the experiments were conducted in an Agie-Charmilles Roboform 350 Sp industrial-type die-sinking CNC EDM machine, by employing a rectangle copper electrode of nominal dimensions 14 × 20 mm. The workpieces were slices of CALMAX steel, while synthetic hydrocarbon oil (kerosene) was utilized as dielectric fluid that was properly channeled into the working tank through a flushing nozzle with nominal pressure of 0.7 MPa. The close circuit voltage and the duty factor were kept constant at 30 V, and 50% respectively, while a 0.5-mm nominal depth of cut was defined. Finally, the machined surfaces were observed by using a Keyence VHX-7000 optical microscope. For the measuring of the crater dimensions, proper software was employed, i.e., Adobe Photoshop® and ImageJ. The material thermophysical properties that were used in simulations are listed in Table 1, whilst the aforementioned machining conditions are listed in brief in Table 2.
The selection of CALMAX steel for the validation of the model was done primarily based on the CALMAX properties and secondary because of the prior authors' knowhow in machining CALMAX with EDM. Namely, CALMAX steel is characterized by high toughness, good wear resistance, and polishability, while it is widely utilized in molding and die applications. Moreover, EDM is a process that is commonly employed for the machining of CALMAX steel and other tool steels. In addition, authors have carried out an extensive experimental study regarding the machining of CALMAX steel with EDM, allowing them to have a more comprehensive view of the relevant topic and the specific behavior of CALMAX steel during its machining with EDM. More specifically, by selecting these specific machining parameters, a wide range of per pulse energies is covered, i.e., from 1.92 mJ per pulse up to 51 mJ per pulse. This range of machining parameters can be reasonably considered representative for CALMAX steel since the MRR and the Ra between the aforementioned machining parameters combinations have been increased approximately by 1096% and 310%, respectively. Thus, and keeping in mind that during the planning of the EDM process, except of the MRR, a number of other parameters have also to be taken into consideration (e.g., surface quality, WL thickness, surface roughness, tool wear), this range of machining parameters can be considered realistic, representative, and with practical value as well.
Regarding the FEM model parameters, axial symmetry was used; hence, 2D axisymmetric models were solved. The control volume height and width were set equal to R 0 pl and 1.2·R 0 pl respectively, ensuring an efficient control volume where boundaries do not affect the obtained results. Triangular mesh of variable density was employed, with a finer discretization on the top boundary, where the plasma heat source acts and on the vertical boundary that handles the geometry deformation due to material erosion. Finally, a time-dependent fully coupled solver, with an absolute tolerance of 0.1 and a relative tolerance of 0.01, was utilized for the numerical solution of the model.

Results and discussion
The simulation results for the energy absorption coefficient are listed in Table 3. Concerning the energy absorption coefficient f w and as it follows by the obtained data that are plotted in Fig. 4, it has an expected behavior (see the average values in Fig. 4)  Table 2 Machining parameters

Machining parameters
No. I P (A) T on (μs) The machining was carried out with hydrocarbon oil (kerosene) as dielectric fluid and with a copper electrode. For the specific machining conditions, the close circuit voltage was kept constant at 30 V, and the duty factor was also constant at 0.5.  in respect to the different machining parameters. More specifically, for low pulse-on current and pulse-on time, low coefficient values are calculated, while for the higher pulseon current and pulse-on time, the absorption is almost doubled. This kind of dependence and behavior of the absorption coefficient from the machining parameter combination is typical in EDM, and it has been previously recorded for other alloys and materials as well [1]. In regard to the influence of the plasma model on the absorption coefficient, it is deduced that the slower the plasma channel growth is considered, the lower absorption coefficient is calculated. The higher f w is defined for the constant plasma channel radius, while the lower is defined for the exponential model with exponent value of 0.25. This can be reasonably attributed to the more efficient material removal for higher power densities. The slower plasma channel growth results in an increased average power density, and thus, the same material removal can be achieved with a lower power, i.e., lower absorption coefficient. However, it is important to be noticed that the difference between the highest and the lowest f w calculated values for the same machining parameters combination (i.e., difference between the constant plasma channel radius model and the exponential model with exponent value of 0.25) is less than 15%, with an average difference of 12.8%. This is reasonable keeping in mind that in fundamental level of physical mechanisms, the whole process is governed by a power balance; thus, the required power to melt a specific volume of material cannot significantly change, since a similar power source and the same physical mechanisms are utilized. Hence, it can be said as a rule of thumb that although the obtained results show a high level of consistency, the slower the plasma channel expansion is considered by each plasma channel radius model, the lower absorption coefficient is expected to be defined. The PFE is directly related to the absorption coefficient, but with inverse trend in respect to the machining conditions combination. An increase of the f w results in a decrease in PFE and vice versa, implying that the material removal mechanism becomes less efficient as the workpiece absorbs a higher portion of the spark power. This can also be supported by the fact that, in general, more intense machining parameters lead to an increased AWLT, a result of the partial only removal of the molten material from the crater's area. These phenomena and mechanisms have been extensively discussed in some previous studies [1,11]. Regarding the impact of the plasma channel radius model on the calculated PFE, and as can be seen from the data in Τable 3 and the respective plot of Fig. 5, the plasma channel radius model has a minor and vague impact on the calculated PFE, as the observed trends are different for different experimental cases. In general, the deviation between the lowest and the highest PFE values for the same machining parameters' combination is low, with a maximum deviation of 12.4% for the 13 A and 50 μs and an overall mean value of 8.8%. Thus, as a general conclusion, it can be deduced that the plasma channel radius model does not have a significant effect on the calculation of PFE, which is mainly affected and correlated with the absorption coefficient and the AWLT.
The plasma channel radius model has a clear impact on the formed crater geometry. Based on the obtained data that are listed in Table 4, it is easily concluded that the slower the plasma channel growth is considered, a narrower and deeper crater is formed. This is reasonable and can be attributed to the higher average power concentration during the discharge simulation. Taking into consideration that the slower plasma channel growth does not result in a sizeable decrease of the power absorption coefficient, a comparable amount of power acts more focused, resulting in the increase of the crater's depth and the simultaneous decrease of the crater's width. Additionally, the higher power concentration for the slower plasma channel radius growth and the consequent formation of a deep and narrow crater creates a "self-amplified" mechanism that further increases the depth of the crater, since the discharge energy is "trapped" and acts on the tapered crater area.
The aforementioned self-amplified mechanism is confirmed by the crater's dimension ratio (i.e., the width to depth ratio) that is presented in Fig. 6. As the plasma channel radius growth rate decreases (i.e., exponential model with higher exponent), the crater's dimension ratio has an exponential increase leading in a drilling-type of crater. This is also clearly inferred by observing the crater geometry for pulse-on current value of 17 A, pulse-on time of 100 μs, and for the different plasma channel models that are presented in Fig. 7. The craters' morphology gradually changes from a shallow and wide crater for constant plasma channel radius to a deep and narrow crater for the exponential plasma radius with an exponent of 0.25.
Hence, a reasonable question is raised, regarding the limit that the exponent may have in order for the simulation results to realistically model the EDM process and the underlying physical mechanisms. Keeping in mind the Higher values of the exponent will result in nonrealistic crater geometries. In order to prove this argument, although the predicted crater dimensions in the aforementioned cases seem realistic, an extreme case with a power law function with an exponent value of 1.0 for the plasma channel radius was also simulated. In Fig. 8, the crater geometry for exponent equal to 1 (i.e., linear increase of the plasma channel radius) is presented, for the lowest discharge energy (i.e., 5 A and 12.8 μs). Clearly, the morphology of the crater is not representative of the EDM process, since an extremely deep and narrow crater is formed. Thus, the argument that only exponents less than or close to 0.25 should be employed is justified, confirming similar results in the relevant literature [35,36]. Finally, to further support the above conclusion and to provide a robust validation of the model, experimental data regarding the craters' characteristics were utilized. In Fig. 9, the machined surfaces for the different machining parameters' combinations are presented. Although it is challenging to define a single crater since the process is chaotic and there is a superposition of the craters, in Fig. 9, a number of distinct craters were defined. The geometrical characteristics of those craters are in line with the simulation results that are Fig. 7 The different crater geometries for I P equal to 17 A and T on equal to 100 μs for a constant plasma radius, b exponential plasma radius with 0.05 coefficient, c exponential plasma radius with 0.15 coefficient, and d exponential plasma radius with 0.25 coefficient  Table 5. The width and depth of the craters gradually increase for higher per pulse energies, remaining though always in the range of the simulation values. It is important to notice that the modeling and simulation regards a typical average crater, while the fluid mechanics phenomena that take place due to the plasma channel gradients of pressure and form the craters rim are not considered. Hence, a small and acceptable deviation is expected between the experimental and the simulation results. Thus, as a general conclusion and taking into the account the experimental results as well, it can be deduced that the presented models for the plasma channel radius are in the realistic range of exponents, and  they also can be utilized for the further study of the process as well as the crater formation mechanisms, by adopting all of the presented plasma channel models under a statistically driven distribution. In Table 6, a comparison between the simulation craters' width and the experimentally measured is presented. At first, it has to be noticed that the experimental measuring of the craters' width is challenging since craters overlap is taking place, and thus, the definition of a single crater geometry is very difficult. Additionally, the stochastic and chaotic nature of the process results to a noticeable deviation in the craters' width. However, despite these inherent difficulties and based on the results of Table 6, it is deduced that the simulation results are in line with the experimental ones, having a less than 16% deviation. Thus, and as an overall conclusion, it can be said that the predicted craters are of comparable width with the experimentally measured ones, while the deviation in the values of the experimental results suggests and implies the feasibility of the appropriate utilization of different plasma models under a statistically driven distribution for more realistic and reliable results.
At this point, it has also to be mentioned and clarified that the deduced conclusions have both a quantitatively and a qualitative interpretation. More specifically, the exact values of the energy absorption coefficient, the plasma flushing efficiency, and the craters' dimensions are defined for the specific alloy (CALMAX steel) and the aforementioned machining conditions and parameters. It has been proved and established in the relevant literature that different alloys have different behavior during the EDM process; thus, these indexes have to be calculated and determined case by case. In general, alloys that belong to the same class or have similar machinability by EDM are expected to have close and comparable behavior; nevertheless, any arbitrary generalization has to be avoided, especially regarding entirely different material classes (e.g., between ferrous alloys and titanium alloys). On the other hand, there are some fundamental and basic qualitative conclusions that have a far more general applicability. Namely, and as it has already been mentioned, the exponent in plasma radius calculation formula must have a low value; otherwise, nonrealistic geometries and results will occur. Moreover, and although the energy absorption coefficient, the plasma flushing efficiency, and the craters' dimensions as absolute values are expected to differ in respect to the different alloys, they are also expected to follow common trends in respect of the machining parameters. Thus, in short, it can be said that the deduced conclusions should be considered as a robust guideline and methodology for further and in-depth study and investigation of the EDM process, and secondarily and with the required caution to be adopted as absolute values.

Conclusions
In this work, a 2D axisymmetrical FE numerical model was used to determine the suitability of various functions regarding the modeling of plasma column radius expansion during die-sinking EDM process. More specifically, for the modeling of expanding plasma column radius, power law functions, with three different exponents, were compared to a constant radius function. The model was based on the experimental MRR and WL formation and comparison with experimental results under various discharge current and pulse-on time values allowed deducing several interesting conclusions.
Regarding the values of energy absorption coefficient, it was found that when a power-law function with high exponent values is adopted for the plasma column radius, lower values are expected, indicating that when the plasma column growth rate is slower, material removal is more efficient due to higher average power density. Moreover, although the difference between the energy absorption coefficient values determined with constant function, which corresponds to a power-law function with an exponent of zero and power-law function with the largest exponent, was less than 15%, the trend of absorption coefficient decrease with exponent value increase was clear.
PFE values showed a minimal variation, namely, below 15%, with power-law exponent increase, and different trends were observed under different process conditions. Thus, it can be concluded that the choice of power-law exponent has generally negligible effect on the PFE. Crater dimensions change considerably with different power-law exponents for the plasma column radius expansion function. More specifically, the higher the exponent is, the lower the width and the higher the depth of the crater becomes, with an exponential growth of the depth/width ratio being also observed, which is unrealistic. These observations were attributed to the higher average power concentration when the growth rate of plasma column radius is lower. In the case of exponent value of 1, implying a linear increase of plasma column radius, the results are similar to those of drilling, with a high value of depth/width ratio, which seems unreasonable.
Based on observations of crater morphology from actual experiments, it was found that the average dimensions of the craters are similar to those predicted from the simulations, something that further proves that the choice of exponent values below 0.25 is reasonable. This suggestion is particularly useful for further investigations of the complex, multiphysics phenomena occurring during EDM.
Author contribution M.L.P., P.K.O., and N.E.K. carried out the experiments and formal analysis of the original draft, while M.T. and A.P.M. as supervisors discussed the methods, conclusions, and provide the final review and editing. Finally, all the authors reviewed the paper.

Declarations
Ethical approval Not applicable.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.