**1. Introduction**

With the development of aerospace and automotive manufacturing and machinery production miniaturization, precision and complex internal characteristics of the small hole structure are more and more widely used, such as small holes in the reaming form of it can be applied to the outer ring of the miniature bearing [1]. GH4169 is a high-temperatureresistant nickel-based alloy material widely used to manufacture miniature bearing parts [2]. However, due to the high hardness of the alloy in the use of conventional machining methods challenging to the process, poor machining accuracy, and surface quality, the service life of miniature bearings is not long [3]. Electric discharge machining (EDM) and laser beam machining (LBM) are both thermal processes with high machining efficiency. The disadvantages are that they produce recast layers, heat-affected areas, and tensile residual stresses that also reduce the service life of micro bearings [4–7]. Electrochemical machining (ECM) is based on the principle of anodic dissolution for metal removal, independent of the hardness of the workpiece. It has become one of the main machining techniques for

**Citation:** Li, Z.; Li, W.; Cao, B. Simulation Analysis of Multi-Physical Field Coupling and Parameter Optimization of ECM Miniature Bearing Outer Ring Based on the Gas-Liquid Two-Phase Turbulent Flow Model. *Micromachines* **2022**, *13*, 902. https://doi.org/ 10.3390/mi13060902

Academic Editors: Youqiang Xing, Xiuqing Hao and Duanzhi Duan

Received: 5 May 2022 Accepted: 6 June 2022 Published: 7 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

machining miniature bearing outer rings on the high-temperature-resistant nickel-based alloy GH4169 because of its advantage of no thermal damage [8–10]. There are still some problems to be solved. For example, the influence of electrolyte temperature distribution and bubble rate distribution on electrolyte conductivity in the ECM area cannot be fully considered, resulting in the simulation model not being able to accurately predict the machining accuracy of the outer ring of the miniature bearing. Therefore, it is difficult to model and predict the optimal process parameters.

In Deconck et al. [11–13], a simulation model for calculating the temperature distribution was developed, and it was pointed out that the electrolyte temperature distribution has an important influence on the machining accuracy. However, the model ignores the effect of electrolyte bubble rate distribution on machining accuracy and uses the laminar flow N–S model, which reduces the heat transfer effect. In Fang et al. [14], a multiphysics field coupled simulation model was proposed to predict the electrolytic machining accuracy. However, they treated the electrolyte as a single-phase flow and neglected the effect of electrolyte bubble rate distribution on machining accuracy. In Gomez-Gallegos et al. [15], a 3D coupled multiphysics field finite element model was developed to predict ECM accuracy. In Li et al. [16], a coupled model of the magnetic field, electric field, and electrolyte flow in ECM was developed to predict the accuracy of ECM. However, both of these researchers neglected the influence of the bubble rate distribution of the electrolyte in the machining area on the ECM accuracy, resulting in the simulation model's inaccurate prediction of the ECM accuracy. Kozak et al. [17] and Mayank et al. [18] pointed out that the ECM accuracy is also affected by the bubble rate distribution of the electrolyte in the machining area. Some researchers have conducted studies on the bubble rate distribution of electrolytes in the machining area. In Shimasaki et al. [19] and Zhang et al. [20], the effect of bubbles generated by the electrolyte in the machining area on the accuracy of the ECM is observed using transparent electrodes. In Chang et al. [21], a two-dimensional two-phase laminar flow quasi steady-state flow model is proposed to predict the accuracy of the ECM. In Klocke et al. [22], the flow field in the machining area is assumed to be a gas–liquid two-phase flow. The bubble rate distribution in the electrolyte is approximated using the laminar bubble flow model. The use of the laminar N–S model reduces the heat transfer effect, making the prediction accuracy of the simulation model low. In Chen et al. [23] and Zhou et al. [24], a multiphysics field coupled simulation model for the ECM of turbine blades was established based on the gas–liquid two-phase turbulent flow k-ε model, and the effects of electrolyte temperature distribution and bubble rate distribution on the accuracy of the ECM of turbine blades were fully considered. The results show that the simulation model is highly accurate in predicting the machining accuracy of turbine blades.

The complexity of the ECM makes it challenging to model and predict the optimal process parameters [25]. In addition, the random selection of process parameters or trialand-error methods is very costly and time-consuming and does not yield the desired results. These problems can be solved by optimization techniques [26,27]. Jain et al. [28] used a genetic algorithm to optimize the ECM process parameters. Three process parameters, namely, tool cathode feed rate, electrolyte flow rate, and machining voltage were selected as input quantities and machining accuracy as output quantities, and the optimization results obtained showed significant improvement in machining accuracy. In Jegan et al. [29], the optimization of the ECM process parameters based on the particle swarm algorithm was investigated. Four process parameters, namely, machining current, machining voltage, electrolyte concentration, and tool cathode feed rate were selected as input quantities and material removal rate and surface roughness as output quantities, and the particle swarm algorithm was determined to be superior to the genetic algorithm in terms of computation time and statistical analysis. In Mehrvar et al. [30], based on the central composite design to optimize the ECM process parameters, four process parameters, namely, machining voltage, tool cathode feed rate, electrolyte flow rate, and electrolyte concentration were selected as input quantities, and the material removal rate and surface roughness were

determined as output quantities. The results show that the proposed optimization method is effective and suitable. and surface roughness were determined as output quantities. The results show that the proposed optimization method is effective and suitable. In summary, the simulation model of ECM of the microbearing outer ring under the

central composite design to optimize the ECM process parameters, four process parameters, namely, machining voltage, tool cathode feed rate, electrolyte flow rate, and electrolyte concentration were selected as input quantities, and the material removal rate

In summary, the simulation model of ECM of the microbearing outer ring under the influence of electrolyte temperature distribution and bubble rate distribution is established based on the gas–liquid two-phase turbulent flow k-ε model to predict its machining accuracy. It is also essential to optimize the process parameters of the ECM of the miniature bearing outer ring through central composite design to improve its machining accuracy. In this paper, a multiphysics field coupled simulation model of electric, flow, and temperature fields during the ECM of the miniature bearing outer ring is established based on the gas–liquid two-phase turbulent flow model. The influence of electrolyte temperature distribution and bubble rate distribution on the accuracy of the ECM is fully considered. The process of the ECM of miniature bearing outer rings was investigated using a central composite design, and a regression equation between the mean error and the process parameters was established. The optimal combination of process parameters for the mean error was predicted. influence of electrolyte temperature distribution and bubble rate distribution is established based on the gas–liquid two-phase turbulent flow k-ε model to predict its machining accuracy. It is also essential to optimize the process parameters of the ECM of the miniature bearing outer ring through central composite design to improve its machining accuracy. In this paper, a multiphysics field coupled simulation model of electric, flow, and temperature fields during the ECM of the miniature bearing outer ring is established based on the gas–liquid two-phase turbulent flow model. The influence of electrolyte temperature distribution and bubble rate distribution on the accuracy of the ECM is fully considered. The process of the ECM of miniature bearing outer rings was investigated using a central composite design, and a regression equation between the mean error and the process parameters was established. The optimal combination of process parameters for the mean error was predicted.

### **2. Simulation Model of ECM of the Miniature Bearing Outer Ring 2. Simulation Model of ECM of the Miniature Bearing Outer Ring**

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 3 of 16

### *2.1. Geometric Model 2.1. Geometric Model*

Since the entire ECM miniature bearing outer ring is modeled as an axisymmetric figure, it is assumed that the multiphysics field coupling is the same for each cross section. Therefore, the model is simplified in two dimensions for the convenience of analysis and calculation. The geometric model of the machining area is shown in Figure 1. The part inside the red dashed box in the figure is taken for analysis, and the shape of the entire machining area can be obtained by rotating it around the tool cathode axis for one week. Boundary Γ<sup>1</sup> is the tool cathode without an insulating layer; boundaries Γ2, Γ3, Γ4, and Γ<sup>5</sup> are the tool cathode with an insulating layer; boundary Γ<sup>6</sup> is the electrolyte inlet; boundary Γ<sup>7</sup> is the electrolyte outlet; boundary Γ<sup>8</sup> is the workpiece anode, i.e., the outer ring of the miniature bearing, and Ω is the ECM area. Since the entire ECM miniature bearing outer ring is modeled as an axisymmetric figure, it is assumed that the multiphysics field coupling is the same for each cross section. Therefore, the model is simplified in two dimensions for the convenience of analysis and calculation. The geometric model of the machining area is shown in Figure 1. The part inside the red dashed box in the figure is taken for analysis, and the shape of the entire machining area can be obtained by rotating it around the tool cathode axis for one week. Boundary Γ<sup>1</sup> is the tool cathode without an insulating layer; boundaries Γ2, Γ3, Γ4, and Γ<sup>5</sup> are the tool cathode with an insulating layer; boundary Γ<sup>6</sup> is the electrolyte inlet; boundary Γ<sup>7</sup> is the electrolyte outlet; boundary Γ<sup>8</sup> is the workpiece anode, i.e., the outer ring of the miniature bearing, and Ω is the ECM area.

**Figure 1.** Geometric model of the machining area. **Figure 1.** Geometric model of the machining area.

### *2.2. Mathematical Model 2.2. Mathematical Model*

2.2.1. Mathematical Model of Electric Field

2.2.1. Mathematical Model of Electric Field According to Ohm's law, the relationship between the current density in the machining area and the electric field strength and potential is:According to Ohm's law, the relationship between the current density in the machining area and the electric field strength and potential is:

$$E = -\nabla \varphi \tag{1}$$

$$
\dot{\iota} = \sigma \mathbf{E} = -\sigma \nabla \varphi \tag{2}
$$

where *E* is the electric field strength; *ϕ* is the electrolyte potential; *i* is the current density; *σ* is the electrolyte conductivity.

In the actual ECM, the speed of the ECM is usually expressed in terms of the dissolution speed in the direction normal to the metal surface of the workpiece anode:

$$
v\_{\boldsymbol{n}} = \eta \omega \mathbf{i} = -\eta \omega \boldsymbol{\sigma} \nabla \boldsymbol{\varphi} \tag{3}$$

where *v<sup>n</sup>* is the ECM speed; *η* is the ECM efficiency; *ω* is the volumetric galvanic equivalent of the workpiece anode.

### 2.2.2. Mathematical Model of the Flow Field

The gas phase and solid phase products generated during the electrolysis process form a three-phase flow of gas, liquid, and solid in the machining area. In contrast, the volume ratio of the solid phase electrolysis products is minimal and has minimal effect on the electrolyte conductivity, so the flow field in the machining area can be simplified to a gas–liquid two-phase flow.

The gas–liquid two-phase flow in the ECM satisfies the conservation of mass:

$$\frac{\partial}{\partial t} \left( \beta\_l \rho\_l + \beta\_\mathcal{S} \rho\_\mathcal{S} \right) + \nabla \cdot \left( \beta\_l \rho\_l u\_l + \beta\_\mathcal{S} \rho\_\mathcal{S} u\_\mathcal{S} \right) = 0 \tag{4}$$

$$\frac{\partial \beta\_{\mathcal{S}} \rho\_{\mathcal{S}}}{\partial t} + \nabla \cdot (\beta\_{\mathcal{S}} \rho\_{\mathcal{S}} u\_{\mathcal{S}}) = m\_{\mathcal{l}g} \tag{5}$$

$$\frac{\partial \beta\_l \rho\_l}{\partial t} + \nabla \cdot (\beta\_l \rho\_l u\_l) = -m\_{l\text{g}}\tag{6}$$

$$
\beta\_l + \beta\_{\mathcal{S}} = 1 \tag{7}
$$

where *β<sup>g</sup>* is the proportion of gas in the total volume of the two-phase flow; *β<sup>l</sup>* is the proportion of liquid in the total volume of the two-phase flow; *ρ<sup>g</sup>* and *ρ<sup>l</sup>* are gas density and liquid density; *u<sup>g</sup>* and *u<sup>l</sup>* are the velocity of the gas phase and liquid phase; *mlg* is the mass transfer rate of the liquid phase into the gas phase.

The gas–liquid two-phase flow satisfies the conservation of momentum:

$$\frac{\partial}{\partial t} \left( \beta\_{\mathcal{S}} \rho\_{\mathcal{S}} u\_{\mathcal{S}} \right) + \nabla \cdot \left( \beta\_{\mathcal{S}} \rho\_{\mathcal{S}} u\_{\mathcal{S}} u\_{\mathcal{S}} \right) = -\beta\_{\mathcal{S}} \nabla p + \nabla \cdot \tau\_{\mathcal{S}} + \beta\_{\mathcal{S}} \rho\_{\mathcal{S}} \mathbf{g} + F\_{\mathcal{U}} \tag{8}$$

$$\frac{\partial}{\partial t}(\beta\_l \rho\_l u\_l) + \nabla \cdot (\beta\_l \rho\_l u\_l u\_l) = -\beta\_l \nabla p + \nabla \cdot \tau\_l + \beta\_l \rho\_l g - F\_m \tag{9}$$

where *p* is the electrolyte pressure; *τ<sup>g</sup>* and *τ<sup>l</sup>* are gas and liquid viscous stress tensors; *F<sup>m</sup>* is the interphase force.

Assuming that hydrogen is produced only on the surface of the tool cathode in the ECM and that the pressure and temperature distributions of the gas and liquid phases are the same, Faraday's law states that:

$$m\_H = k\_H I t = k\_H iSt \tag{10}$$

where *m<sup>H</sup>* is the mass of hydrogen produced; *k<sup>H</sup>* is the electrochemical equivalent of the hydrogen mass; *I* is the current; *i* is the current density; *S* is the area of the tool cathode.

The mass transfer rate *mlg* for the transformation of liquid phase into gas phase is the mass flux of hydrogen gas produced on the cathode per unit width of the tool, which combined with Equation (10) can be obtained after finishing:

$$Lm\_{l\mathcal{K}} = ik\_H \tag{11}$$

where *L* is the width of the tool cathode unit.

The density of hydrogen is calculated from the ideal gas equation of state as:

$$
\rho\_{\mathcal{S}} = \frac{p}{RT} \tag{12}
$$

where *R* is the gas constant; *T* is the electrolyte temperature.

The electrolyte in the ECM region is in a turbulent state, and considering the effect of bubbles, the RANS *k*-*ε* turbulence model is used in this paper. The electrode surface near the wall is solved by the wall function as follows:

$$\frac{\partial k}{\partial t} + \nabla \cdot \left[ku - \left(\mu + \frac{\mu\_T}{\sigma\_k}\right) \nabla k\right] = P\_k + S\_k - \varepsilon \tag{13}$$

$$\frac{\partial \varepsilon}{\partial t} + \nabla \cdot \left[ \varepsilon u - \left( \mu + \frac{\mu\_T}{\sigma\_\varepsilon} \right) \nabla \varepsilon \right] = \frac{\varepsilon}{k} (\mathsf{C}\_1 P\_k + \mathsf{C}\_\varepsilon \mathsf{S}\_k - \mathsf{C}\_2 \varepsilon) \tag{14}$$

$$P\_k = \frac{\mu\_T}{2} \left| \nabla u + \left( \nabla u \right)^T \right|^2 \tag{15}$$

$$\mathcal{S}\_k = -\beta \mathcal{C}\_k |\nabla p|^2 \tag{16}$$

where *k* is the turbulent kinetic energy; *ε* is the turbulent dissipation rate; *u* is the electrolyte flow rate; *µ* is the electrolyte dynamic viscosity; *µ<sup>T</sup>* is the turbulent viscosity coefficient; *C*1, *C*2, *C<sup>k</sup>* , *C<sup>ε</sup>* , *σ<sup>k</sup>* , and *σ<sup>ε</sup>* are the model constants.

### 2.2.3. Mathematical Model of the Temperature Field

During ECM, fluid heat transfer occurs in the machining area. The heat generated in the machining area is carried away by the thermal convection of the electrolyte in the turbulent state. According to the law of energy conservation, the expression of the convective heat transfer equation in the machining area is:

$$
\rho \mathbb{C}\_p \frac{\partial T}{\partial t} + \rho \mathbb{C}\_p u \cdot \nabla T = \nabla(\lambda \nabla T) + Q \tag{17}
$$

$$Q = i\nabla\varphi\tag{18}$$

where *ρ* is the density of the electrolyte; *C<sup>p</sup>* is the specific heat capacity of the electrolyte; *λ* is the thermal conductivity of the electrolyte; *Q* is the heat generated in the ECM.

### 2.2.4. Multiphysics Field Coupling Model

The conductivity of the electrolyte in the actual ECM as affected by temperature and bubble rate can be expressed as:

$$
\sigma = \sigma\_0 (1 - \beta)^m [1 + \gamma (T - T\_0)] \tag{19}
$$

where *σ<sup>0</sup>* is the initial conductivity of the electrolyte; *β* is the bubble rate; *m* is the bubble rate influence index; *γ* is the temperature correlation coefficient; *T*<sup>0</sup> is the initial temperature of the electrolyte.

Substituting Equation (19) into Equation (3), the coupling equations for the electric field, flow field, temperature field, and ECM speed are obtained:

$$
\sigma\_{\rm tr} = -\eta \omega \sigma\_0 (1 - \beta)^m [1 + \gamma (T - T\_0)] \nabla \varphi \tag{20}
$$

## *2.3. 2D COMSOL Multiphysics Field Coupling Simulation Model*

As shown in Figure 2, the electric field, flow field, temperature field, and deformation geometry modules on the software were selected for multiphysics field coupling simulation. The coupling method is that the electric and temperature field modules were coupled through an electromagnetic heat source module, and the flow and temperature field modules were coupled through a non-isothermal flow module. The actual ECM of the outer ring of the miniature bearing can be reproduced to the maximum extent.

extent.

**Figure 2.** Simulation model and simulation module. **Figure 2.** Simulation model and simulation module.

*2.3. 2D COMSOL Multiphysics Field Coupling Simulation Model*

As shown in Figure 2, the electric field, flow field, temperature field, and deformation geometry modules on the software were selected for multiphysics field coupling simulation. The coupling method is that the electric and temperature field modules were coupled through an electromagnetic heat source module, and the flow and temperature field modules were coupled through a non-isothermal flow module. The actual ECM of the outer ring of the miniature bearing can be reproduced to the maximum

### 2.3.1. Boundary Condition Setting 2.3.1. Boundary Condition Setting

The boundary conditions of the electric field module were set as follows: the boundary Γ<sup>8</sup> was connected to the voltage U; the boundary Γ<sup>1</sup> was grounded, and the boundaries Γ2, Γ3, Γ4, Γ5, Γ6, and Γ<sup>7</sup> were electrically insulated. The boundary conditions of the flow field module were set as follows: boundary Γ<sup>1</sup> gas mass flux, electrolyte wall function; boundary Γ<sup>6</sup> electrolyte normal inflow velocity *u0*, no gas flux; boundary Γ<sup>7</sup> electrolyte and gas outlet; and boundaries Γ2, Γ3, Γ4, Γ5, and Γ<sup>8</sup> no gas flux, electrolyte wall function. The temperature module boundary conditions were set as follows: boundaries Γ1, Γ2, Γ3, Γ4, Γ5, and Γ<sup>8</sup> thermal insulation; boundary Γ<sup>6</sup> electrolyte initial temperature *T0*; boundary Γ<sup>7</sup> electrolyte outflow. The boundary conditions of the deformation geometry module were set as follows: boundary Γ<sup>8</sup> normal mesh moving speed, i.e., electrochemical processing speed *vn*. The boundary conditions of the electric field module were set as follows: the boundary Γ<sup>8</sup> was connected to the voltage U; the boundary Γ<sup>1</sup> was grounded, and the boundaries Γ2, Γ3, Γ4, Γ5, Γ6, and Γ<sup>7</sup> were electrically insulated. The boundary conditions of the flow field module were set as follows: boundary Γ<sup>1</sup> gas mass flux, electrolyte wall function; boundary Γ<sup>6</sup> electrolyte normal inflow velocity *u*0, no gas flux; boundary Γ<sup>7</sup> electrolyte and gas outlet; and boundaries Γ2, Γ3, Γ4, Γ5, and Γ<sup>8</sup> no gas flux, electrolyte wall function. The temperature module boundary conditions were set as follows: boundaries Γ1, Γ2, Γ3, Γ4, Γ5, and Γ<sup>8</sup> thermal insulation; boundary Γ<sup>6</sup> electrolyte initial temperature *T*0; boundary Γ<sup>7</sup> electrolyte outflow. The boundary conditions of the deformation geometry module were set as follows: boundary Γ<sup>8</sup> normal mesh moving speed, i.e., electrochemical processing speed *vn*.

### 2.3.2. Material Parameter Setting of the Simulation Model

2.3.2. Material Parameter Setting of the Simulation Model The workpiece anode material of the simulation model was GH4169 high-temperature-resistant nickel-based alloy; the tool cathode material was titanium alloy electrode (insulating film covered with PTFE), and the electrolyte was NaNO<sup>3</sup> The workpiece anode material of the simulation model was GH4169 high-temperatureresistant nickel-based alloy; the tool cathode material was titanium alloy electrode (insulating film covered with PTFE), and the electrolyte was NaNO<sup>3</sup> solution of a given concentration. The material parameters of the specific simulation model are shown in Table 1.

solution of a given concentration. The material parameters of the specific simulation


model are shown in Table 1. **Table 1.** Material parameter settings of the simulation model.

### **3. Simulation Analysis of ECM of the Miniature Bearing Outer Ring**

The simulation conditions were as follows: machining voltage was 18 V; the electrolyte was NaNO<sup>3</sup> solution with 12% concentration, and the electrolyte inlet flow rate was 9 m/s.

### *3.1. Electrolyte Temperature Distribution in the Machining Area* As shown in Figure 3, the maximum electrolyte temperature in the machining area is

*3.1. Electrolyte Temperature Distribution in the Machining Area*

**3. Simulation Analysis of ECM of the Miniature Bearing Outer Ring**

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 7 of 16

Gas density (kg/m<sup>3</sup>

GH4169 volumetric electrochemical equivalent

(cm<sup>3</sup>

rate was 9 m/s.

As shown in Figure 3, the maximum electrolyte temperature in the machining area is 293.6 K, which is 0.45 K higher than the initial electrolyte temperature of 293.15 K. The maximum temperature of the electrolyte in the machining area from the initial stage of machining to the machining time of 20 s increased by 0.1 K, and the area of the thermally affected area was about 45% of the machining area. This is because more Joule heat is generated at this stage by higher current density, and the machining gap is small, resulting in poor electrolyte circulation and the accumulation of Joule heat in the direction of the electrolyte flow, causing the electrolyte temperature to increase and the area of the heataffected area to expand. From the machining time of 20 s to the machining time of 60 s, the maximum temperature of the electrolyte in the machining area is reduced by 0.2 K, and the area of the heat-affected area is about 35% of the machining area. This is because the current density at this stage is low; the Joule heat generated is low, and the machining gap is large. The electrolyte flow is smooth, so that the electrolysis temperature is lower, and the area affected by heat is reduced. 293.6 K, which is 0.45 K higher than the initial electrolyte temperature of 293.15 K. The maximum temperature of the electrolyte in the machining area from the initial stage of machining to the machining time of 20 s increased by 0.1 K, and the area of the thermally affected area was about 45% of the machining area. This is because more Joule heat is generated at this stage by higher current density, and the machining gap is small, resulting in poor electrolyte circulation and the accumulation of Joule heat in the direction of the electrolyte flow, causing the electrolyte temperature to increase and the area of the heat-affected area to expand. From the machining time of 20 s to the machining time of 60 s, the maximum temperature of the electrolyte in the machining area is reduced by 0.2 K, and the area of the heat-affected area is about 35% of the machining area. This is because the current density at this stage is low; the Joule heat generated is low, and the machining gap is large. The electrolyte flow is smooth, so that the electrolysis temperature is lower, and the area affected by heat is reduced.

) 8.99 × 10−<sup>2</sup>

Air bubble diameter (m) 1 × 10−<sup>5</sup> Bubble rate impact index 1.5

/A/min) 0.00178

The simulation conditions were as follows: machining voltage was 18 V; the electrolyte was NaNO<sup>3</sup> solution with 12% concentration, and the electrolyte inlet flow

**Figure 3.** Electrolyte temperature distribution in the machining area. **Figure 3.** Electrolyte temperature distribution in the machining area.

### *3.2. Distribution of Electrolyte Bubble Rate in the Machining Area*

*3.2. Distribution of Electrolyte Bubble Rate in the Machining Area* As shown in Figure 4, the bubble rate of electrolytes in the machining area gradually increases as the ECM proceeds. The highest bubble rate of electrolytes in the machining area reaches 18.5% at 60 s of machining time. The highest bubble rate of electrolytes in the machining area from the initial machining stage to 20 s of machining time reached 17.3%, and the electrolyte bubbles were distributed around the boundary Γ<sup>1</sup> with a minimal area. This is because the high flow rate of the electrolyte at this stage can carry away the hydrogen produced at the boundary Γ<sup>1</sup> in time. The highest bubble rate of electrolytes in the machining area increased by 1.2% from 20 s of machining time to 60 s of machining time, and the area of electrolyte bubble distribution was about 55% of the machining area. It is because the low flow rate of electrolyte at this stage is not able to take away the As shown in Figure 4, the bubble rate of electrolytes in the machining area gradually increases as the ECM proceeds. The highest bubble rate of electrolytes in the machining area reaches 18.5% at 60 s of machining time. The highest bubble rate of electrolytes in the machining area from the initial machining stage to 20 s of machining time reached 17.3%, and the electrolyte bubbles were distributed around the boundary Γ<sup>1</sup> with a minimal area. This is because the high flow rate of the electrolyte at this stage can carry away the hydrogen produced at the boundary Γ<sup>1</sup> in time. The highest bubble rate of electrolytes in the machining area increased by 1.2% from 20 s of machining time to 60 s of machining time, and the area of electrolyte bubble distribution was about 55% of the machining area. It is because the low flow rate of electrolyte at this stage is not able to take away the hydrogen produced by the boundary Γ<sup>1</sup> in time, so the hydrogen accumulates along the electrolyte flow direction, causing the electrolyte bubble distribution area to expand.

hydrogen produced by the boundary Γ<sup>1</sup> in time, so the hydrogen accumulates along the

**Figure 4.** Distribution of electrolyte bubble rate in the machining area. **Figure 4.** Distribution of electrolyte bubble rate in the machining area. **Figure 4.** Distribution of electrolyte bubble rate in the machining area.

#### *3.3. Electrolyte Flow Rate Distribution in the Machining Area 3.3. Electrolyte Flow Rate Distribution in the Machining Area 3.3. Electrolyte Flow Rate Distribution in the Machining Area*

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 8 of 16

increases, while the electrolyte flow rate remains constant.

As shown in Figure 5, the electrolyte flow rate in the machining area at the initial stage is high, up to 106 m/s. As the ECM progresses, the electrolyte flow rate gradually decreases, and the highest electrolyte flow rate in the machining area is 55.2 m/s when the machining time is 60 s. The electrolyte flow rate decreases as the machining gap As shown in Figure 5, the electrolyte flow rate in the machining area at the initial stage is high, up to 106 m/s. As the ECM progresses, the electrolyte flow rate gradually decreases, and the highest electrolyte flow rate in the machining area is 55.2 m/s when the machining time is 60 s. The electrolyte flow rate decreases as the machining gap increases, while the electrolyte flow rate remains constant. As shown in Figure 5, the electrolyte flow rate in the machining area at the initial stage is high, up to 106 m/s. As the ECM progresses, the electrolyte flow rate gradually decreases, and the highest electrolyte flow rate in the machining area is 55.2 m/s when the machining time is 60 s. The electrolyte flow rate decreases as the machining gap increases, while the electrolyte flow rate remains constant.

Initial stage The 20 s The 40 s The 60 s **Figure 5.** Distribution of electrolyte flow rate in the machining area. **Figure 5.** Distribution of electrolyte flow rate in the machining area.

### **Figure 5.** Distribution of electrolyte flow rate in the machining area. *3.4. Electrolyte Current Density Distribution and Workpiece Anode Profile Changes in the Machining Area 3.4. Electrolyte Current Density Distribution and Workpiece Anode Profile Changes in the Machining Area*

*3.4. Electrolyte Current Density Distribution and Workpiece Anode Profile Changes in the Machining Area* As shown in Figure 6a, the distance between the tool cathode and the workpiece anode gradually becomes larger. The electrolyte current density in the machining gap gradually decreases as the workpiece anode is dissolved during ECM. As shown in Figure 6b, the machining depth and height of the workpiece anode increase progressively as the ECM progresses, while the machining volume simultaneously tends As shown in Figure 6a, the distance between the tool cathode and the workpiece anode gradually becomes larger. The electrolyte current density in the machining gap gradually decreases as the workpiece anode is dissolved during ECM. As shown in Figure 6b, the machining depth and height of the workpiece anode increase progressively as the ECM progresses, while the machining volume simultaneously tends to decrease. The main reason is that as the current density decreases during ECM, the speed of ECM becomes slower, and the machining volume declines simultaneously. As shown in Figure 6a, the distance between the tool cathode and the workpiece anode gradually becomes larger. The electrolyte current density in the machining gap gradually decreases as the workpiece anode is dissolved during ECM. As shown in Figure 6b, the machining depth and height of the workpiece anode increase progressively as the ECM progresses, while the machining volume simultaneously tends to decrease. The main reason is that as the current density decreases during ECM, the speed of ECM becomes slower, and the machining volume declines simultaneously.

to decrease. The main reason is that as the current density decreases during ECM, the speed of ECM becomes slower, and the machining volume declines simultaneously.

**Figure 6.** Electrolyte current density distribution in the machining area and workpiece anode profile change curve. **Figure 6.** Electrolyte current density distribution in the machining area and workpiece anode profile change curve.

### **4. Process study of ECM of the Miniature Bearing Outer Ring 4. Process study of ECM of the Miniature Bearing Outer Ring**

### *4.1. Process Evaluation Index for ECM of the Miniature Bearing Outer Ring 4.1. Process Evaluation Index for ECM of the Miniature Bearing Outer Ring*

process evaluation index of the ECM miniature bearing outer ring is constructed:

As shown in Figure 7, it is assumed that the equation of the standard curve of the outer ring section of the miniature bearing is: *x* <sup>2</sup>+ (*y* − 1.2)2 = 4, among them: 0 ≤ *x* ≤ 0.8. The cross-sectional curve of the ECM miniature bearing outer ring does not precisely coincide with the standard curve, and there is a specific error, so the average error of the As shown in Figure 7, it is assumed that the equation of the standard curve of the outer ring section of the miniature bearing is: *x* <sup>2</sup> + (*<sup>y</sup>* <sup>−</sup> 1.2)<sup>2</sup> = 4, among them: <sup>0</sup> <sup>≤</sup> *<sup>x</sup>* <sup>≤</sup> 0.8. The cross-sectional curve of the ECM miniature bearing outer ring does not precisely coincide with the standard curve, and there is a specific error, so the average error of the process evaluation index of the ECM miniature bearing outer ring is constructed:

$$\delta = \frac{1}{n} \sum\_{1}^{n} |\mathbf{x}\_{b} - \mathbf{x}\_{a}| \tag{21}$$

Where *δ* is the average error; *x<sup>a</sup>* is the horizontal coordinate of the standard curve; *x<sup>b</sup>* is the horizontal coordinate of the machining curve; *n* is the number of points taken uniformly along the vertical coordinate. where *δ* is the average error; *x<sup>a</sup>* is the horizontal coordinate of the standard curve; *x<sup>b</sup>* is the horizontal coordinate of the machining curve; *n* is the number of points taken uniformly along the vertical coordinate.

1

*n*

**No.**

**Factor A** 

**U/V**

**Figure 7.** Standard curve and machining curve. **Figure 7.** Standard curve and machining curve.

### *4.2. ECM Miniature Bearing Outer Ring Center Composite Design Solution 4.2. ECM Miniature Bearing Outer Ring Center Composite Design Solution*

### 4.2.1. Design Solutions and Simulation Results 4.2.1. Design Solutions and Simulation Results

In this design scheme, three process parameters, namely, machining voltage, electrolyte concentration, and electrolyte inlet flow rate were selected as input quantities. One process evaluation index, namely, average error was chosen as the output quantity. The process parameters were coded with −1, 0 with 1 representing the different level values of each process parameter, where "0" represents the center point of the level value; "−1" represents the low-level value; "1" represents the high-level value. The actual and coded values of the central composite design scheme are shown in Table 2. The center composite design scheme and simulation results are shown in Table 3. In this design scheme, three process parameters, namely, machining voltage, electrolyte concentration, and electrolyte inlet flow rate were selected as input quantities. One process evaluation index, namely, average error was chosen as the output quantity. The process parameters were coded with −1, 0 with 1 representing the different level values of each process parameter, where "0" represents the center point of the level value; "−1" represents the low-level value; "1" represents the high-level value. The actual and coded values of the central composite design scheme are shown in Table 2. The center composite design scheme and simulation results are shown in Table 3.


**Table 2.** Actual and coded values of the central composite design scheme.

### **Table 3.** Center composite design scheme and simulation results. 4.2.2. Establishing the Regression Equation

**Machining Voltage Factor B Electrolyte Concentration C/% Factor C Electrolyte Inlet Flow Rate V/m/s Average Error δ/mm** According to the results given in Table 3, the regression equation between the mean error of the process evaluation index of ECM miniature bearing outer ring and the process parameters was established through multiple quadratic orthogonal regression analysis:

$$y = a\_0 + \sum\_{i=1}^{3} a\_i A\_i + \sum\_{i=1}^{3} a\_{ii} A\_i^2 + \sum\_{i$$

5 12 8 12 0.0165 6 24 8 12 0.01567 7 12 16 12 0.01635 where *y* is the output of the model; *α<sup>0</sup>* is a constant; αi is the primary term regression coefficient; *αii* is the quadratic term regression coefficient; *αij* is the interaction term regression coefficient; *γ* is the error estimate; *A<sup>i</sup>* is the coded value of the process parameters.

### 8 24 16 12 0.03715 9 12 12 9 0.01573 4.2.3. Analysis of Variance and Significance of Each Factor

 18 16 9 0.02791 18 12 6 0.01734 18 12 12 0.01767 18 12 9 0.01756 18 12 9 0.01835

4 24 16 6 0.03747

10 24 12 9 0.02558 11 18 8 9 0.01602 The regression analysis of the mean error of the outer ring of the ECM miniature bearing yielded the variance and significance analysis as shown in Table 4.


**Table 3.** Center composite design scheme and simulation results.

**Table 4.** Variance and significance of the mean error.


Where at *p* < 0.001, the factor is highly significant; at *p* < 0.05, the factor is significant; at *p* ≥ 0.05, the factor is not significant.

As shown in Table 4, the F-value of the model is 149.77 with a *p*-value less than 0.0001; the miss drafting F-value is 3.06, and the miss drafting *p*-value is 0.1226, indicating that the model is highly significant. The miss drafting term is not significant, showing that the model is meaningful and plausible. Simultaneously, R<sup>2</sup> = 0.9926, Adjusted R <sup>2</sup> = 0.9860, and Predicted R<sup>2</sup> = 0.9718. The three values are similar and less different from 1, indicating that the model fits relatively well throughout the regression region. Adequate precision = 40.3255 and much greater than 4 indicates the model is more realistic and reliable. The established regression model has a good response. However, factors C, AC, and BC are not significant, so by excluding these insignificant factors, the regression model can be optimized.

4.2.4. Optimization of Regression Models

As shown in Table 4, the impact of factors C, AC, and BC on the mean error was not significant, so the regression model was optimized using a stepwise elimination of insignificant factors to rerun the regression analysis. The variance and significance analysis of the mean error of the regression model after optimization are shown in Table 5.


**Table 5.** Variance and significance of the mean error of the optimized regression model.

Where at *p* < 0.001, the factor is highly significant; at *p* < 0.05, the factor is significant; at Podel was optimized using a stepwise.

As shown in Table 5, the F-value of the model is 290.62 with a *p*-value less than 0.0001; the miss drafting F-value is 1.92, and the miss drafting *p*-value is 0.2441, indicating that the model is highly significant. The miss drafting term is not significant, showing that the model is meaningful and plausible. Simultaneously, R<sup>2</sup> = 0.9926, Adjusted R <sup>2</sup> = 0.9892, and Predicted R<sup>2</sup> = 0.9813. The three values are similar and less different from 1, indicating that the model fits relatively well throughout the regression region. Adequate precision = 54.6205 and much greater than 4 indicates that the model is more realistic and reliable. The response of the established regression model is good. All the factors are significant at this time.

4.2.5. The Regression Equation of the Mean Error with the Normal Probability Distribution of the Residuals

According to the variance and significance analysis of the mean error of the optimized regression model, the calculation of each coefficient of the mean error regression equation was carried out, and the regression equation of the optimized mean error was obtained as:

$$\begin{aligned} \delta &= 0.0182 + 0.0050II + 0.0055C + 0.0054U \cdot C + \\ &0.0017U^2 + 0.0030C^2 - 0.0014V^2 \end{aligned} \tag{23}$$

As shown in Figure 8, the residuals of each factor are distributed around a straight line. The residuals of each factor conform to a normal distribution, so the prediction of the mean error using Equation (23) is reliable.

**Figure 8.** Normal probability distribution of residuals. **Figure 8.** Normal probability distribution of residuals. **Figure 8.** Normal probability distribution of residuals.

#### 4.2.6. Response Surface Analysis 4.2.6. Response Surface Analysis 4.2.6. Response Surface Analysis

Since the effect of electrolyte inlet flow rate on the average error of the ECM of the miniature bearing outer ring is not significant, the inlet flow rate of electrolyte is taken as 9 m/s in the response surface analysis, and the effect of machining voltage and electrolyte concentration on the average error is shown in Figure 9. Since the effect of electrolyte inlet flow rate on the average error of the ECM of the miniature bearing outer ring is not significant, the inlet flow rate of electrolyte is taken as 9 m/s in the response surface analysis, and the effect of machining voltage and electrolyte concentration on the average error is shown in Figure 9. Since the effect of electrolyte inlet flow rate on the average error of the ECM of the miniature bearing outer ring is not significant, the inlet flow rate of electrolyte is taken as 9 m/s in the response surface analysis, and the effect of machining voltage and electrolyte concentration on the average error is shown in Figure 9.

**Figure 9.** Effect of machining voltage and electrolyte concentration on the average error **Figure 9.** Effect of machining voltage and electrolyte concentration on the average error **Figure 9.** Effect of machining voltage and electrolyte concentration on the average error.

As shown in Figure 9, the average error increases with the rise of processing voltage and electrolyte concentration. This is because both the increase of processing voltage and electrolyte concentration will increase the dissolution rate of the metal material on the anode surface of the workpiece, and when the processing depth is 0.8 mm, the radius of curvature of the cross-sectional curve of the outer ring of the miniature bearing produced by electrolysis increases, resulting in a rise in the average error between the standard curve and the standard curve. As shown in Figure 9, the average error increases with the rise of processing voltage and electrolyte concentration. This is because both the increase of processing voltage and electrolyte concentration will increase the dissolution rate of the metal material on the anode surface of the workpiece, and when the processing depth is 0.8 mm, the radius of curvature of the cross-sectional curve of the outer ring of the miniature bearing produced by electrolysis increases, resulting in a rise in the average error between the standard curve and the standard curve. As shown in Figure 9, the average error increases with the rise of processing voltage and electrolyte concentration. This is because both the increase of processing voltage and electrolyte concentration will increase the dissolution rate of the metal material on the anode surface of the workpiece, and when the processing depth is 0.8 mm, the radius of curvature of the cross-sectional curve of the outer ring of the miniature bearing produced by electrolysis increases, resulting in a rise in the average error between the standard curve and the standard curve.

#### 4.2.7. Prediction of the Best Combination of Process Parameters for the Average Error 4.2.7. Prediction of the Best Combination of Process Parameters for the Average Error 4.2.7. Prediction of the Best Combination of Process Parameters for the Average Error

The optimal combination of process parameters for the average error of the ECM miniature bearing outer ring was obtained using Design-Expert 12 software analysis as shown in Figure 10. The optimal combination of process parameters for the average error of the ECM miniature bearing outer ring was obtained using Design-Expert 12 software analysis as shown in Figure 10. The optimal combination of process parameters for the average error of the ECM miniature bearing outer ring was obtained using Design-Expert 12 software analysis as shown in Figure 10.

**Figure 10.** Optimal combination of process parameters for the average error. **Figure 10.** Optimal combination of process parameters for the average error.

As shown in Figure 10, the average error of ECM of the miniature bearing outer ring can reach the minimum value of 0.014 mm under a machining voltage of 16.20 V, an electrolyte concentration of 9.29%, and an electrolyte inlet flow rate of 11.84 m/s. As shown in Figure 10, the average error of ECM of the miniature bearing outer ring can reach the minimum value of 0.014 mm under a machining voltage of 16.20 V, an electrolyte concentration of 9.29%, and an electrolyte inlet flow rate of 11.84 m/s.

### **5. Conclusions 5. Conclusions**

In this paper, a multiphysics field coupled simulation model of electric, flow, and temperature fields during ECM of the miniature bearing outer ring is established based on the gas–liquid two-phase turbulent flow model. The process of ECM of miniature bearing outer rings is investigated using a central composite design. Then some conclusions were drawn from this study, which are summarized as follows: In this paper, a multiphysics field coupled simulation model of electric, flow, and temperature fields during ECM of the miniature bearing outer ring is established based on the gas–liquid two-phase turbulent flow model. The process of ECM of miniature bearing outer rings is investigated using a central composite design. Then some conclusions were drawn from this study, which are summarized as follows:


concentration of 9.29%, and an electrolyte inlet flow rate of 11.84 m/s.

**Author Contributions:** Conceptualization, Z.L. and W.L.; methodology, Z.L. and W.L.; validation, Z.L. and W.L.; formal analysis, Z.L. and W.L.; investigation, Z.L. and W.L.; resources, Z.L. and W.L.; data curation, Z.L.; writing—original draft preparation, Z.L.; writing—review and editing, Z.L.; visualization, Z.L.; supervision, Z.L.; project administration, Z.L. and W.L.; funding acquisition, Z.L. and W.L.; examine and proposal, B.C.; All authors have read and agreed to the published version of the manuscript. **Author Contributions:** Conceptualization, Z.L. and W.L.; methodology, Z.L. and W.L.; validation, Z.L. and W.L.; formal analysis, Z.L. and W.L.; investigation, Z.L. and W.L.; resources, Z.L. and W.L.; data curation, Z.L.; writing—original draft preparation, Z.L.; writing—review and editing, Z.L.; visualization, Z.L.; supervision, Z.L.; project administration, Z.L. and W.L.; funding acquisition, Z.L. and W.L.; examine and proposal, B.C.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China, grant number 52075134.

**Data Availability Statement:** Data are contained within the article.

**Acknowledgments:** The authors would like to thank the Key Laboratory of Advanced Manufacturing Intelligent Technology of the Ministry of Education for helpful discussions on topics related to this work. The authors are grateful to Bingren Cao for his help with the preparation of figures in this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.
