In Section 3, the promotion influence of a flow field induced by highpressure inner flushing on a fast EDmilling process is proven by the analysis of crater morphology. However, the dynamic evolution process of molten material in the combined crater under the effect of a fluid field remains unclear. Without this it is impossible to explain the directionality of the combined craters and the variation trend of Vr with respect to Pin. Therefore, on the platform of COMSOL Multiphysics 6.0, a novel thermalfluid coupling model, which incorporates a highpressure inner flushing and the special geometric properties, is developed. This thermalfluid model is used to investigate a complete process of generation, movement, and solidification of molten material that generated by a single discharge in a multipulse discharge experiment. A conjugate heat transfer module, a turbulent twophase flow module, a phase change module, and a levelset method are combined to simulate an evolution process of molten material. Based on the temperature, velocity, and pressure fields from the model, the moving tendency of molten material and the phenomenon occurring during the abovementioned experiments can be explained. The details of the model are illustrated as follows.
4.1 Assumptions
To establish the coupling model, some reasonable assumptions and declarations are needed to be made first to simplify the calculation.
 The workpiece material is homogeneous and isotropic, and the liquid metal is a Newtonian laminar liquid that cannot be compressed.
 The vaporization of the metal is not considered, when the temperature of the solid metal exceeds the melting point, it is completely transformed into molten material after the latent heat of fusion is absorbed.
 Similar to Ref. [11], in this simulation, the radius of a discharge plasma channel remains constant during the whole pulse duration time.
 The discharge maintaining voltage and discharge current remain constant during the entire pulse duration time.
 The influence of side flushing is not considered.
4.2 Geometric modeling
According to the crosssection profiles of craters measured in Section 3, a simplified geometric model which schematically represents the multipulse discharge experiments carried out in fast EDmilling can be developed, as shown in Fig. 17 (a). For the sake of simplifying the calculation, a 2D axially symmetrical structure as shown in Fig. 17 (b) is adopted in the thermalfluid coupling model. Dominant parameters, such as gap width, machining depth, etc., are measured in actual machining, as listed in Table 2.
Table 2
The geometric dimension of the axisymmetric model.
Parameter

Electrode outer diameter

Electrode inner diameter

Gap width

Machining depth

Uneroded height

Value (mm)

0.6

0.2

0.05

0.025

0.025

4.3 Heat flux
In terms of the mechanism of material removal, fast EDmilling is in line with other EDM processes. During a discharge process, heat flux is generated by a discharge plasma channel and transferred from the surface of a workpiece to the gap and the body of the workpiece [19]. As noted by Descoeudres et al. [20] and Kojima [21], the heat flux in a discharge plasma channel obeys a Gaussian distribution. Moreover, it has been proven that as compared to a volume heat source, a surface heat source can better explain the phenomenon occurring during in a discharge, and the diameter of a heat source is the same as that of a plasma channel [22]. Therefore, based on the assumptions made in Section 4.1, a Gaussian surface heat source adopted in the coupling model is described as
$$Q(r)=\frac{{4.57P}}{{\pi {R^2}}}\exp (  \frac{{4.5{{(r  {r_0})}^2}}}{{{R^2}}})$$
3
where r represents the coordinate of a point under the cylindricalcoordinate system (m), \({r_0}\)denotes the coordinate of a plasma channel center, R is the radius of the plasma channel (m), under the conditions used in this paper. R is found to be about 50µm. P denotes the discharge power (W), which can be calculated by the product of the discharge maintaining voltage (18V) and the discharge current (10A).
To simulate the pulse duration and pulse interval, a piecewise function is utilized to control the action time of the heat flux, as described in
$$h(t)=\left\{ \begin{gathered} 1{\text{ }}0<t \leqslant {T_{on}} \hfill \\ 0{\text{ }}{T_{on}}<t \leqslant {T_{on}}+{T_{off}} \hfill \\ \end{gathered} \right.$$
4
In a pulse discharge process, the generated discharge energy is distributed among the anode, cathode, and gap. According to Xia et al. [23], when a negative tool polarity is adopted, the energy is distributed with an approximate ratio of 1:1.4:1 among the tool electrode, workpiece, and gap. Therefore, the heat flux exposed on the workpiece can be expressed as
$${Q_w}(r,t)=\frac{{846.794}}{{\pi {R^2}}}\exp (  \frac{{4.5{{(r  {r_0})}^2}}}{{{R^2}}})h(t)$$
5
A boundary heat source in the heat transfer in solids and fluids module is applied to implement the heat flux. During a discharge pulse duration, the solid workpiece is continuously heated. When the temperature exceeds the melting point, the heated part is completely transformed from solid superalloy to molten metal. To determine the melting boundary of the workpiece, a phase change module is also applied. Thermophysical properties of nickelbased superalloys are listed in Table 3.
Table 3
Thermophysical properties of nickelbased superalloys.
Property (unit)

Value

Density (kg/m3)

8180

Melting point (K)

1798

Thermal conductivity (W/m/K)

11.5

Heat capacity at constant pressure (J/kg/K)

446

Latent heat of fusion (J/kg)

3.33×105

Dynamic viscosity of liquid phase (Pa·s)

0.001 [24]

Ratio of specific heats

1

4.4 Molten material tracking
To simulate a movement process of the generated molten material under the effect of a flushing fluid, the type of a flow pattern when flowing through the gap channel needs to be determined first as laminar or turbulent. The type of flow pattern can be determined with the value of the Reynolds number as
$$\operatorname{Re} =\frac{{\rho vd}}{\mu }$$
6
where\(\rho\)(kg/m3),(m/s), and\(\mu\)(Pa·s) represent the density, velocity, and viscosity coefficient of the fluid, respectively, and d (m) denotes the gap width. Under a flushing pressure of 0.48Mpa, the average flow velocity is measured to be about 25m/s, and the corresponding Reynolds number is about 4950. Therefore, in this paper, a Reynoldsaveraged Navier–Stokes (RANS) model is adopted to simulate the flow field within a discharge gap and the dynamics of the turbulence are described by a standard twoequation kε model with realizability constraints.
To track the movement of molten material, a levelset method [25] is used. The basic idea of this method is to track moving interfaces in a fluidflow model by solving a transport equation for the level set function\(\phi (p,t)\). In the coupling model developed in this paper,\(\phi (p,t)\)is defined as
$$\phi (p,t)\left\{ \begin{gathered} 1{\text{ }}r \in flushing{\text{ }}liquid \hfill \\ 0{\text{ }}r \in molten{\text{ }}material \hfill \\ \end{gathered} \right.$$
7
where p denotes the position coordinate (m), and t represents the time (s). To update\(\phi (p,t)\)by calculating the velocity field of molten material, the transport equation can be written as
$$\frac{{\partial \phi }}{{\partial t}}+\mathbf{u} \cdot \nabla \phi =\gamma \nabla \cdot \left( {{\varepsilon _{ls}}\nabla \phi  \phi (1  \phi )\frac{{\nabla \phi }}{{\left {\nabla \phi } \right}}} \right)$$
8
where\(\phi\)is the levelset function introduced above, u denotes the moving speed of the molten material (m/s),\(\gamma\)is the reinitialization speed (m/s), and\({\varepsilon _{ls}}\)is the parameter controlling interface thickness.
4.5 Governing equations
In the developed thermalfluid coupling model, heat transfer in solids and fluids, turbulent twophase flow with a level set method, and phase change are comprehensively considered. To achieve a combined solution of the temperature, flow, and velocity field, the RANS equation for conservation of momentum, the continuity equation for conservation of mass, and the heat transfer equation, as described in Eqs. (9) through (11), need to be simultaneously satisfied.
$$\rho \frac{{\partial \mathbf{u}}}{{\partial t}}+\rho (\mathbf{u} \cdot \nabla )\mathbf{u}=\nabla \cdot \left[ {  p\mathbf{I}+\mu (\nabla \mathbf{u}+{{(\nabla \mathbf{u})}^T}} \right]+\mathbf{F}$$
9
$$\rho \nabla \cdot (\mathbf{u})=0$$
10
$$\rho {C_P}\frac{{\partial T}}{{\partial t}}+\rho {C_P}\mathbf{u} \cdot \nabla T=k\nabla \cdot \nabla T+Q$$
11
where\(\rho\)is density (kg/m3), t is time (s), u denotes the velocity filed (m/s), p is the pressure (Pa), I is the identity matrix,\(\mu\)is dynamic viscosity (Pa·s),\({(\nabla \mathbf{u})^T}\)represents the transposed matrix of\(\nabla \mathbf{u}\), F denotes the volume force,\({C_P}\)heat capacity at constant pressure (J/kg/K), T is the absolute temperature (K), k denotes thermal conductivity (W/m/K), and Q is the heat source (W/m).
4.6 Boundary conditions
Boundary conditions for the thermalfluid coupling model are schematically illustrated in Fig. 18. The computational domain of this model consists of three parts, namely the electrode (brass), gap channel (water), and workpiece (nickelbased superalloy). As for the thermal boundary conditions, boundaries in the model are classified into four categories and each category is shaded with different colors enclosed with solid lines. As shown in Fig. 18, a purple line represents a heat flux across boundary, a black line represents a thermal insulation boundary, a red line represents a heat conduction boundary, and a blue line represents a heat convection boundary. A heat conduction boundary where the boundary heat source is applied can be defined as
$$k\mathbf{n} \cdot \nabla T={Q_w}(r,t)$$
12
A thermal insulation boundary is defined by Eq. (13). A heat convention boundary and a heat flux across boundary are defined by the same Eq. (14) with different heat transfer coefficients.
$$k\mathbf{n} \cdot \nabla T=0$$
13
$$k\mathbf{n} \cdot \nabla T={h_c}(T  {T_0})$$
14
where hc denotes the heat transfer coefficient (W/m2/K), T0 is the ambient temperature (293.15K), n represents the normal direction to the boundary, and other parameters are the same as described in Section 4.5.
With regard to the fluid boundary conditions, as demonstrated in Fig. 18, the inlet and outlet of a flushing fluid can be respectively defined as
$${\mathbf{n}^T}\left[ {  p\mathbf{I}+\mu (\nabla \mathbf{u}+{{(\nabla \mathbf{u})}^T}} \right]\mathbf{n}=  p_{0}^{{in}}$$
15
$$\left[ {  p\mathbf{I}+\mu (\nabla \mathbf{u}+{{(\nabla \mathbf{u})}^T}} \right]\mathbf{n}=  p_{0}^{{out}}\mathbf{n}$$
16
where\({\mathbf{n}^T}\)is the transpose of a boundary normal vector n pointing out of the domain,\(p_{0}^{{in}}\)is the flushing pressure (0.48MPa), and\(p_{0}^{{out}}\)represents the ambient pressure (1atm). The remaining boundaries of the flow field except for the symmetrical axis are defined as walls (depicted with blue lines), which can be formulated by