An approach to modelling evaporation pulsed laser drilling and its energy efficiency

Laser drilling is a thermal process with relatively low energy efficiency since the material removal mechanism is mostly based either on melting or on vaporization. Aiming at the investigation of the laser drilling efficiency, a theoretical both analytical and numerical study of evaporation pulsed laser drilling is presented. The analysis is based on a linear approximation of the temperature profile and separates the process into three phases, those of the heating, the melting and the vaporization. Based on these models, the energy efficiency and its dependence on the process parameters have been investigated and selection of relevant process variable guidelines, towards improving energy efficiency, are given. Moreover, the physical mechanisms responsible for most of the energy losses are analysed and classified according to their importance.

consumed in the industry is wasted. An action plan for the reduction of energy consumption [3] requires the study of manufacturing processes' energy efficiency, especially for high-energy-consuming processes. Under this prism, manufacturing processes have to be assessed/optimized in terms of their energy efficiency [4,5].
Laser drilling [6,7] is a rather accurate, non-conventional process with a variety of applications. The major advantage of laser drilling, in comparison with the conventional methods, is the small size and aspect ratio (up to 1:20) of the hole that can be created. Either continuous or pulsed laser sources can be used in such applications. The pulsed laser sources present some advantages, such as less plasma generation [7]. The major physical phenomenon resulting in the creation of the hole is the conversion of light power into heat. The current study focuses on the case of high intensity values, resulting in the evaporation of bulk material during percussion drilling, also referred to as pulsed laser drilling [6].
Percussion drilling delivers successive laser pulses to the same spot. A major drawback of percussion drilling is the formation of a recast layer. Additional energy is being used in order to evaporate the re-solidified material contributing to the low energy efficiency of the process.
Several models of the laser drilling process have been reported in the literature. In order to capture the physical phenomena, it is required that the temperature field [8] induced inside the workpiece by the laser source be modelled. Several attempts are reported in the literature, either analytical [9][10][11] or numerical methods [8][9][10][11][12][13][14][15][16], with the relevant advantages and disadvantages. A three-dimensional analytical model for laser grooving is presented in [17], while a one-dimensional analytical model for drilling is presented in [18][19][20]. In [18], the authors focus on the shape formation of the hole, while in [19], it is assumed that drilling occurs due to material evaporation neglecting the intermediate liquid state. In [20], they assume that the hole is formed from the ejection of melted material utilizing an assisting gas jet.
In order for the energy efficiency of a manufacturing process to be increased, a definition of the term "energy efficiency" is required [21,22]. The majority of the energy efficiency studies focus on factory level [23][24][25][26][27][28][29]. Most of these studies make general assumptions of each machines average energy consumption. The energy consumption of processes is rarely known due to the insufficient existing infrastructure and the missing measuring devices [30]. Thus, the outcome of these models is very difficult to be transformed into a strategy towards increasing energy efficiency. At process level, few studies exist for conventional manufacturing processes, such as an empirical approach for grinding process in [31] or the energy efficiency investigation of turning in [32].
The laser processes are generally characterized by energy efficiency lower than that of conventional processes; however, as indicated in [33,34], they pose a different ecological advantage due to the limited use of consumables. Other studies [35,36] conclude that the energy efficiency of such processes can be significantly improved as a result of their extremely short pulse duration; the heat diffusion is confined, and the heat-affected zone (HAZ) is rather limited. This localized heating, in each laser pulse, results in more precise machining results compared with the ones obtained from longer laser pulses. Furthermore, different studies investigate the dependence of energy efficiency on the laser beams geometrical characteristics, such as the spot geometry [37].
The current study presents an analytical and numerical approach of evaporation laser drilling, towards the theoretical specification of the process parameters that maximize energy efficiency. The phenomenon is described as a function of the processing time. The analytical approach, based on some simplifying approximations, provides simple expressions that describe the process and its energy efficiency. The numerical approach is mainly used for justifying that the approximations introduced in the analytical model do not have a significant effect on the accuracy of the results.
The analysis of the process is separated in three phases, those of the heating, the melting and the vaporization. During the heating phase, a laser pulse hits the surface of the workpiece, thus increasing the surface temperature until it reaches the melting temperature. The melting phase begins as soon as the melting temperature has been reached. During the melting phase, more and more material is melted, until the surface temperature has reached the evaporation point. Finally, during the vaporization phase some material is evaporated and thus, removed from the workpiece.

Time profile of the laser beam intensity
It is assumed that the laser turns on instantly at the beginning of every pulse and turns off instantly at the end of every pulse as expressed in Eq. (2.1).
where (t) is Heavisides step function, t p is the pulse duration and τ is the pulsing period. I 0 is the laser beam intensity when the laser is on, or where P is the mean laser beam power and r 0 is the beam radius at the surface, where laser intensity is considered. The time profile of the laser beam intensity is shown in Fig. 1.
A useful definition is also the absorbed power density, which equals where R is the workpiece material reflectivity at the laser wavelength.

Laser defocusing
When drilling starts, the surface that the laser beam hits is moving with time. Thus, defocusing of the laser beam becomes an issue. The radius of the laser beam obeys a hyperboloid hourglass distribution [7] given by where r 0 is the beam radius at the focal plane, M is the beam quality parameter, λ is the laser wavelength and δf is Fig. 1 The time profile of the laser beam intensity the z coordinate of the focal plane. Thus, the absorbed laser intensity when the erosion front is at depth z is given by (2.5) To sum up, the absorbed power density as a function of time and depth is given by where I 0a (z) is given by (2.5).

The problem
During the laser drilling process, the depth of the region of the workpiece that has been affected by the heat flow, at any given time, is very small. The reason is the laser beams high intensity that causes the rapid evaporation of the greatest part of this region. Thus, the diffusion of heat in directions perpendicular to the laser beam axis is limited and can be neglected. This assumption allows the process to be considered taking place in one dimension.

Heating phase
The one-dimensional heating equation is where a is the thermal diffusivity of the material. Thermal diffusivity is a function of temperature; however, in this modelling approach, it will be considered being constant, for the simplicity of the model. At the beginning of every pulse, there is a short phase when the surface temperature has not reached the melting point, meaning that no phase transition is taking place. For the study of this phase, the appropriate initial condition for Eq. (3.1) is where T 0 (z) is the temperature field after the cooling period following the previous pulse. Typical values of process parameters used for laser drilling allow for a large cooling time comparable with the pulse duration. This means that the profile of the temperature after a cooling period will be flat in comparison with the temperature profile during the time period the pulse is on. Thus, it can be considered that the temperature field T 0 (z) is constant. This approximation has to be neglected in case that the pulse duration is a considerable percentage of the pulsing period or if continuous laser drilling is considered. There are two basic mechanisms of heat transfer inwards and outwards the workpiece, conduction and convection. As this study focuses on the processing of metals, convection is a negligible factor. However, if laser treatment of other materials is considered, convection should probably be taken into account. Thus, the boundary conditions for the heat equation are purely defined by the heat flux of the laser pulses where k is the thermal conductivity of the material and z 0 is the depth of the workpiece surface when the considered pulse begins, which may be different from zero if it is not the first pulse. Thermal conductivity may vary with temperature; however, since the acquired solution should be as simple as possible, it will be considered being constant.

Melting phase
When the surface temperature reaches the melting temperature, the two phases of the material, solid and liquid coexist and interchange heat. Then, in both regions, the heat equation is satisfied where the index l stands for liquid and the index s stands for solid. The two regions are separated by a planar surface at z = Z m (t). The temperature on this surface equals the melting temperature T m The velocity of the moving boundary between the two regions is set by the latent heat of fusion L f . The relevant boundary condition is the so-called Stefan condition, which is the energy conservation at the phase-separating surface.
Finally, the boundary condition is again given by Eq. (3.3) and the initial condition for the melting phase is provided by the outcome of the previous heating phase.

Vaporization phase
The vaporization phase begins when the surface temperature reaches the vaporization point. In the following approach, all materials that turn to a gaseous state is considered to be instantaneously escaping. Then, in the liquid and solid regions, the heat equation has to be satisfied as described by Eqs. (3.4) and (3.5).
The solid and liquid regions are separated by a planar surface at z = Z m (t), while the liquid and gaseous regions are separated by another planar surface z = Z v (t). The temperature at these surfaces equals the melting temperature T m and the vaporization temperature T v respectively, resulting in condition (3.6) and Similarly to the melting phase, the velocity of the moving boundaries among the three regions is set by the latent heat of fusion L f and the latent heat of vaporization L v . The relevant boundary conditions are Eq. (3.7) and Finally, the initial condition for the vaporization phase is provided by the outcome of the previous melting phase.

Cooling phase
When the pulse ends, a cooling phase follows, until the next pulse begins. The treatment of the cooling phase is identical to that of the heating and melting phases, in the sense that at the beginning of the cooling phase, the liquid and solid phases coexist and interchange energy. The melting phase will cease to exist very fast and cooling will continue in the presence of merely the solid phase. The only difference in the treatment is that the surface separating the liquid and the gaseous states stops moving and the boundary condition at this surface becomes: since there is no influx of energy during the cooling period. The appropriate initial condition is given by the outcome of the previous vaporization phase.

The numerical model
The finite differences method has been used to model laser drilling numerically. A restriction imposed by this approach is that the spatial dimension of the problem cannot be considered infinite. For the analytical approach, it is easier to consider the spatial dimension as infinite. However, if the length of the spatial dimension in the numerical approach is larger than the depth of the region affected by the laser, any reflection phenomena will be limited and therefore the difference between the two approaches may not be large. In this study, the discrete depths are considered uniformly spaced in the interval 0 ≤ z ≤ D, where D is the depth of the workpiece. There are nodes at z = 0 and z = D, thus z i = (i − 1) z, where i takes values from 1 to N z and z = D/ (N z − 1). Similarly, the discrete times are also considered uniformly spaced, thus t j = (j − 1) t, where j takes values from 1 to N t and t = t max / (N t − 1). The temperature field is a finite set of the elements T ij = T z i , t j .
A discrete version of the heat equation can be expressed and solved time slice by time slice. However, such an approach has the disadvantage that phase transitions have to be dealt manually, as described by Eqs. (3.7) and (3.9). Additionally, the heat equation has to be solved separately in the region of each phase.
There is an alternative way of bypassing this obstacle in a numerical method. The second spatial derivative that appears in the heat equation has the meaning of the infinitesimal difference of the incoming and outgoing heat currents, in an infinitesimal region. Hence, the second spatial derivative (multiplied with the thermal conductivity to recover the heat current) is the amount of heat rate that settles in this infinitesimal element. This is proportional (with constant of proportionality equal to the density times the specific heat) to the temperature rate of change at this infinitesimal element, resulting in the heat equation. The non-smoothness of the temperature field at a phase-separating surface is due to the fact that during a phase transition, the incoming heat does not correspond to an increase of temperature because of the existence of the latent heat. This phenomenon can be dealt more easily with the introduction of a new field corresponding to the energy density. Then, the second derivative of the temperature will result in a rate of change for the energy density and then the energy density will be connected to the local temperature by taking into account not only the specific heat but also the latent heat of fusion and evaporation. In other words, the heat equation, including phase transitions, can be expressed as where energy density and temperature are connected through where U m and U v are the required energy densities for melting and vaporization respectively. They equal The relation between energy density and temperature is indicated in Fig. 2. It is clear from Eqs. (3.11) and (3.12) that since the U i1 and T i1 are given in the form of the initial conditions the problem can be integrated time slice by time slice. Equation (3.11) allows the calculation of the energy density values in the next step, which can be translated into the temperatures in the next step using Eq. (3.12) and so on. The last detail deals with the fact that the definition of the second spatial derivative is problematic at i = 1 and i = N z or alternatively the way the boundary conditions are implemented need to be determined. A relevant concern is the fact that when evaporation starts, the boundary of the region studied is moving since material is removed. Thus, at each time slice, the position of the left boundary has to be specified. Evaporation is occurring when energy density is exceeding the critical value U v , given by Eq. (3.15). Thus, at each time slice, the boundary has to be specified as the minimum i for which the energy density is smaller than the critical value Finally, the boundary condition, described by Eq. (3.3) and a similar one for the opposite boundary describing vanishing heat flow can be expressed as where I a z b j , t j is given by Eq. (2.6). The problem is completely described by Eqs.  (3.20), which can be integrated time slice by time slice, providing the temperature field as a function of depth and time. This method can serve as a standalone numerical model for evaporation laser drilling. In this study, it also serves as a benchmark for the validity of the approximations used for the development of the analytical model in Section 3.3.

The analytical model
The phase transitions that occur during the specific process make the mathematical nature of the problem nonlinear. The solutions of the heat equation are characterized by an exponential decay towards a temperature field with no curvature; if there is no curvature of the temperature field, then it is also time independent. The approach followed was based on the assumption that the absorbed power density was so high that it locally "stretched" the temperature field sufficiently enough to assume it had a linear profile with depth. Once, such an assumption is made, the partial differential equation is transformable into an ordinary one for the parameters of the linear "ansatz" plugged into the equations.

Heating phase
According to the above-mentioned, in the current approximation, the temperature field during the heating phase is considered being equal to This assumption can be visualized as the red curve in Fig. 3.

Continuity demands that
It is assumed that the inclination of the aforementioned linear function in the heated zone z < Z max (t) is set by the heat flow boundary condition (3.3). Thus, (3.23) heating phase ansatz melting phase ansatz evaporation phase ansatz Fig. 3 The approximation of the temperature field during the heating, melting and vaporization phase Finally, there is the constraint of the overall energy conservation. The energy that has been transferred into the material has resulted in a temperature increase. This energy is simply cρ times the surface under the red curve in Fig. 3, where c is the specific heat. Thus, the energy conservation per area can be written as Equations (3.22), (3.23) and (3.24) are sufficient for determining the time evolution of the three unknown functions of the problem, namely Z max (t), A s (t) and B s (t).

Melting phase
During the melting phase, like during the heating phase, it is assumed that the temperature field is a linear function of depth. However, it may be a different linear function in the liquid region from that in the solid region. So, the temperature field is given by This assumption can be visualized as the green curve in Fig. 3. This "ansatz" contains six unknown functions of time. However, the temperature field has to be continuous in space. The temperature has to be equal to the melting temperature at the phase-separating surface z = Z m (t) and the environment temperature at z = Z max (t). These demands result in the following constraints . (3.28) After the continuity constraints, the problem has three unknown parameters that depend on time, the temperature inclination in the liquid region Al(t), the depth of the phaseseparating surface z = Z m (t) and the maximum depth that has been affected by the heat flow z = Z max (t). So, three equations are needed in order to specify the above parameters. The first of those is the same as the one used in the heating phase. The inclination of the aforementioned linear function in the liquid region, which is directly hit by the laser beam, is set by the heat flow boundary condition, −kA l (t) = I a . (3.29) The second equation is the Stefan condition (3.7), which in the linear "ansatz" is expressed as Finally, the third relation is simply the overall energy conservation. The energy that has been transferred into the material has been spent in two ways. The first way is the phase change latent heat, which is L f ρ multiplied by the volume that has changed phase. The second one is the energy corresponding to the temperature change, which is cρ multiplied by the surface under the green curve of Fig. 3. Thus, energy conservation (per area) can be written as Out of the three equations describing the problem, only the second one is a differential equation. Thus, Eqs. (3.29) and (3.31) can be solved algebraically for Z max (t) and B l (t) and the latter can be substituted in Eq. (3.30) to get an easily solvable first-order differential equation for Z m (t). Then, the other two parameters of the solution can be specified by Z m (t) with the use of Eqs. (3.29) and (3.31). The aforementioned differential equation can be written as .
The appropriate initial condition for this differential equation is a vanishing Z m at the time melting starts.

Vaporization phase
In the vaporization phase, similarly to the previous two phases, it is assumed that the temperature field is a linear function of depth; however, it may be a different linear function in the liquid region from that in the solid region. The temperature field is given by This assumption is depicted as the blue curve in Fig. 3. This "ansatz" contains seven unknown functions of time. The temperature field has to be continuous in space. The temperature equals the vaporization point at z = Z v (t), the melting point at the phase-separating surface z = Z m (t) and the environment temperature at z = Z max (t). These demands result in the following constraints . (3.38) After the continuity constraints, the problem has three unknown parameters that depend on time, the depth of the vapour-liquid separating surface z = Z v (t), the depth of the liquid-solid separating surface z = Z m (t) and the maximum depth that has been affected by the heat flow z = Z max (t). Three equations are required to specify the above parameters. Two of these are the Stefan conditions at the two phase-separating surfaces namely Eq. (3.30) and The third relation is simply the overall energy conservation. The energy that has been transferred into the material has been spent in three ways. The first is the vaporization latent heat, which equals L v ρ times the volume that has been vaporized. The second is the melting latent heat, which is L f ρ times the volume that has been melted. Finally, the third is the energy corresponding to the temperature change. This is simply cρ times the surface under the blue curve in Fig. 3, where c is the specific heat. Thus, energy conservation (per area) can be written as Out of the three equations regarding the problem, only two are differential equations. Therefore, Eq. (3.40) can be solved algebraically for Z max (t) and the latter can be substituted into Eqs. (3.30) and (3.39) to get an easily solvable first-order system of two differential equations for Z m (t) and Z v (t) where The appropriate initial conditions for this system of differential equations are a vanishing Z v at the time that evaporation starts and the outcome of the melting phases Eq. (3.32) for Z m .

Cooling phase
As described in Section 3.1.4, at the beginning of the cooling phase, the solid and liquid phases coexist. Soon enough, the liquid phase transmits adequate amount of heat to the solid phase and ceases to exist. Then, heat diffusion continues in the solid phase, homogenizing the temperature distribution into the workpiece, resulting in the cooling of the surface. As the numerical model indicates, the coexistence of the two phases lasts for a very short period of time thus, it is negligible. It is desirable to find an appropriate form for the temperature field during the cooling phase in order to study the cooling of the surface, using the same method as the one applied to the heating, melting and evaporation phases. Since only the solid state is involved, a natural guess would be the "ansatz" is used for the study of the heating phase. However, the boundary condition is now different, not allowing the use of a linear profile for the temperature field.
The physics of the problem and the numerical model solution indicate that a good candidate for the temperature field profile is a Gaussian profile of the form T (z, t) = C (t) e − (t )z 2 + T 0 . (3.44) Energy conservation relates the two unknown functions of time C (t) and (t). The integral of the temperature field in the whole z semi-axis times the density times the specific heat should equal the amount of energy left into the workpiece after the removal of the evaporated material and the end of the evaporation phase. The latter equals where Z v , Z m and Z max are calculated at the time instance the pulse ends. So, conservation of energy results in the following relation between C (t) and (t) Substituting the temperature profile in the heat equation and using Eq. (3.46) to eliminate C (t) results in the following equation which stands if The solution to the latter is simply , (3.49) where t 0 is the integral constant to be specified by condition (3.46) and the demand that C (0) = T v − T 0 .

The first pulse solution
The first step towards the solution of the model developed is the acquisition of the solution for the first pulsing period. In this case, the initial condition is the temperature field being equal to the environment temperature. The results for the first pulse will be the basis for the description of the advancement of drilling during many pulses in Section 3.5.
In the following graphs, the solutions of the analytical and numerical models are presented and a comparison is made between the two. For these graphs, the material of the workpiece is assumed to be described by the thermal properties given in Table 1. These numbers correspond to the properties of mild steels [38]. For the laser source, the variables given in Table 2 are used. These correspond to values characterizing CO 2 laser sources.

Heating phase
The solution of the analytical model for the heating phase is obtained by solving Eqs.   Thus, using the fact that diffusivity, conductivity, specific heat and density are connected as a = k/cρ, the solution for the heating phase (3.21) can be written as (3.52) A correct prediction of the approximation used is the fact that surface temperature increases with time like √ t [39]. The time that the surface temperature reaches the melting point equals (3.53) In Fig. 4, the temperature field during the heating phase of the first pulse is shown. In the displayed diagram, time stops when the melting point is reached. For Fig. 4, the material variables and the laser variables given in Tables 1  and 2 were used. It is noticeable that the time necessary to reach the melting point is very short, for the parameters used; it equals t m = 0.016 μs.

Melting Phase
In order to acquire the analytical model solution for the melting phase, Eq. (3.32) has to be integrated, so that the depth of the phase-separating surface as a function of time is specified. Once this is done, the full temperature field solution can be found using Eqs. (3.29) and (3.31). With the variables of Tables 1 and 2, the solution for the temperature field has been calculated and is shown in Fig. 5. It can Fig. 4 The temperature field solution for the heating phase The temperature gradient in the liquid phase remains constant because of the boundary condition. For the same reason, the temperature gradient in the solid phase also remains constant before the melting starts. However, as the melting starts, a maximum deviation of the latter from the value given by the boundary conditions is acquired very fast, and then it slowly converges towards this value. This effect is also visible in Fig. 6, which shows the depth of the phase-separating surface and the depth of the heated zone as functions of time. As Eq. (3.7) describes, the difference between the two inclinations is proportional to the rate that the boundary between the liquid and the solid regions moves. In conclusion, this velocity acquires a maximum value very fast once the melting starts and then it slows down.
When melting starts, there appears an extra obstacle in the heat conduction, the latent heat of fusion, resulting in an initially vanishing velocity of the phase-separating surface. Once the obstacle is surpassed, and there is a pool of liquid material on top of the solid phase, the advancement of the front of the phase transition increases, acquiring a maximum rate and then it slows down as it moves away from the heat source. This slowing down occurs because the phase-separating surface is getting further away from the workpiece surface, thus making heat transfer more and more difficult, as more time is required for the heat to travel this distance and additionally, because there is more material

Vaporization phase
In order for the solution of the temperature field, during the vaporization phase to be acquired, the differential system, Eqs. (3.41) and (3.42) has to be solved. In Figs. 7, 8, 9 and 10, the solution of these equations is presented, for the parameters given in Tables 1 and 2. The analytical model is in good agreement with the numerical model.
The temperature field for the whole pulse duration Fig. 9 The temperature field during heating, melting and the beginning of vaporization phases; the solid lines correspond to the two phase-separating surfaces and the boundary of the heated zone; the dashed lines separate the heating, melting and evaporation phases Unlike the melting phase, the temperature gradient in the liquid phase does not remain constant. This is clearly visible in Fig. 7. This is happening due to the latent heat of evaporation. Alike the melting phase, in the beginning of the vaporization phase, the velocity of the gas-liquid phase-separating surface is low; since at the time evaporation starts, an extra obstacle in heat conduction that of the latent heat of evaporation appears. This is clearly visible in Figs. 11 and 12, where the depth of the phase-separating surfaces and the heated zone as function of time are plotted.
Most importantly, after some time, the phase-separating surfaces acquire a constant velocity, the same for all of them

Cooling phase
The analytical solution for the cooling phase has already been provided in Section 3.3.4. The results of the analytical and numerical models for the parameters given in Tables 1  and 2 are shown in Fig. 13. The surface temperature is shown in Fig. 14. The numerical and analytical models are in good agreement, especially for time scales larger than that of the pulse duration. Since the typical process variables used in laser drilling allow for a cooling period much larger than the pulse duration, the analytical model presented here is a good approximation for the time evolution of the temperature field during the cooling period.
As it can be seen in Fig. 14, the temperature reaches that of the environment very fast. This, combined with the fact that the heating phase does not actually last long, as shown in Section 3.4.1, results in the remaining heat from one pulse being insignificant for the next one. Thus, it is a reasonable approximation to consider that the initial condition for any pulse is the temperature field being equal to the environment temperature.

Time evolution of drilling
In previous sections, the time evolution of the erosion front as a function of time has been modelled. An approximation of the differential equation solutions is going to be  Figure 12 suggests that both phase-separating surfaces and the depth of the heated zone can be very well approximated by three parallel straight lines, which is expected from energy conservation arguments, and also The linear approximation for the advancement of the erosion front during one pulse is a good approximation when the pulse duration is sufficient. However, as drilling advances, the laser defocusing will reduce the incident power density, resulting in lower drilling speed, meaning that evaporation does not really start early and the linear  Fig. 15 The depth of the phase-separating surface and the heated zone as function of time and the approximations approximation may not be sufficient. Thus, the linear approximation requires to be improved as it concerns the temporally advanced stages of laser drilling. The aforementioned will be analytically developed so as for the used function to reach the asymptotic linear form, for times larger than, z v /v for the evaporation front and be simultaneously enforced to vanish at the time that evaporation starts. Such an appropriate form is , (3.61) where t v and the time instant evaporation starts and v and z v are given by Eqs. (3.57) and (3.58). The evaporation starting time t v can be approximated, assuming that the inclination in both solid and liquid regions is specified by the boundary conditions during the melting phase. In this case, it is given by (3.62) In Fig. 16, it can be seen that the approximation (3.61) is much better than the linear one described by Eq. (3.54), especially for short times.
Since the evaporated depth is well approximated by Eq. (3.61), the advancement of drilling with time can be studied with the help of this equation. A key point of the process is the fact that the absorbed laser intensity is decreasing as drilling advances due to laser defocusing as described by Eq. (2.5). The constants v and z v depend on the absorbed intensity, as indicated in Eqs. (3.57) and (3.58). Thus, the drilling rate will decrease to reach asymptotically zero at a finite depth.
In the following approximation, the depth that is drilled in one pulse equals to: , (3.63) where v and z v are considered to depend on the drilled depth as explained above and Z v (t) is given by Eq.  16 The approximation for the drilling front there is no more drilling during the cooling phase, this is the advancement of drilling occurring in the entire pulsing period τ . Thus, the drilling process can be described by the differential equation The solution of the differential equation above, for the parameters given in Tables 1 and 2, is presented in Fig. 17.
The drilling rate decreases as drilling advances, going asymptotically to zero at a finite depth. As long as drilling has not advanced too much, so that the linear approximation of Fig. 15 is accurate, it can be said that the drilling rate is proportional to the absorbed intensity, thus formula (2.5) also describes the decrease of the drilling rate with depth.

Maximum drill depth
As indicated in Fig. 17, the drilling depth reaches an asymptotic value. This value could be calculated with Eq. (3.62) as being the depth for which the time for the starting of evaporation equals the pulse duration. The minimum absorbed energy density resulting in evaporation is An appropriate definition for the energy efficiency of a process is the ratio of an appropriate "result" and the energy required in order for this result to be attained [5]. The "result" can be selected as a quantity with energy units, typically the theoretical minimum energy required for the same "result" [40][41][42][43], and then, the energy efficiency is defined as being dimensionless. This kind of definition, although it can be implemented for a specific single process, it is difficult to be generalized for the entire class of processes that can provide the same "result", given that the theoretical minimum energy for the required "result" may be different for each kind of process. Thus, such a definition may not be appropriate for an energy efficiency that has to be used as a comparative measure among different processes. On the contrary, if the "result" is quantified by some geometric characteristic, these disadvantages do not appear, and additionally, an easy connection between the geometry described in a CAD file and the prediction of energy consumption can be provided by an IT system. For these reasons, the energy efficiency definition used in this study is given by the following relation [5]: where V removed is the removed volume and E is the energy required to remove this volume. In the following approach, the removed volume can be predicted as a function of the processing time, using the results of Section 3. The required energy can be simply calculated as where P is the laser beam power and t is the processing time.

Dependence of energy efficiency on process parameters
Since the speed of the erosion front is decelerating as shown in Fig. 17, the material removal rate also decreases. On the other hand, the energy consumption rate is constant; thus, the energy efficiency decreases as the process time increases. However, the most energy efficient selection of parameters for the opening of a hole with a specific depth has still to be determined. The previous argument suggests that the most energy efficient selection of parameters is the one that minimizes the process time too. In the following, three parameters will be considered, namely laser power, pulsing frequency and duty percentage; the latter being the ratio between the pulse duration and the pulsing period. The dependence of energy efficiency on the aforementioned process parameters is shown in Figs. 18, 19 and 20 for desired hole depth equal to 1, 2 and 4 mm. The parameters used for these graphs are given in Table 3, with the exception of the corresponding horizontal axis variable for each graph.
The trend of the plotted curves can be explained as follows: the increase in laser power as well as the decrease in duty results in higher laser beam intensity, thus in faster drilling and less time for heat lost due to conduction. The laser beam intensity does not depend on the pulsing frequency; however, reduction in the pulsing frequency simply results in less pulses being required for the process. In Section 3, it was shown that a specific part of each pulse, that depended on the laser beam intensity was lost due to conduction as described by Eq. (3.45). Consequently, lower pulsing frequency results in fewer losses and higher energy efficiency.
The selection of the process variables leading to higher energy efficiency is made clear in Figs. 18, 19 and 20. The most efficient laser power and duty selection correspond to higher laser intensities. It is well known [44,45] that intensities higher than 10 8 W/cm 2 result in multiphotonic and thermionic emission processes, both responsible for the primary electron generation from the irradiated metal target (plasma formation). When plasma is generated, the energy provided by the source is not consumed for processing purposes, but for plasma conservation. As a result, the energy efficiency of the process will decrease. Consequently, it is The dependence of energy efficiency on duty expected that the energy efficiency will reach a maximum value, as the laser intensity increases up to the point of decrease. The use of lower pulsing frequencies can prevent the generation of plasma [46]. Since lower pulsing frequencies prevent plasma formation and correspond to higher energy efficiencies, aiming at the increase of the latter, regarding the pulsed laser drilling process, it is the lowest available pulsing frequency of a machine that has to be applied. Then, the highest possible power and lowest possible duty that do not lead to plasma formation should be used. Plasma generation has not been considered in the current modelling approach; however, the process variables used for Figs. 18, 19 and 20 correspond to laser beam intensities smaller or marginally equal to the limit of 10 8 W/cm 2 , implying that there is negligible plasma formation.
As indicated in Figs. 18, 19 and 20, the energy efficiency takes a wide range of values, depending on the selection of process variables. Thus, proper adjustment of the latter can provide a significant improvement on energy efficiency, resulting in the corresponding gains in energy cost. As shown in Fig. 18, the energy efficiency for the drilling of a 4mm hole is doubled, when the laser power is being increased from 1.5 to 4 kW.

Physical analysis of energy efficiency
If an imaginary perfect machine is considered, driving all its energy flow towards heating only the material to be removed, from room temperature to evaporation temperature and then providing the latent heat in order for this Table 3 The laser variables used for the energy efficiency graphs P 2 kW f 1 kHz d 10 % material to be evaporated and disposed, its energy efficiency would be For the parameters used in this study, this corresponds to ideal energy efficiency about equal to Eef ideal = 16 mm 3 /kJ. In Figs. 18, 19 and 20, it is shown that the best selection of process variables corresponds to energy efficiency equal to about Eef best = 1.1 mm 3 /kJ. Thus, the best achieved energy efficiency is about 6.5 % of the theoretical ideal. Reflectivity is responsible for the loss of a percentage of the in-falling energy equal to Q r = R. (4.4) The rest of the losses are due to heat being diffused into the material and to the laser defocusing. Other factors may certainly intervene. However, these are the most important ones and those considered in the present study. The percentage of losses because of laser defocusing can be estimated as the complementary of the ratio of the volume of the cylinder with radius equal to the minimum beam radius, over the volume of the hourglass-shaped solid that is described by Eq. (4.5) The percentage of losses due to conduction is indirectly calculated in Section 3 since the energy efficiency calculated based of the results of this section is equal to (4.6) In Fig. 21, it is shown how energy efficiency is determined by these factors, as drilling advances. The parameters for this graph are the ones given in Table 3. It is evident that the most important factor determining the laser drilling energy efficiency, at process level, is the materials reflectivity.

Conclusions
In the current study, an analytical and numerical approach of evaporation laser drilling, towards the theoretical specification of the process parameters that maximize energy efficiency, has been developed. The process has been described as a function of the processing time, a critical factor, for the retrieval of a prediction of energy consumption. Based on these models, the energy efficiency and its dependence on the process parameters have been investigated. The investigation has revealed that the material removal rate decreased as the drilling/processing time increased, reaching a zero value at a given depth. The maximum drill depth was laser parameter dependent, increasing with the increase in laser power and/or duty. Furthermore, the energy efficiency of the pulsed laser drilling process was also found dependent on the laser source parameters, indicating an increase with the increase of the laser power and maintaining this trend with the decrease of the pulsing frequency and duty. It can be stated that there is enough evidence to prove that the energy losses during the laser drilling process, determining the energy efficiency, are occurring mainly due the material reflectivity, laser beam defocusing and heat conduction, with the most detrimental being reflectivity.
A future study will focus on the adaptation of the analytical model developed with experimental data. Improvements on the model are possibly necessary in order to include other factors, such as plasma formation, variations of reflectivity due to the formed geometry or heat losses due to convection. Another improvement of the analytical model that can be developed is the substitution of the last linear sector of the solution "ansatz" with an exponential one as it is suggested by the results of the numerical model. The introduced model can also be used for the study of the drilled hole's shape and for providing guidelines as to the way that effects such as barreling can be prevented.