Direct particle–fluid simulation of flushing flow in electrical discharge machining

Abstract

The efficient removal of material debris by flushing is critical for the performance and the surface quality of the electrical discharge machining (EDM) manufacturing process. The particle concentration alters the thermal, mechanical, and electrical properties of the particle–dielectricum suspension and affects the manufacturing process and the surface quality by inducing high thermal loads. For the first time, we perform a direct particle–fluid simulation of the flushing cycle in a generic EDM cavity using a hierarchical Cartesian sharp-interface cut-cell method. The flow around each debris particle is completely resolved and the heat transfer between the particles and the fluid is taken into account. The rate of material removal and the cooling performance of the flushing mechanism are studied for three-volume loadings. The results show that the flow around the particles has a pronounced impact on the heat transfer between the dielectricum and the workpiece. A particle loading of 14% increases the mean heat transfer by approximately 16%. The particle loading also has a pronounced impact on the rate at which particles are flushed out of the cavity. Increasing the initial volume loading from 6% to 14% decreases the amount of particles that get flushed out of the cavity by 14%.

Keywords:

• Direct particle–fluid simulation
• cut-cell method
• electrical discharge machining
• mono-hole inlet flushing
• multiphase flow

1. Introduction

Electrical discharge machining (EDM) (Singh et al., 2004) is a manufacturing process frequently used in the automotive, aerospace, microelectromechanical systems (Michel et al., 2000; Uhlmann et al., 1999), and bio-medical industry (Ho & Newman, 2003). It produces high-quality surfaces while allowing precise material removal. The process shapes a workpiece using a large number of high-frequency electrical discharges produced by an electrode submerged in a liquid insulating dielectricum under the action of an electric current. After removal, the material solidifies in the form of very small spherical hollow particles with a diameter in the range of $\mathcal{O}\left[10-100\right]\text{μ}\mathrm{m}$ depending on the process parameters (Haas et al., 2013; Murray et al., 2016). The particle concentration alters the thermal, mechanical, and electrical properties of the particle–dielectricum suspension and affects the manufacturing process and the resulting surface quality by inducing high thermal loads (Kumar et al., 2009). The heat transfer between the workpiece and the fluid is responsible for the temperature distribution in the workpiece's rim zone and critical for the resulting surface integrity. That is, due to thermal expansion effects, the heat flux defines the machining precision and is to be determined in detail to evaluate the EDM process. This study focuses on die-sink EDM where the material is removed in a small cavity within the workpiece. The stability of the process is dominated by the capacity of the flushing mechanism to cool the electrode and the workpiece and to pump the removed material through the lateral gaps of the electrode out of the cavity (Goodlet & Koshy, 2015; Uhlmann et al., 2005). Therefore, the die-sink EDM process uses a hybrid cycle with two phases: (I) a material removal phase in which the workpiece is eroded by a large number of high-frequency sparks and (II) a flushing phase in which the particles generated from the eroded material are flushed out from the working cavity. In the pressure flushing technique, the fresh dielectricum is pumped through a drilled machining electrode right into the cavity, to remove the material symmetrically through the lateral gap between the electrode and the workpiece. The lateral gap thicknesses between the electrode and the workpiece and the cavity depth are on the order of $\mathcal{O}\left(1000\right)\text{μ}\mathrm{m}$. These small gaps combined with the geometrical constraints of the problem and the physical properties of the EDM process, i.e. intensive electrical discharges with small time scales, make measurements of the particle–laden flow in the working cavity extremely challenging. Therefore, primarily numerical simulations were conducted to study the flushing mechanisms in EDM. Haas et al. (2013) performed numerical simulations of the fluid flow in wire-cut EDM using the Reynolds-Averaged Navier–Stokes (RANS) equations of a two-dimensional single-phase fluid to optimize the design of flushing nozzles. The analysis resulted in significant improvements of the cooling efficiency and wire stability of the process. The fluid flow induced by single electrical discharge was simulated by Shabgard et al. (2013) using a RANS method to study how the pulse time affects the flushing efficiency of molten material. Z. Wang et al. (2018) used steady two-dimensional axisymmetric simulations to develop a flushing mechanism that combines electrode-internal with electrode-external flushing and significantly improves the process efficiency by reducing the material debris concentration in the machining gap. Kliuev et al. (2018) performed three-dimensional RANS simulations using the k-ϵ turbulence model to evaluate the flushing efficiency of multi-channel electrode configurations in deep hole drilling.

Two-phase RANS simulations were performed by Xu et al. (2015) and Ji et al. (2013) to evaluate the performance of new flushing mechanisms. Xu et al. (2015) designed a multi-hole electrode with variable hole diameter, where the numerical simulations were used to optimize the diameter distribution and evaluate the performance of the electrode with respect to the produced material removal rate. Ji et al. (2013) simulated a new hybrid coaxial-orbital flushing technique and the material removal rate performance was evaluated against that of non-hybrid coaxial and orbital methods. Liu et al. (2018) studied the impact of ultrasonic electrode vibrations on the flushing of material debris from the machining gap using three-dimensional fluid simulations with a Lagrangian-particle model. The influence of the gas-bubble dynamics on the material removal and its role on the cavity flow field of the EDM process were analyzed by J. Wang and Han (2014) using a three-phase RANS approach. Recently, the heat transfer and the heat dissipation rate of the dielectric fluid were investigated by Das et al. (2019) using a three-dimensional fluid simulation with a state-of-the-art bubble-dynamics model that takes material vaporization into account. The erosion, ejection, and solidification processes in EDM were simulated by Zhang et al. (2019). Modeling the heat flux density, vapor flow, and the stages of the eroded material prior to particles, i.e. melting, ejection, and solidification, the EDM process of a single–pulse discharge experiment could be analyzed.

Note that these studies and numerous other analyses of the flushing cycle found in the literature are based on the RANS approach. That is, the accuracy of the results relies on turbulence models developed for aerodynamic slender body flows over smooth surfaces. However, in the analysis of EDM flushing impinging jets, boundary layer separation, and complex vortical structures occur. Furthermore, the application of one-way coupling between fluid and particles is an extremely simple model since it does not account for the interaction between the phases (Elghobashi, 1994). Even the Euler–Lagrange approach for multi-phase models determines an uncertainty since a fundamental condition that has to be satisfied by the numerical method for such analyses is the conservation of mass, momentum, and energy at the particle–fluid interface. Euler-Lagrangian techniques are generally not well suited, since to satisfy the geometric conservation law expensive and continuous remeshing is required due to the arbitrarily large particle displacements which are not known a priori. Non-boundary conforming methods avoid remeshing by superimposing the particle geometry on a background mesh (Saliha et al., 2019). In the original immersed boundary method (Peskin, 1972), the effect of the moving boundaries on the flow is simulated by a discrete forcing term in the momentum equation. The known jumps across the interface are taken into account via correction terms in the numerical discretization (Z. Li & Lai, 2001). This approach, however, suffers from being not strictly conservative and smearing the sharp interface over several cell layers. A superior approach to overcome the limits of modeling are direct particle–fluid simulations (DPFS).

Direct particle–fluid simulations, i.e. numerical simulations of multi-phase flows in which turbulence and the flow field around each particle are fully resolved, are computationally a challenging multi-physics multi-scale problem when the particle diameter is in the range of the Kolmogorov length scale (Balachandar & Eaton, 2010). Schneiders et al. perform DPFS to investigate particle–laden isotropic turbulence for particle and turbulence parameters that forbid any modeling approach (Schneiders et al., 2019, 2017a, 2017b). They developed and used a massively parallel (Lintermann et al., 2014) adaptive hierarchical Cartesian cut-cell method (Hartmann et al., 2011) for multi-phase flows (Meinke et al., 2013; Schneiders et al., 2016) and complex geometries (Günther et al., 2014).

Based on the method by Schneiders et al. (2017a), we perform DPFS of the mono-hole pressure flushing flow in a generic die-sink EDM cavity, i.e.the flow around each debris particle is fully resolved and the heat transfer between the particles and the fluid is also taken into account. Hence, the effect of the particle-particle interaction in the dielectricum and the fluid–particle interactions are accounted for. We address the question whether fluid–particle interactions can be neglected or should be modeled in the numerical analyses of the EDM flushing cycle to predict the thermal loads that need to be sustained by the workpiece. The necessity of incorporating the debris particles in numerical simulations of the flushing cycle of EDM is studied by quantifying their impact on the total heat exchange.

The paper is organized as follows. In Section 2, the generic die-sink EDM cavity and the mathematical model to describe the multi-phase flow are introduced. In Section 3, the numerical method is presented. Then, the results are discussed in Section 4 and finally, some conclusions are drawn in Section 5.

2. Mono-hole pressure flushing in die-sink EDM

The generic axisymmetric die-sink EDM cavity is schematically shown in Figure 1(a), in which the cavity geometry and boundary conditions of the problem are described. A cut through the drilled electrode, depicted by dashed lines, reveals the single pipe-like hole of diameter $d_{\mathrm{i}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{t}}$ through which the fresh dielectricum enters the cavity. The cylindrical cavity has a diameter $d_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}$, a height $h_{\mathrm{g}\mathrm{a}\mathrm{p}}$, and is initially filled with a particle suspension at rest. The particle–laden fluid exits the system through the lateral cylindrical gap of thickness $s_{\mathrm{g}\mathrm{a}\mathrm{p}}$ between the electrode and the workpiece. On the bottom boundary between the fluid and the workpiece, a temperature distribution is used to model the latent heat stored in the workpiece during the material removal phase of the EDM process.

Figure 1. (a) Axisymmetric geometry of the die-sink EDM cavity including boundary notations. (b) Axisymmetric temperature distribution from thermographic camera measurements (Brito Gadeschi, Schneider, et al., 2017).

Next, the mathematical model to describe the two-phase fluid flow is introduced. This includes the model for the unsteady flow, the rigid body particle dynamics, the heat transfer between the particles and the fluid, and the boundary conditions. The parameters of the fluid and the solid particles are given in Tables 2 and 4 in Section 4.

2.1. Mathematical model

This section briefly introduces the mathematical model used to describe the physics of the two-phase solid–fluid flow in the EDM gap once all gas bubbles created by the electrical discharges have left the cavity, i.e. a third gaseous phase is not considered in this study. The Navier–Stokes equations are used to describe the motion of the fluid, whereas the debris particles are modeled as rigid spheres.

2.1.1. Navier–Stokes equations

The integral non-dimensional form of the unsteady conservation equations of mass, momentum, and energy for a moving control volume $V\left(t\right)$ reads (1) $\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\begin{array}{c}V\left(t\right)\end{array}}\mathbf{Q}\mathrm{d}V+{\oint }_{\mathrm{\partial }V\left(t\right)}\mathbf{F}\left(\mathbf{Q}\right)\cdot \mathbf{n}\mathrm{d}A=0.$(1) The conserved quantities $\mathbf{Q}=(\rho ,\rho {\mathbf{u}}^T,\rho E{)}^{\mathrm{T}}$ are defined by the density ρ, pressure p, velocity vector $\mathbf{u}$, and specific total energy E. The nonlinear flux vector (2) $\mathbf{F}\left(\mathbf{Q}\right)=\left(\begin{array}{c}\rho (\mathbf{u}-{\mathbf{u}}_{\mathrm{\partial }V})\\ \rho \mathbf{u}(\mathbf{u}-{\mathbf{u}}_{\mathrm{\partial }V})+p\\ \rho E(\mathbf{u}-{\mathbf{u}}_{\mathrm{\partial }V})+p\mathbf{u}\end{array}\right)+\frac{1}{R{e}_0}\left(\begin{array}{c}0\\ \overline{\mathit{\tau }}\\ \overline{\mathit{\tau }}u+\mathbf{q}\end{array}\right)$(2) Reynolds number $R{e}_0$, shear stresses $\mathit{\tau }$, and heat flux $\mathbf{q}$. The quantity ${\mathbf{u}}_{\mathrm{\partial }V}$ denotes the velocity of the control volume. Since the fluid can be considered to be incompressible, a low Mach number of M = 0.1 is assumed. Further information on the governing equations can be found in Brito Gadeshi, Schneiders, et al. (2017) and Schneiders et al. (2017a, 2017b).

2.1.2. Boundary conditions

The boundary conditions for the inflow, outflow, and solid wall boundaries are introduced next. Let Γ denote a boundary surface, the inflow boundary condition on ${\mathrm{\Gamma }}_{\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}}$ is given by (3) $\rho ={\rho }_{\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}},\mathbf{u}=\mathbf{n}\cdot u_{\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}}(r,t),\mathrm{a}\mathrm{n}\mathrm{d}\frac{\mathrm{\partial }p}{\mathrm{\partial }\mathbf{n}}=0,$(3) where $\mathbf{n}$ is the vector normal to the boundary surface pointing inside the fluid domain. The quantity $r\in [0,d_{\mathrm{i}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{t}}/2]$ is the radial coordinate of the cylindrical coordinate system shown in Figure 1(a). The axisymmetric parabolic laminar velocity profile $u_{\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}}(r,t)$ prescribed at the inflow reads (4) $u_{\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}}(r,t)=u_{max}\left(t\right)\left(1-\frac{4r^2}{d_{\mathrm{i}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{t}}^2}\right),$(4) where the maximum velocity at r = 0 is $u_{max}\left(t\right)$. The maximum velocity is $u_{max}\left(t\right)=u_{max}{\chi }_{[a,b]}\left(t\right)$ using the rectangular function χ that is $\chi =1$ for $a\le t\le b$ and $\chi =0$ otherwise. The mass flow is (5) $\dot{m}\left(t\right)={\rho }_{\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}}{\int }_0^{d_{\mathrm{i}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{t}}/2}{u}_{\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}}(r,t)2\pi r\mathrm{d}r.$(5) On the outflow boundary ${\mathrm{\Gamma }}_{\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}}$, (6) $\frac{\mathrm{\partial }\rho }{\mathrm{\partial }\mathbf{n}}=0\mathrm{a}\mathrm{n}\mathrm{d}\frac{\mathrm{\partial }\mathbf{u}}{\mathrm{\partial }\mathbf{n}}=0$(6) are prescribed. The pressure is defined by the hydrostatic pressure $p_{\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}}=p_{\mathrm{h}\mathrm{y}\mathrm{d}\mathrm{r}}$ at the outflow of the lateral cylindrical gap. On the bottom of the cavity ${\mathrm{\Gamma }}_{\mathrm{w}\mathrm{o}\mathrm{r}\mathrm{k}\mathrm{p}\mathrm{i}\mathrm{e}\mathrm{c}\mathrm{e}}$, a no-slip boundary condition with a prescribed temperature distribution (7) $\mathbf{u}=0,T=T_{\mathrm{s}}\left(r\right)\mathrm{a}\mathrm{n}\mathrm{d}\frac{\mathrm{\partial }p}{\mathrm{\partial }\mathbf{n}}=0$(7) is imposed, where the pressure is determined by the approximation of a vanishing normal gradient. The surface temperature $T_{\mathrm{s}}$ on the workpiece which is sketched in Figure 1(a) is obtained from thermographic camera measurements by Brito Gadeschi, Schneider, et al. (2017) the distribution of which is shown in Figure 1(b). It is prescribed using a steady Gaussian temperature profile with a peak temperature $T_{max}={42}^{\circ }\mathrm{C}$.

2.2. Heat-conducting moving solid particles

The debris particles are modeled as rigid spheres of diameter $d_{\mathrm{p}}$ and their boundary denoted by ${\mathrm{\Gamma }}_{\mathrm{p}}$ is illustrated in Figure 2(a). First, the governing equations are briefly discussed. Then, the boundary conditions are defined.

Figure 2. (a) Exemplary two-dimensional hierarchical Cartesian grid in which the cells intersected by the particle–fluid interface ${\mathrm{\Gamma }}_{\mathrm{p}}$ are reshaped to conform to the moving boundary. (b) Modeling of the collision between the particle i and a solid wall at distance $d_{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}}=d_{\mathrm{c}\mathrm{o}\mathrm{l}}/2$ using a mirror particle j such that the short-range repulsive force ${\mathbf{F}}_{\mathrm{c}\mathrm{o}\mathrm{l}}$ points in the wall-normal direction towards the particle i.

2.2.1. Governing equations

The governing equations of translational and rotational motion are (8) $\begin{array}{rl}& m\frac{\mathrm{d}{\mathbf{u}}_{\mathrm{p}}}{\mathrm{d}t}=\mathbf{F}+{\mathbf{F}}_{\mathrm{c}\mathrm{o}\mathrm{l}}+m\mathbf{g}\mathrm{a}\mathrm{n}\mathrm{d}\end{array}$(8) (9) $\begin{array}{rl}& \mathbf{I}\frac{\mathrm{d}\hat{\mathit{\omega }}}{\mathrm{d}t}+\hat{\mathit{\omega }}\times \left(\mathbf{I}\hat{\mathit{\omega }}\right)=\hat{\mathbf{T}}\end{array}$(9) with the particle mass m, the particle velocity ${\mathbf{u}}_{\mathrm{p}}$, the hydrodynamic force $\mathbf{F}$, the collision force ${\mathbf{F}}_{\mathrm{c}\mathrm{o}\mathrm{l}}$, the gravitational acceleration $\mathbf{g}$, the moment of inertia $\mathbf{I}$, the angular velocity $\hat{\mathit{\omega }}$, and the hydrodynamic torque $\hat{\mathbf{T}}$. That is, the particles move in response to the fluid forces and alter the flow field, which in turn changes the fluid forces acting on the particles. Hence, a fully coupled fluid-structure interaction problem is considered.

Furthermore, the particles exchange heat with the embedding fluid altering the particle temperature $T_{\mathrm{p}}$. The change of $T_{\mathrm{p}}$ is modeled using a first-order approach for small particles (10) $\frac{\mathrm{d}{T}_{\mathrm{p}}}{\mathrm{d}t}=\frac{-1}{{\rho }_{\mathrm{p}}{c}_{\mathrm{p}}{V}_{\mathrm{p}}}{\oint }_{{\mathrm{\Gamma }}_{\mathrm{p}}}{\mathbf{q}}_{\mathrm{p}}\cdot \mathbf{n}\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{p}}.$(10) The material parameters of the particle density ${\rho }_{\mathrm{p}}$, specific heat capacity $c_{\mathrm{p}}$, and thermal conductivity ${\lambda }_{\mathrm{p}}$ are listed in Table 4. Further computational details can be found in Brito Gadeshi, Schneiders, et al. (2017) and Schneiders et al. (2017a, 2017b).

2.2.2. Boundary conditions

The momentum and energy exchange between the particle and the fluid is prescribed using the following boundary condition for the fluid phase on the particle boundary ${\mathrm{\Gamma }}_{\mathrm{p}}$ (11) $\begin{array}{rl}\mathbf{u}& ={\mathbf{u}}_{{\mathrm{\Gamma }}_{\mathrm{p}}},\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathbf{n}}={\mathbf{q}}_{\mathrm{p}}\mathrm{a}\mathrm{n}\mathrm{d}\frac{\mathrm{\partial }p}{\mathrm{\partial }\mathbf{n}}\left(\frac{\gamma }{T_{\mathrm{f}}}\frac{D\mathbf{u}}{Dt}\cdot \mathbf{n}\right)\\ & =(\mathrm{\nabla }\cdot \mathit{\tau })\mathbf{n}.\end{array}$(11)

2.3. Initial condition

The flow field for a quiescent fluid is initialized by $\mathbf{u}=\mathbf{0}$, $T=T_{\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}}$, and $p=p_{\mathrm{h}\mathrm{y}\mathrm{d}\mathrm{r}}$. Initially, the particles are at rest and the particle temperature equals the maximum temperature obtained from the thermographic camera measurements by Brito Gadeschi, Schneider, et al. (2017). The particles are distributed homogeneously and axisymetrically at a distance of one $d_{\mathrm{p}}$ over the cavity bottom surface.

3. Numerical method

The cut-cell finite-volume method used in this work is a conservative multi-phase method for compressible flows (Hartmann et al., 2008, 2011) that has been extended to moving boundary problems (Schneiders et al., 2016, 2013), complex moving boundaries (Günther et al., 2014), and massively parallel computing (Lintermann et al., 2014). It has been successfully applied to study particle–laden flows using Euler–Lagrange (Siewert et al., 2014) and fully resolved particle models (Schneiders et al., 2019, 2017a, 2017b). In Brito Gadeschi et al. (2015) and Brito Gadeshi, Schneiders, et al. (2017), also the heat transfer between particles and fluid is taken into account. As for direct-numerical simulations (DNS) of single-phase flows, all scales of particle–laden flows are resolved in the DPFS, i.e. no turbulence modeling is used. A brief description of the numerical method, i.e. the spatial discretization, the temporal integration, and the representation of particles, is given in the remainder of this section.

The fluid equations described in Section 2.1 are solved on an unstructured hierarchical Cartesian mesh (Schneiders et al., 2016). For the spatial discretization of the inviscid flux terms an advection upstream splitting method (AUSM) is used (Liou & Steffen, 1993; Meinke et al., 2002). The viscous flux is discretized by central differences.

The boundary of the particles ${\mathrm{\Gamma }}_{\mathrm{p}}$ in Figure 2(a) is described by a level-set function on the hierarchical Cartesian grid (Osher & Sethian, 1988). The function is the scalar field of the shortest distance between point $\mathbf{x}$ and the particle surface ${\mathrm{\Gamma }}_{\mathrm{p}}$. The level-set field is computed analytically using the particle position ${\mathbf{x}}_{\mathrm{p}}$ and the radius $r_{\mathrm{p}}$ (Günther et al., 2014). The fluid cells intersected by the solid boundary are reshaped to yield a linear boundary representation. The pressure and stress forces are integrated over the boundary of the fully resolved particles and the right-hand side of Equations (8(8) $\begin{array}{rl}& m\frac{\mathrm{d}{\mathbf{u}}_{\mathrm{p}}}{\mathrm{d}t}=\mathbf{F}+{\mathbf{F}}_{\mathrm{c}\mathrm{o}\mathrm{l}}+m\mathbf{g}\mathrm{a}\mathrm{n}\mathrm{d}\end{array}$(8) ) and (9(9) $\begin{array}{rl}& \mathbf{I}\frac{\mathrm{d}\hat{\mathit{\omega }}}{\mathrm{d}t}+\hat{\mathit{\omega }}\times \left(\mathbf{I}\hat{\mathit{\omega }}\right)=\hat{\mathbf{T}}\end{array}$(9) ) are determined to move the particles. Analogously, the heat flux is computed to integrate Equation (10(10) $\frac{\mathrm{d}{T}_{\mathrm{p}}}{\mathrm{d}t}=\frac{-1}{{\rho }_{\mathrm{p}}{c}_{\mathrm{p}}{V}_{\mathrm{p}}}{\oint }_{{\mathrm{\Gamma }}_{\mathrm{p}}}{\mathbf{q}}_{\mathrm{p}}\cdot \mathbf{n}\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{p}}.$(10) ). The forces of particle–particle and particle–wall collisions ${\mathbf{F}}_{\mathrm{c}\mathrm{o}\mathrm{l}}$ in Equation (8(8) $\begin{array}{rl}& m\frac{\mathrm{d}{\mathbf{u}}_{\mathrm{p}}}{\mathrm{d}t}=\mathbf{F}+{\mathbf{F}}_{\mathrm{c}\mathrm{o}\mathrm{l}}+m\mathbf{g}\mathrm{a}\mathrm{n}\mathrm{d}\end{array}$(8) ) are based on a model by Glowinski et al. (1999) illustrated in Figure 2(b). For the latter case, a mirror image of the particle outside the physical domain is considered such that the direction of the short-range repulsive force is normal to the solid wall.

The heat flux ${\mathbf{q}}_{\mathrm{p}}$ in Equation (10(10) $\frac{\mathrm{d}{T}_{\mathrm{p}}}{\mathrm{d}t}=\frac{-1}{{\rho }_{\mathrm{p}}{c}_{\mathrm{p}}{V}_{\mathrm{p}}}{\oint }_{{\mathrm{\Gamma }}_{\mathrm{p}}}{\mathbf{q}}_{\mathrm{p}}\cdot \mathbf{n}\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{p}}.$(10) ) is based on Fourier's law and the temperature difference between the particle and the surface of the particle. While the particle temperature $T_{\mathrm{p}}$ represents an average, the temperature at the surface of the particle is defined by the local boundary conditions at ${\mathrm{\Gamma }}_{\mathrm{p}}$.

For the temporal integration of the flow field governed by Equation (1(1) $\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\begin{array}{c}V\left(t\right)\end{array}}\mathbf{Q}\mathrm{d}V+{\oint }_{\mathrm{\partial }V\left(t\right)}\mathbf{F}\left(\mathbf{Q}\right)\cdot \mathbf{n}\mathrm{d}A=0.$(1) ) and the particle motion described by Equations (8(8) $\begin{array}{rl}& m\frac{\mathrm{d}{\mathbf{u}}_{\mathrm{p}}}{\mathrm{d}t}=\mathbf{F}+{\mathbf{F}}_{\mathrm{c}\mathrm{o}\mathrm{l}}+m\mathbf{g}\mathrm{a}\mathrm{n}\mathrm{d}\end{array}$(8) ) and (9(9) $\begin{array}{rl}& \mathbf{I}\frac{\mathrm{d}\hat{\mathit{\omega }}}{\mathrm{d}t}+\hat{\mathit{\omega }}\times \left(\mathbf{I}\hat{\mathit{\omega }}\right)=\hat{\mathbf{T}}\end{array}$(9) ), a 5-stage Runge–Kutta scheme is used (Schneiders et al., 2016). The predictor–corrector scheme ensures stability up to a Courant–Friedrichs–Levy (CFL) number of 4. In this study, the CFL numbers for the inviscid and viscous terms are $C_f=1.0$ and $C_f^{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c}}=0.5$.

4. Results

The numerical method has previously been validated for multi-phase flows (Schneiders et al., 2016, 2017a) and conjugate heat transfer (Brito Gadeschi et al., 2014, 2015). In the following, first, the single-phase formulation of the numerical method is validated for an internal flow problem and the accuracy of the particle–wall collision model is assessed. Next, the heat transfer and particle transport of mono-hole pressure flushing is analyzed. The single-phase flushing flow is studied first, and then these results are compared against the particle–laden flushing flow.

4.1. Channel flow

The findings of the single-phase formulation of the numerical method for channel flow are compared with the direct numerical simulation results of Moser et al. (1999). The flow and numerical parameters are given in Table 1 and the simulation setup is shown in Figure 3. The quantity H is the channel height, $L=\pi H$ the channel length, and W = 1.75H the cannel width. Symmetry boundary conditions are used on the lateral sides of the channel and the inflow bulk velocity $\overline{u}$ and the outflow pressure $p_o$ are prescribed according to the pressure gradient required to obtain $R{e}_{\tau }=180$ (Moser et al., 1999). The Mach number is M = 0.1 and the simulated time equals the time required for a fluid particle to traverse the channel length 27 times according to the fluid bulk velocity. The number of cells is approximately $1.5\cdot {10}^6$ and the number of time steps is 200,000.

Figure 3. Domain and boundary conditions of channel flow.

Table 1. Flow and numerical parameters for channel flow.

Parameter	Value
Mach number: M	0.1
Prandtl number: Pr	0.72
Ratio of specific heats: γ	1.4
$R{e}_{H/2}$	2800
$R{e}_{\tau }$	180
Channel height: H	1
Length in the streamwise direction: L	$\pi H$
Length in the spanwise direction: W	1.75H
Number of cells: $n_c$	1,572,864
Number of time steps: $n_{\mathrm{\Delta }t}$	200,000
Simulated time: $\mathrm{\Delta }{t}_s$	$27L/\overline{\mathbf{u}}$

Table 2. Geometric, material, boundary and initial condition parameters for generic flushing-flow in die-sink EDM using mono-hole pressure flushing with the Oelheld IME-68 dielectric (Oelheld dielectrics, https://www.oelheld.com/).

Cavity geometric parameter	Value	Unit
Lateral gap: $s_{\mathrm{g}\mathrm{a}\mathrm{p}}$	0.4	[mm]
Electrode diameter: $d_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{d}\mathrm{e}}$	9.2	[mm]
Cavity diameter: $d_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}$	10	[mm]
Cavity height: $s_{\mathrm{h}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}$	1.35	[mm]
Material parameters at $T={20}^{\circ }\mathrm{C}$	Value	Unit
Fluid density: ${\rho }_{\mathrm{f}}$	0.77	[g cm${}^{-3}$]
Fluid kinematic viscosity: η	1.8	[mm${}^2$ s${}^{-1}$]
Fluid thermal conductivity: λ	0.1	[Wm${}^{-1}$ K${}^{-1}$]
Fluid specific heat capacity: $c_{\mathrm{p}}$	1.7	[kJ kg${}^{-1}$ K${}^{-1}$]
Flow parameters	Value	Unit
Prandtl number: Pr	23.4	[–]
Reynolds number: $R{e}_{\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}}$	94.5	[–]
Mach number: M	0.1	[–]
Boundary and initial conditions	Value	Unit
Max. inlet velocity: $u_{max}$	170	[mm s${}^{-1}$]
Inlet velocity period: $\mathrm{\Delta }{t}_{\mathrm{i}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{t}}$	0.5	[s]
Fluid inlet temperature: $T_{\mathrm{i}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{t}}$	20	[${}^{\circ }\mathrm{C}$]

The instantaneous ${\lambda }_2$-contours (Jeong & Hussain, 1995) colored by the streamwise velocity component in Figure 4 illustrate the turbulent structures of the flow field. The root-mean-square (RMS) of the velocity components spatially averaged in the streamline and spanwise direction are shown in Figure 5(a). The temporal averaging was computed over 27 convective units. A good agreement with the results from Moser et al. (1999) is obtained for $〈v^{\mathrm{\prime }}{〉}^{+}$, $〈w^{\mathrm{\prime }}{〉}^{+}$, and the turbulent shear stress $〈u^{\mathrm{\prime }}{v}^{\mathrm{\prime }}{〉}^{+}$. A slight overshoot of $〈u^{\mathrm{\prime }}{〉}^{+}$ in the range $y/H\in [0.09,0.3]$ is observed. The shear balance in Figure 5(b) shows a minor deviation of the total shear stress τ in the range $y/H\in [0.05,0.15]$. Finally, the comparison of the distributions of the streamwise velocity component $u^{+}$ normalized by the friction velocity in Figure 5(c) shows an excellent agreement with the data of Moser et al. (1999).

Figure 4. Instantaneous ${\lambda }_2$-contours colored by the non-dimensional streamwise velocity component $u/\overline{u}$.

Figure 5. (a) Root-mean-square (RMS) distributions of the velocity components; (b) shear balance, laminar and turbulent shear stress, ${\tau }_l$ and ${\tau }_t$, and the total shear stress $\tau ={\tau }_l+{\tau }_t$; (c) streamwise velocity component $u^{+}$ normalized by the friction velocity $u_{\tau }$ on a semi-logarithmic scale.

In conclusion, these results for the single-phase channel flow show a convincing match with the results of Moser et al. (1999). The overshoot in $〈u^{\mathrm{\prime }}{〉}^{+}$ and the deviations observed in the shear stress τ are small and can be attributed to the usage of different meshes between the current simulations and the reference. The findings demonstrate the validity of the method to accurately analyze internal flows such as EDM flushing.

4.2. Sphere–wall collision

The study of the dynamics of a sphere immersed in a viscous fluid impinging on a solid wall is an area of active research (Izard et al., 2014; X. Li et al., 2011). In this section, the results of the numerical method for particle-wall collisions based on Glowinski's model (Glowinski et al., 1999) are compared with findings from the literature. The experimental and numerical study of X. Li et al. (2011) considers a configuration, ‘Case 7’, with a ratio of the particle density and the fluid density ${\rho }_{\mathrm{p}}/{\rho }_{\mathrm{f}}=7780/1203=6.47$, which is on the same order of magnitude as the density ratio in the subsequent sections of this study (${\rho }_{\mathrm{p}}/{\rho }_{\mathrm{f}}=10.17$ ).

The simulation setup is shown in Figure 6(a). The stationary solid wall is located at $x_{\mathrm{w}}=0.0$ and free-stream boundary conditions are imposed on all other outer boundaries. The fluid is a mixture of Glycerol and water with a viscosity of $50.2\cdot {10}^3\mathrm{P}\mathrm{a}\mathrm{s}$. The time is non-dimensionalized by $\sqrt{d_{\mathrm{p}}/g}$. The particle is initially at rest and its bottom is located at a distance of 3.75 $d_{\mathrm{p}}$ above the wall and at rest. The impact Reynolds number is $Re\approx 127$. Note that the particle does not reach its terminal velocity before impact.

Figure 6. Impact of a sphere against a solid wall: (a) domain and boundary conditions; (b) computational grid; (c) comparison of the position of the particle bottom $x_{\mathrm{b}}$ as a function of time $t/\sqrt{d_{\mathrm{p}}/g}$ computed using the Glowinski model with the numerical results of X. Li et al. (2011).

The three-dimensional computational domain is discretized using the adaptive hierarchical Cartesian grid is shown in Figure  6(b). The grid contains approximately 800,000 cells and the particle diameter is discretized by approximately 25 cells, that is, the smallest cell length is $\mathrm{\Delta }{x}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}=0.04d_{\mathrm{p}}$.

Figure 6(c) shows the distance between the particle bottom and the wall as a function of time $x_{\mathrm{b}}\left(t\right)$. The time span in which the first two collisions with the wall occur is shown. The first collision occurs at t = 3.6 and the second one at $t=6.9$. During the initial fall, the Glowinski model results agree well with the experimentally calibrated numerical results of X. Li et al. (2011) and the time instant at which the first collision with the wall occurs is captured by the scheme. Due to the grid in conjunction with the discretization, the Glowinski model starts repulsing the particle before the actual collision occurs, i.e. the particle's minimum distance to the wall is $x_{\mathrm{b}}=0.26d_{\mathrm{p}}$ which corresponds to 6.5 grid cells. Furthermore, the Glowinski model transfers less energy into the particle resulting in a smaller peak height at t = 5.1 after the first bounce.

Note that the collision model used by X. Li et al. (2011) is computationally very costly and models the lubrication and elastic force terms with experimentally calibrated coefficients. For the purpose of this study considering EDM flushing, in which numerous particles and particle–particle and particle–wall collisions are considered, the simpler Glowinski model is a reasonable compromise between accuracy and computational efficiency.

4.3. Flushing flow

Next, the flow field in a generic die-sink EDM cavity using mono-hole pressure flushing described in Section 2 is studied. First, the geometric and flow parameters are introduced and an appropriate spatial discretization for the problem is defined. Then, the single-phase flow, i.e. the flushing flow without particles, is analyzed in Section 4.3.2. Subsequently, the particle–laden flushing flow is investigated and the single-phase and multi-phase flow results are compared in Section 4.3.3.

4.3.1. Geometry and flow parameters

The parameters governing the geometry of the die-sink EDM cavity illustrated in Figure 7(left) are the cavity height $s_{\mathrm{h}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}=1.3\mathrm{m}\mathrm{m}$ and the cavity and electrode diameters $d_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}=10\mathrm{m}\mathrm{m}$ and $d_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{d}\mathrm{e}}=9.2\mathrm{m}\mathrm{m}$. The difference $d_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}-d_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{d}\mathrm{e}}$ defines the lateral gap of the cavity $2\cdot s_{\mathrm{g}\mathrm{a}\mathrm{p}}=0.8\mathrm{m}\mathrm{m}$. The aspect ratios between the lateral cavity gaps and the cavity diameter and between the cavity height and the cavity diameter are $2s_{\mathrm{g}\mathrm{a}\mathrm{p}}:d_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}=1:12.5$ and $s_{\mathrm{h}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}:d_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}=1:7.7$.

Figure 7. (a) Radial axisymmetric slice of the hierarchical Cartesian grid $g_{10}$; (b) turbulent kinetic energy $E_k$ of the fluid phase over time for the grids $\{g_8,g_9,g_{10},g_{11}\}$.

The dielectric fluid Oelheld IME 68 is used. Its density and kinematic viscosity at $T={20}^{\circ }\mathrm{C}$ are ${\rho }_{\mathrm{f}}=0.77\mathrm{g}{\mathrm{c}\mathrm{m}}^{-3}$ and $\eta =1.8{\mathrm{m}\mathrm{m}}^2{\mathrm{s}}^{-1}$. The Prandtl number of the flow is Pr = 23.4 and the Mach number is M = 0.1. The Reynolds number of the flow is $R{e}_{\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}}=94.5$. A parabolic velocity profile with a maximum inlet velocity $u_{max}=170\mathrm{m}\mathrm{m}{\mathrm{s}}^{-1}$ is applied at the inflow boundary condition for a period of $\mathrm{\Delta }{t}_{\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}}=0.50\mathrm{s}$, i.e. before $t=0\mathrm{s}$ and after $t=0.5\mathrm{s}$ the inlet velocity is zero. The fluid temperature at the inflow corresponds to the ambient temperature in the experiments, $T_{\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}}={20}^{\circ }$ C. These parameters are summarized in Table 2.

4.3.2. Single-phase flushing flow

Next, the mesh resolution for the flushing flow analysis is defined and the single-phase flow is investigated. A grid convergence study based on four grids is performed and the parameters of the grids are given in Table 3. A radial slice of the grid $g_{10}$ is shown in Figure 7(a). The parameters in Table 3 are the length of the smallest isotropic Cartesian cell $h_{\mathrm{c},\mathrm{m}\mathrm{i}\mathrm{n}}$, the total number of cells $n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{s}}$, and the number of time steps to perform the simulation satisfying the stability constraints on the time step. The unstructured grids $\{g_8,g_9,g_{10},g_{11}\}$ have the same ‘pattern’ as that shown in Figure 7(a), i.e. every mesh is based on 4 levels. The subscript m in the mesh notation $g_m$ denotes the number of isotropic refinements performed from an initial square Cartesian cell, which is large enough to encompass the whole computational domain, to the cell with the smallest cell length. That is, the smallest cell in $g_9$ has twice the length as the smallest cell in $g_{10}$.

Table 3. The parameters of the spatial discretizations for the grid convergence study are the smallest cell length $h_{\mathrm{c},\mathrm{m}\mathrm{i}\mathrm{n}}$, the total number of cells $n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{s}}$, and the total number of time steps $\mathrm{\Delta }{t}_n$.

Notation	$h_{\mathrm{c},\mathrm{m}\mathrm{i}\mathrm{n}}$ (µm)	$n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{s}}$	$\mathrm{\Delta }{t}_n$
$g_8$	80	$\approx 0.16\cdot {10}^6$	57,500
$g_9$	40	$\approx 1.2\cdot {10}^6$	115,000
$g_{10}$	20	$\approx 10\cdot {10}^6$	230,000
$g_{11}$	10	$\approx 80\cdot {10}^6$	460,000

Table 4. Geometric, material, and initial condition parameters of the particles.

Material parameters	Value	Unit
Solid density: ${\rho }_{\mathrm{p}}$	7.83	$\left[{\mathrm{g}\mathrm{c}\mathrm{m}}^{-3}\right]$
Solid thermal conductivity: ${\lambda }_{\mathrm{p}}$	50.2	$\left[{\mathrm{W}\mathrm{m}}^{-1}{\mathrm{K}}^{-1}\right]$
Solid specific heat capacity: $c_{\mathrm{p},\mathrm{p}}$	500	$\left[{\mathrm{J}\mathrm{k}\mathrm{g}}^{-1}{\mathrm{K}}^{-1}\right]$
Initial solid temperature: $T_{p,0}$	42.0	$\left[{}^{\circ }\mathrm{C}\right]$
Particle distribution	Value	Unit
Particle diameter: $d_{\mathrm{p}}$	100	[µm]

The criterion for the grid resolution study is the volume integrated turbulent kinetic energy (12) $E_k={\int }_{\mathrm{\Omega }}\frac{1}{2}\rho 〈u_i^{\mathrm{\prime }2}〉\mathrm{d}\mathrm{\Omega }$(12) which is shown as a function of time for each grid in Figure 7(b). The turbulent kinetic energy increases until $t\approx 0.19\mathrm{s}$ and remains constant until the jet is turned off at $t=0.5\mathrm{s}$. Then, $E_k$ decays rapidly until $t\approx 0.65\mathrm{s}$, where it achieves almost zero. Subsequently, the rate of dissipation is reduced and the kinetic energy continues to decay very slowly. A small amount of kinetic energy is still present in the flow at $t=0.8\mathrm{s}$.

The global temporal dynamics of the production, saturation, and dissipation of the kinetic energy is satisfactorily predicted by the meshes $g_9$, $g_{10}$, and $g_{11}$. However, only the two finest grids $g_{10}$ and $g_{11}$ hardly show any difference in the growth rate, the decay rate, and the plateau of the kinetic energy of the velocity fluctuations. Note that although the $g_{10}$ resolution suffices, the grid $g_{11}$ is used throughout this study.

The velocity magnitude and vorticity fields for 8 distinct time instants in the time interval $[0.0\mathrm{s},0.91\mathrm{s}]$ are shown in Figures 8 and 9. During the time interval $[0.0\mathrm{s},0.5\mathrm{s}]$, the fluid is injected into the cavity. The inlet jet impinges upon the bottom of the cavity and generates a tip vortex and a stagnation flow. At $t=0.05\mathrm{s}$, the boundary layer extends over approximately $r=D/4$ where also the weak tip vortex occurs. At the same time, large vortices form next to the core of the jet. While the fluid accelerates, vortex stretching occurs extending the structures to the lateral cavity gap at $t=0.455\mathrm{s}$. During vortex stretching, the boundary layer covering the cavity bottom interacts with the cylindrical cavity walls at $t=0.136\mathrm{s}$ and a torus-like zone of recirculation occurs. The blockage caused by this recirculation decreases the cross-section to transport fluid and particles towards the lateral gaps. This narrow section causes the flow to be accelerated in the time interval $[0.227\mathrm{s},0.5\mathrm{s}]$. The stream is surrounded by the vortex next to the jet and the zone of recirculation.

Figure 8. Velocity magnitude of the flushing flow at the time instants $t=[0.05\mathrm{s},0.136\mathrm{s},0.227\mathrm{s},0.318\mathrm{s},0.455\mathrm{s},0.545\mathrm{s},0.682\mathrm{s},0.910\mathrm{s}]$.

Figure 9. Nondimensional vorticity of the flushing flow at the time instants $t=[0.05\mathrm{s},0.136\mathrm{s},0.227\mathrm{s},0.318\mathrm{s},0.455\mathrm{s},0.545\mathrm{s},0.682\mathrm{s},0.910\mathrm{s}]$. Blue and red regions are counter–rotating.

The structure of the jet and the recirculation region at $t=0.455\mathrm{s}$ can be observed in the streamlines of the projection of the velocity field in a two-dimensional cross section at $t=0.455\mathrm{s}$ in Figure 10. It evidences that the recirculation region on the cavity bottom corner is formed by two small and elongated vortex cores at $r=d_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}/2$ and $r\approx 0.33d_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}/2$ and a larger core at $r\approx 0.66d_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}/2$ which grows till it reaches a height of $s_{\mathrm{h}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}/2$ normal to the wall. The recirculation region on the electrode base is adjacent to the jet and also consists of three vortices. A large core at $r\approx 0.25d_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}/2$ and two small cores close to the jet at $r\approx 0.15d_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}/2$ and $r\approx 0.1d_{\mathrm{c}\mathrm{a}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}/2$ occur.

Figure 10. Streamlines of the two-dimensional axisymmetric projection of the flushing flow at the time instant $t=0.455\mathrm{s}$.

The heat flux through the workpiece surface controls the temperature in the rim zone. This distribution is critical for the resulting surface integrity and the machining precision of the electrical discharge machining process. The average heat flux per surface area between the solid and the fluid is determined by (13) $\overline{q}\left(t\right)=\frac{1}{A}{\int }_A\mathbf{q}(\mathbf{x},t)\cdot \mathbf{n}\mathrm{d}A$(13) where $\mathbf{n}$ is positive in the wall-normal direction. That is, a positive heat flux denotes that heat is transferred from the solid into the fluid. The cooling effect of the single-phase flow on the workpiece is shown in Figure 14, where the average heat flux over a full flushing cycle is shown. The heat flux increases monotonically as the thin boundary layer shown in Figures 8 and 9 develops. When the recirculation region in the cavity corner begins to form at $t\approx 0.2\mathrm{s}$, the area of the cavity bottom that interacts with the cold fluid stream decreases. The heat flux reduces slowly when the recirculation cores begin to form at $t=0.3\mathrm{s}$. This process is accelerated by the growth of the recirculation region. When the jet is turned off at $t=0.5\mathrm{s}$, the heat flux decays to approximately zero due to the dissipation of the kinetic energy.

4.3.3. Multi-phase flushing flow

In the following, the initial condition of the single-phase flushing flow is modified by introducing solid spherical heat conducting particles in the fluid. Three different volume loadings 6%, 9%, and 14% are studied. The particle diameter observed by Brito Gadeschi, Schneider, et al. (2017) is $d_{\mathrm{p}}=100\text{μ}\mathrm{m}$ which is of the same order of magnitude as the particle diameters observed by Tanjilul et al. (2018). Since the total volume in the working gap is $106{\mathrm{m}\mathrm{m}}^3$, this results in 1500, 2300, and 3500 particles in the computational domain.

To quantify the influence of the material debris on the heat transfer through the workpiece during flushing, three homogeneous particle distributions corresponding to a volume loading of 6%, 9%, and 14% are studied. These volume loadings correspond to high material-removal rates ${\dot{V}}_{\mathrm{m}\mathrm{r}}$ in the range of $40{\mathrm{m}\mathrm{m}}^3/\mathrm{m}\mathrm{i}\mathrm{n}\le \dot{\phantom{a}}{V}_{\mathrm{m}\mathrm{r}}\le 140{\mathrm{m}\mathrm{m}}^3/\mathrm{m}\mathrm{i}\mathrm{n}$ (Murray et al., 2016) which corresponds to a range from 1000 to 3500 particles. The particles are radially initialized at a distance of $d_{\mathrm{p}}$ over the cavity bottom surface. The heat transferred from the solid into the fluid during the solidification of the particles is taken into account by using an initial particle temperature of ${42}^{\circ }$C which corresponds to the maximum temperature obtained from the thermographic camera measurement shown in Figure 1(b). Although this initial particle temperature introduces a heat source in the initial condition that is not present in the single-phase flow case, the overall impact on the heat transfer through the workpiece is negligible.

The material parameters of the particles are listed in Table 4. The solid density is ${\rho }_{\mathrm{p}}=7.83\mathrm{g}{\mathrm{c}\mathrm{m}}^{-3}$ which corresponds to a solid-to-fluid density ratio ${\rho }_{\mathrm{p}}/{\rho }_{\mathrm{f}}=10.17$. The thermal conductivity and heat capacity of the particles ${\lambda }_{\mathrm{p}}=50.2$ W/(m K) and $c_{\mathrm{p},\mathrm{p}}=500$ J/(kg K) are required to compute the heat exchange between the particles and the fluid.

The flow around the debris particles is resolved by adaptively refined mesh $g_{11}$ shown in Figure 11. The signed-distance field that captures moving boundaries is used to refine the grid around moving boundaries (Schneiders et al., 2016). Shear layers and steep velocity gradients are detected by the magnitude of the vorticity vector which drives the physics-based adaptation algorithm (Hartmann et al., 2008).

Figure 11. The adaptive mesh around a debris particle in regions of high vorticity. Blue and red regions are counter-rotating.

Next, the particle motion and the effect of the particles on the flushing flow are analyzed. Before the jet impinges on the cavity bottom, the temperature of the particles approximately takes on the local temperature of the fluid. It is shown in Figure 12 that the particles are lifted-off the bottom of the cavity by the initial vortex tip and are transported out of the cavity through the fluid stream into the lateral gap. Figure 13 illustrates the particle distribution for $6\mathrm{\%}$ and $14\mathrm{\%}$ volume loading at the time instant $t=0.455\mathrm{s}$. The particles are colored by the temperature and their positions are projected onto a two-dimensional plane. In the cavity region near the axis, i.e. for the radius $r\le 1.0\mathrm{m}\mathrm{m}$, most particles have already been flushed away at $t=0.455\mathrm{s}$. The remaining particles in this region are located close to the electrode surface at the top of the cavity. These particles were lifted into this region by the tip vortex that is formed in the beginning of the flushing process and are trapped in the vortex ring between the surface of the electrode and the main flow. The pressurized jet continuously cools this region with fresh fluid and the particles attain a temperature of 20–22${}^{\circ }\mathrm{C}$. The temperature of the fluid above the main flow increases monotonously towards the lateral gap up to 25–26${}^{\circ }\mathrm{C}$. For both volume loadings, the particles located close to the workpiece surface possess a higher temperature of 35–36${}^{\circ }\mathrm{C}$ in the area up to the radius $r\le 3.5\mathrm{m}\mathrm{m}$, where the prescribed workpiece temperature reaches ${42}^{\circ }\mathrm{C}$. The flow field underneath the main flow contains two recirculation regions centered at (I) $r=3.25\mathrm{m}\mathrm{m}$ and (II) $r=4.75\mathrm{m}\mathrm{m}$. For both volume loadings, the particle temperature and particle distribution show these recirculation regions. However, for the larger volume loading 14%, the recirculation regions are more pronounced since there are more particles and their average temperature is higher. On average, the particle temperature of the second recirculation region is lower than that of the first for both volume loadings. Finally, the temperature of the particles in the main flow rises from ${25}^{\circ }$C at the center of the working gap to 28–29${}^{\circ }\mathrm{C}$ in the lateral gap. In the lateral gap, the particle temperature decreases since the particles are flushed away due to the higher fluid velocity in the gap.

Figure 12. The fluid temperature (color map) and the particles (grey) located within a cylindrical region of radius $r=2.5\mathrm{m}\mathrm{m}$ of the domain at the time instant $t=0.305\mathrm{s}$.

Figure 13. Planar projection of the particle distribution colored by the particle temperature for a volume loading of (a) 6% and (b) 14% at the time instant $t=0.455\mathrm{s}$. (a) Volume loading 6% and (b) Volume loading 14%.

Figure 14. The distribution of the total heat flux on the surface of the workpiece as a function of time for the single phase flow and the three multi-phase flow cases.

The effect of the particle loading on the spatially averaged heat flux $\overline{\mathbf{q}}$ (13(13) $\overline{q}\left(t\right)=\frac{1}{A}{\int }_A\mathbf{q}(\mathbf{x},t)\cdot \mathbf{n}\mathrm{d}A$(13) ) is shown in Figure 14. During the flushing of the cavity, positive values of the heat flux occur, i.e. the workpiece is cooled by the flow. The average heat flux through the cavity bottom of the three multi-phase flushing flow cases show a similar trend as that of the single-phase flow. Nevertheless, the distribution is significantly more perturbed. Initially, the heat flux increases up to the time $t=0.19\mathrm{s}$. During this phase, the average heat flux is lower than that of the single phase flow. This is due to the smaller temperature gradients that result from the heat released by the initially hot particles close to the cavity bottom. Due to their higher heat capacity, each solid particle stores the same thermal energy as an equivalent isothermal spherical fluid volume of radius $r=144\text{μ}\mathrm{m}$ which has a three times larger volume. When this energy is released into the fluid, its temperature increases and reduces the heat flux through the cavity bottom. Depending on the volume loading, the heat flux reaches its maximum and plateaus in the time interval $[0.19\mathrm{s},0.39\mathrm{s}]$. It oscillates strongly and shows a somewhat chaotic behavior. These oscillations result from increased heat transport in the recirculation regions containing entrained solid particles. Due to their higher heat capacity they are able to transport heat more efficiently into the main fluid stream since they are lifted up by the vortex cores. As the particles are flushed away from the cavity, the heat flux of all cases starts to converge. Finally, when the jet is turned off at $t=0.5\mathrm{s}$, the differences in the heat flux between the multi-phase and the single-phase flow are caused by the inertia of the particles in the recirculation regions which slow down the decay of total kinetic energy in the system.

To quantify the difference between the cases considered, the total average heat flux during the flushing cycle, i.e. the time integral of the surface averaged heat flux, is given in Table 5 for the single-phase and the three multi-phase flow cases. The findings show that the solid particles have a pronounced impact on the heat flux through the workpiece. As the volume loading increases, the total heat flux increases with respect to the single-phase flow configuration.

Table 5. Total heat flux through a flushing cycle.

Flow-type	Value	Unit
Single phase:	50.48	$\left[\mathrm{J}{\mathrm{m}}^{-2}\right]$
Multi-phase 6%:	55.97	$\left[\mathrm{J}{\mathrm{m}}^{-2}\right]$
Multi-phase 9%:	58.13	$\left[\mathrm{J}{\mathrm{m}}^{-2}\right]$
Multi-phase 14%:	58.82	$\left[\mathrm{J}{\mathrm{m}}^{-2}\right]$

The particle removal statistics for a single flushing cycle of $\mathrm{\Delta }t=0.5\mathrm{s}$ are provided in Table 6. Here, $p_{\mathrm{r}\mathrm{e}\mathrm{m}}$ is the total amount of particles flushed, $V_{\mathrm{p},\mathrm{r}\mathrm{e}\mathrm{m}}$ is their volume, and $V_{\mathrm{p},\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}$ is the volume of the initial amount of particles in the cavity. The findings show that the initial volume loading has a pronounced impact on the amount of debris flushed out of the cavity. As the initial volume loading increases, the amount of debris particles flushed $p_{\mathrm{r}\mathrm{e}\mathrm{m}}$ increases. However, this increase grows slower than the initial volume loading as shown by the ratio of the removed and initial volume loadings $V_{\mathrm{p},\mathrm{r}\mathrm{e}\mathrm{m}}/V_{\mathrm{p},\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}$, which decreases as the initial volume loading increases. That is, an increase in the initial volume loading from $6\mathrm{\%}$ to $14\mathrm{\%}$ reduces the relative amount of particles that are flushed from the cavity and therefore the efficiency of the flushing mechanism by $14.2\mathrm{\%}$. This phenomena is caused by the constant surface area of the lateral cavity gap for which the particles need to compete for to escape from the cavity. Assuming a spacing of one particle diameter between the particles when exiting the gap, for the current configuration, the surface area of the lateral cavity gap is the same as that of 768 debris particles. This phenomena is observed in Figure 13(a) at the location where the lateral gap and the working gap meet, where the forming of a particle cluster characteristic for particle choking is found.

Table 6. Debris removal statistics: $p_{\mathrm{r}\mathrm{e}\mathrm{m}}$ is the total number of debris particles removed during a flushing cycle, $V_{\mathrm{p},\mathrm{r}\mathrm{e}\mathrm{m}}$ is their volume, and $V_{\mathrm{p},\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}$ is the volume of the initial amount of particles in the cavity.

Flow-type	$p_{\mathrm{r}\mathrm{e}\mathrm{m}}$	$V_{\mathrm{p},\mathrm{r}\mathrm{e}\mathrm{m}}$	$V_{\mathrm{p},\mathrm{r}\mathrm{e}\mathrm{m}}/V_{\mathrm{p},\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}$
Multi-phase 6%:	986	$516.2\text{μ}{\mathrm{m}}^3$	$65.7\mathrm{\%}$
Multi-phase 9%:	1444	$756.1\text{μ}{\mathrm{m}}^3$	$62.8\mathrm{\%}$
Multi-phase 14%:	1801	$943.0\text{μ}{\mathrm{m}}^3$	$51.5\mathrm{\%}$

5. Conclusions

The efficient removal of material debris by flushing is critical for the performance and resulting surface quality of the electrical discharge machining (EDM) manufacturing process. In this study, a numerical model for quantifying the heat transfer between the workpiece and the flushing flow in a generic die-sink EDM cavity has been introduced. A hierarchical Cartesian cut-cell method that accounts for the fluid-particle heat transfer has been presented and applied to study the single-phase and multi-phase flushing flow in the generic EDM cavity by direct particle–fluid simulations. Three particle volume loadings have been considered. The analysis shows that the dynamics of the recirculation generates a high-speed stream which transports most of the fluid out of the cavity, while simultaneously creating recirculation regions in which particles are entrained. The comparison of the heat transfer from the workpiece into the fluid for the single-phase and the three multi-phase flow cases shows that the debris particles have a pronounced impact on the heat flux through the workpiece. To be more precise, a particle volume loading of 6%, 9%, and 14% increases the mean heat flux by 10%, 15%, and 16%. Furthermore, the analysis shows that the initial volume loading has a pronounced impact on the relative amount of particles that are flushed out of the working gap during a flushing cycle. Since the lateral gap surface remains constant, the particles choke trying to exit the cavity, reducing the overall efficiency of the flushing mechanism. Increasing the initial volume loading from 6% to 14% decreases the relative amount of particles that get flushed out of the cavity by 14%.

This study presents a fully resolved two-phase flow model to analyze the flushing flow in electrical discharge machining. This work is currently extended along two research directions, to experimentally validate the two-phase flow model presented in this work and to further develop the present method such that three-phase solid–fluid-gas flushing flows can be analyzed by fully-resolved simulations.

Acknowledgements

This work has been funded by the German Research Foundation (DFG) within the framework of the SFB/Transregio 136 ‘Process Signatures’ (subproject M04). This support is gratefully acknowledged. Computing resources were provided by the High Performance Computing Center Stuttgart (HLRS) and by the Jülich Supercomputing Center (JSC) within a Large-Scale Project of the Gauss Center for Supercomputing (GCS).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by Deutsche Forschungsgemeinschaft [grant number 223500200].