Modeling of solid–liquid coupling and material removal in robotic wet polishing

In this paper, the flow characteristics of the polishing fluid between the polishing pad and the workpiece are studied for the robotic wet polishing process. The distribution of the polishing fluid radial velocity Ur and the liquid film thickness z at different rotating radii r is revealed. The two-dimensional computational domain consisting of the polishing pad surface, the workpiece wall, and the polishing fluid is established. The particle-liquid two-phase flow simulation is carried out in Fluent. The influence of different rotation rates ω of the polishing pad and different robot swing speeds v2 on the change and distribution of polishing fluid flow rate and temperature is elaborated. For polishing fluids with different average abrasive diameters dp, the position distribution of the abrasive particles in the wet polishing process is simulated; the velocity distribution of particles in the x and y directions impacting on the workpiece surface is analyzed. The three-dimensional calculation domain for wet polishing is established; the workpiece surface erosion is simulated in Fluent. Considering different combinations of polishing fluid properties Ci and polishing kinematics Pi, the material removal rate MRR and standard deviation of material removal σ on the workpiece surface are calculated. Under the same process parameters, the material removal rate test value MRRT and the standard deviation of material removal test value σT are compared with the simulated values, respectively. The results show that under the combination of 64 groups of physical parameters C1-C64 of the polishing fluids, the error between the test value (MRRT, σT) and the simulation value (MRR, σ) is within 5%. With 64 sets of polishing kinematic parameters P1-P64, the average error between the test value MRRT and the simulated value MRR is 4.19%. However, when the polishing pad rotation rate ω is high, there is an inefficient polishing area in the smaller radius from the polishing pad rotation center, which results in a lower MRRT in some tests than that in simulation, with a maximum error of 8.1%. The average error between the test value σT and the simulation value σ is 3.77%. When the pressure P of the polishing pad is high, the large particles embedded in the polishing pad surface follow its rotation, causing deep scratches on the workpiece surface, which results in a larger σT in some tests, with a maximum error of 7.8%. In conclusion, the material removal principle and the influence of different process parameters in the robotic wet polishing process are revealed in this paper through modeling and simulation of the particle-liquid two-phase flow, giving an accurate estimation of the material removal rate of the robotic wet polishing process.


Introduction
The development of science and technology in the fields of aerospace, national defense, and medical brings a higher standard for the surface quality of critical parts which requires precision polishing to achieve the required surface finish [1,2]. Wet polishing uses polishing fluid as the polishing medium. During the polishing process, the polishing pads are used to drive the fluids and the abrasive particles to polish the surface of the workpiece and realize material removal. It can effectively reduce the temperature of polishing operations, so as to reduce the surface burns and deformation of the workpiece. It is also helpful for achieving stable surface quality of the workpiece and reducing the harm to human body due to the polishing dust [3,4]. Therefore, wet polishing is considered an important processing technology to improve the dimensional accuracy and reduce the surface roughness of components.
Currently, different scholars have constructed empirical formulas and established mathematical models through experimental data. Some of them also use artificial intelligence to explore the principle of material removal in wet polishing process.
Yan Y [5] proposed a material removal model with the Preston equation and verified the effectiveness of the model experimentally with single-variable tests. Orthogonal experiments were conducted on process parameters of rough and fine polishing, material removal depth and surface roughness were set as evaluation indexes, and the better combination of process parameters was selected according to the experimental results.
Zhang P [6] et al. build an accurate mathematical model of orthogonal rotary regression test of tri-factor quadratic of MRR. The model about double-faced mechanical polishing was established to obtain the relationship between any point on SiC substrate and polishing pads. In order to improve the material removal rate (MRR), on the premise of ensuring the surface roughness requirements, the best optimized material removal rate parameters were obtained. Liu S [7] studied the wet polishing method of carbide tools and established a contact model between the tool, abrasive particles, and polishing pads considering the physical properties and motions of the polishing pads and tool. Besides, they proposed a surface roughness model of the carbide tool from the statistics and tribological aspects. The influence of particle volume, polishing pressure, and particle diameter on surface roughness of carbide tools is studied in their research by experiments. The variation curves of surface roughness with polishing parameters were obtained in high agreement with the theoretical prediction curves.
Fan S [8] proposed a BP-NCABC algorithm. They established a model between the components of polishing fluid, polishing process parameters, polishing efficiency, and within-sheet non-uniformity. Based on the proposed approach, the polishing efficiency and within-sheet nonuniformity are well predicted and optimized.
In the above researches, it is assumed that the particle diameters were consistent and the distribution was uniform when modeling the polishing process. The dynamic effects of liquid flow and particle collision on the trajectory and velocity of abrasive particles in wet abrasive polishing operations were not considered. The transient variability of the abrasive particle trajectory and velocity as well as the uneven distribution of the abrasive particle position and diameter at different times was not considered. It leads to a large deviation between the theoretical model and the actual results of material removal. The method of predicting material removal characteristics by combining experimental data with artificial intelligence algorithm required a large amount of data to build empirical formulas. It cannot explain the influence of wet polishing process parameters on workpiece surface material removal from the polishing principle side. In addition, when the polishing object or processing conditions change, the prediction model derived from the original experiment is no longer applicable.
The wet polishing is a multi-input, multi-output kinetic process. In addition to the controllable processing kinematic parameters such as the rotation rate, pressure, and trajectory of polishing pad, there are other factors such as the material type, speed, trajectory, and volume fraction of abrasive particles, as well as the flow property and temperature change of the polishing fluid. These factors will affect the final surface finish in polishing process significantly [9]. The diameter of the particles is microscopic; and the number of particles in the polishing fluid is too huge to count. It is difficult to record the particles' physical information and the interaction with liquids during the polishing process with available test equipment.
In recent years, with the development in computer software and hardware, numerical simulation has become a third scientific research method alongside with experimental and theoretical methods [10]. Polishing fluid mainly consists with abrasive particles and liquids, which is a solid-liquid two-phase mixture fluid. The research method of computational fluid dynamics for solid-liquid two-phase fluids allows numerical modeling of some difficult-to-observe quantities. The results obtained by this method are highly consistent with the actual results and have perfect repeatability [11][12][13]. Relevant work has been conducted by different researchers.
Nguyeu N [14] et al. explored the distribution of the abrasive particles in the polishing fluid using the discrete phase model in computational fluid dynamics combined with finite element simulation. The results showed that with the increased distance between the polishing pad and the workpiece surface, and with the increased rotation rate of the polishing pad or the workpiece, the number of abrasive particles in the gap between the workpiece and the polishing pad tended to increase. However, the distribution of particles in the gap is not uniform. Runnels S [15] et al. are the first to analyze the flow behavior of polishing fluid in wet polishing. They solved the steady three-dimensional Navier-Stokes equation with the finite element method and calculated the liquid film thickness and the angle for attack condition of the workpiece in the wet polishing process. The limitation is that many mandatory assumptions were made in the analysis process for this mathematical model. For example, both the workpiece and the polishing pad were rigid bodies with smooth surfaces. The rotation of the workpiece relative to the polishing pad was not considered; the workpiece surface was regarded as a spherical surface with large radius. Therefore, the liquid thickness can only be estimated qualitatively.
Cho C [16] et al. analyzed the thickness, pressure distribution, and contact shear stresses of the polishing fluid. The results showed that the normal stress on the workpiece surface was more uniform when the load is small while the rotation rate and viscosity ratio are large. Besides, this leads to a reduced edge shear stress. However, it assumed that the surface of the polishing pad was smooth, so the results were only applicable to qualitative prediction.
Hua C [17] deduced the flow state of the polishing fluid under the inviscid assumption by using the hydrodynamic equation. The flow behavior of the polishing fluid is analyzed by using the finite element method. Finally, he verified the correctness of the simulation results by using the particle image velocimetry experiment. Combined with the finite element simulation of multi-phase flow model and particle image velocimetry experiment, the flow characteristics of the polishing fluid on the polishing surface were explored. However, it assumed that the liquid is nonviscous and the polishing pad surface is smooth. Only the flow characteristics of the polishing fluid were explored, so there was an error between the model and the actual conditions.
The above computational fluid dynamics (CFD) researches only explored the flow characteristics of the fluid or the distribution of abrasive particles. The material removal rate for wet polishing process cannot be directly simulated in the above methods. Therefore, the effectiveness of their proposed methods cannot be verified experimentally. Besides, the properties of the polishing pad and the polishing fluid were not considered sufficiently in the calculation process, leading to a large deviation in the simulation results compared with the actual case.
The above researches have revealed the material removal principle of wet polishing to a certain extent. However, the theoretical models and simulation models proposed in the above researches are too idealistic to predict the material removal rate accurately. There are few studies that fully consider the interaction among the fluids, abrasive particles, the surfaces of polishing pad, and the workpiece. Besides, the theoretical relation between different process parameters and material removal theory is hardly revealed. In most researches, the polishing pad is simplified as a smooth surface, without considering the influence of the polishing pad roughness on the fluid flow and the movement of the abrasive particles. In this paper, based on the fluid motion equation, the flow characteristics of the polishing fluid between the polishing pad and the workpiece are studied to reveal the distribution of the polishing fluid radial velocity and the liquid film thickness at different rotating radius. The two-dimensional computational domain consisting of the polishing pad surface, workpiece surface, and the polishing fluid is established. Based on the discrete phase model (DPM) and the volume of fluid (VOF), the two-phase flow model of liquid and abrasive particles is developed to simulate the velocity and temperature distribution of the liquid in the computational domain, as well as the transient position and velocity distribution of abrasive particles. The threedimensional computational domain consisting of workpiece surface, polishing pad surface, and polishing fluid is build. Then, the erosion state of workpiece surface is simulated based on computational fluid dynamics material removal rate model. The influence of the different process parameters on the material removal rate and the standard deviation of material removal of the workpiece surface are clarified. Through the material removal test on the workpiece surface, the effectiveness of the proposed solid-liquid coupled twophase flow simulation model is verified, providing a theoretical basis for wet polishing operations.

Methodology
In order to model the solid-liquid flow in wet polishing process and study the material removal principle, the methodology of this research is introduced in Fig. 1, and can be explained from the following four aspects.
1) The flow characteristics of polishing fluid in the process of robotic wet polishing are first described from the macro-point of view. Then, considering the polishing pad motion condition at the boundary, the relation between the radial velocity u r of the polishing fluid and the radius of rotation r at typical working condition is analyzed with N-S equation. Furthermore, the relationship between the polishing fluid liquid film thickness z and the rotation radius r under typical working conditions is obtained with the continuity equation. 2) In the two-dimensional computational domain consisting of the workpiece surface, polishing pad surface, and the polishing fluid, the two-phase flow model of liquid and abrasive particles is established to simulate and elaborate the influence of wet polishing process parameters on the solid-liquid coupling characteristics. The twodimensional fractal geometry W-M model is selected to describe the micro-profile of the polishing pad roughness. Then, the influence of different polishing pad rotation rates ω, robot swing speed v 2 on the velocity distribution, and temperature change of polishing fluid is simulated based on VOF and DPM models in Fluent. At the same time, for different mean abrasive particle diameters d p , the position and the velocity distributions of the abrasive particles on the workpiece surface in the x and y directions were simulated and analyzed for the robotic wet polishing process.
3) The three-dimensional computational domain consisting of workpiece surface, polishing pad surface, and polishing fluid is constructed; and the erosion of workpiece surface is simulated. Besides, the influence of the polishing fluid parameters C i and the polishing kinematic parameters P i on the material removal rate MRR and the standard deviation of the material removal σ are analyzed with finite element software. 4) Experiments are carried out to further verify the influence of the parameters C i and P i on MRR and σ. Sixtyfour different sets of C i and P i are considered in the experiments. The experimentally obtained MRR and σ are compared with the values calculated from simulation to prove the effectiveness of our proposed modeling methods.

Radial velocity analysis of the polishing fluids
Industrial robots have the advantages of large moving range, flexible configuration, and low cost comparing to a 5-axis machine tool; and they gradually become an important equipment in wet polishing [18]. In the wet polishing process, the polishing pad adheres to the end of the polishing spindle. The normal force and tangential force in the polishing process are implemented on the workpiece through the polishing fluid under the polishing pad area. Also, the waste and heat generated during the polishing process are took away from the working area by the polishing fluid when it leaves the polishing pad. Therefore, the motion of the polishing fluid plays an important role in wet polishing process [19,20]. In this paper, the motion of the polishing fluid between the polishing pad and the workpiece is theoretically deduced to reveal the flow characteristics. Polishing fluid is a kind of fluid which has inconspicuous viscosity-temperature and viscosity-pressure properties. The properties of the polishing fluid can be considered not changed over time throughout the machining process. According to the Navier-Stokes equation [21], it is known that: where ρ denotes the density of the in-compressible fluid; V represents the absolute motion velocity of the micro-element; t is time; F represents the body force on the microelement (in this paper, F mainly refers to gravity); ▽ is the gradient operator; P represents the pressure on the microelement; µ is the viscosity of the polishing fluid; and ▽ 2 is the Laplace operator. The left side of the equation represents where u, v, and w are the velocity components of the fluid at the point (x, y, z) at time t; and F x , F y , and F z are the components of the body force on the micro-element in the directions of x, y, and z. From Eq. (2), it can be concluded that the change of the momentum of the fluid is mainly caused by the body force and pressure on it.
In robotic wet polishing, the polishing fluid is injected from the center hole of the polishing pad; and it flows out of the polishing pad area due to the centrifugal force while the polishing pad rotates. The velocity of the polishing fluid along the radial direction of the polishing pad determines the motion of the polishing fluid and the thickness of the liquid film. In order to better express this fluid motion, the u component of the polishing fluid in Eq. (2) is expressed in a scalar form in the cylindrical coordinate system as: where r, θ, and z are the radial, circumferential, and axial coordinates in the cylindrical coordinate system, respectively; u r , u θ , and u z are the radial, circumferential, and axial components of the fluid velocity, respectively; and F r is the body force in the radial direction.
In Eq. (3), △ u r = ∂ 2 u r/ ∂r 2 + ∂u r/ ∂r*r −1 + ∂ 2 u r/ ∂θ 2 *r −2 + ∂ 2 u r/ ∂z 2 . Because the fluid flows from the center to the edge of the polishing pad evenly, the circumferential velocity of the fluid at the same radius is considered to be the same at different angles. Therefore, ∂u r /∂θ = 0. Since the thickness of the liquid film is much smaller than the diameter of the polishing pad, the change of the radial velocity of the fluid in the axial direction is negligible (∂u r /∂z = 0). Besides, the rotational velocity of polishing pad is a constant during the polishing process, which leads to ∂u θ /∂θ = 0. When the flow field is stable, the parameters in the flow field do not change with time, and the material derivative of the flow field is zero. So the left term of Eq. (3) is zero. Since there is no gravity acting in the radial direction, F r is zero. Therefore, Eq. (3) can be simplified as: Figure 2 shows the motion of the polishing pad during robotic wet polishing process. The rotation rate ω of the polishing pad is controlled by a spindle, while the trajectory of the polishing pad is controlled by the robot. The polishing pad moves in multiple "V" shape pattern on the workpiece surface, as shown in Fig. 2. The translational velocity of the polishing pad can be split into two components, which are swing speed v 2 and step speed v 1 .
During the wet polishing process, polishing fluid flows in the micro-channel formed by the workpiece and the polishing pad. At the edge of the polished pad, r = r p , the velocities of the polishing fluid can be expressed as follows.
where r p is the radius of the polishing pad; v q is the flow rate of the polishing fluid; and u rp , u θp , and u zp are the radial, circumferential, and axial velocities at the edge of the polishing pad.
The polishing fluid is a type of Newtonian fluid, so the viscosity µ is a constant. The following equations can be obtained from Eq. (4) by integration of r: After the quadratic integration, we can get: The boundary condition is obtained from Eq. (5). When r = r p , the radial velocity component of the liquid film at the edge of the polishing pad is u rp = v 1 + v 2 + v q + ω rp . Substituting Eq. (5) into Eq. (7), the integral constant C can be calculated. Substituting the integral constant C back into Eq. (7), In order to show the variation rules of the radial velocity u r in Eq. (8) intuitively, the relationship curves between the radial flow rate u r of the polishing fluid and the rotation radius r are plotted as shown in Fig. 3. Multiple sets of parameters (as given in Table 1) in robotic wet polishing are considered in the plot. The viscosity of polishing fluids is measured by the viscometer of three typical alumina polishing fluids, and the average value is taken as the final viscosity µ. Figure 3 shows how the radial velocity u r of the polishing fluid changes with the rotation rate of the polishing pad ω, the swing speed v 2 , and the flow rate v q at different rotation radii r. It can be seen that the changes of the rotation rate ω affect the relation between the radial velocity of the polishing fluid u r and the rotation radius r more comparing to the changes of the swing speed v 2 and the flow rate v q . This is because that the centrifugal force loaded on the polishing fluid is proportional to rotation rate ω and the square of the rotation radius r 2 .

Liquid film thickness analysis of the polishing fluids
In the robotic polishing process, the polishing fluid is brought into the polishing area on the workpiece surface through the center hole of the rotating polishing pad (Fig. 4). The polishing fluid between the workpiece and the polishing pad is flowing with a certain speed during the wet polishing process. Due to the wedge effect in the hydrodynamics theory [22], when the system reaches equilibrium, the workpiece and the polishing pad are not parallel. There is a very small tilt angle in between, called hydrodynamic pressure inclination angle HPI. In this case, the liquid film thickness of polishing fluid is thicker at the inlet and is thinner at the outlet. Due to the existence of HPI, the polishing fluid can generate a stable hydrodynamic pressure to withstand the external load while flowing circularly. The thickness of the polishing liquid film between the workpiece and the polishing pad at different radii maintains stable when the wet polishing process is in a steady state. At this time, the load on the polishing pad is equal to the load on the polishing liquid film. Figure 4 shows the flow diagram of the polishing fluid during the robotic wet polishing process. It is assumed that during the wet polishing process, the polishing fluid  distributed between the polishing pad and the workpiece is in a thin liquid film shape, with a thickness z at different positions. The distribution of the radial velocity u r along the radius of rotation r has been obtained in Eq. (8). Considering the relation between the flow rate and the cross-sectional area explained in the continuity equation, the relation between the liquid film thickness z and the radial velocity u r can be obtained: where Q is the flow volume of the polishing fluid; and z is the liquid film thickness. Substituting Eq. (8) into Eq. (9), we have: Several typical working conditions are simulated and the results are shown in Fig. 5. The thickness of the polishing fluid decreases rapidly from the center to the periphery in the area where the radius is less than 50 mm. The decreasing speed slows down in the area where the radius is between 50 and 125 mm. As can be seen from Fig. 5a, when the rotation rate ω of the polishing pad increases, the distribution of the film thickness z becomes more nonuniform along the radial directions, with a larger thickness difference between the center and the edge position.  Fig. 5b, we can see that the overall thickness of the liquid film between the polishing pad and the workpiece increases. This is because that the liquid flow from the center to the edge of the polishing pad increases when the flow rate v q goes up.

Establishment of polishing pad surface micro-profile model
The polishing pad commonly used in industry is polyurethane or non-woven material. Its surface is a rough area composed of many irregular and complex micro-unit [23]. W-M model is a combination of deterministic periodic function and random function, which has anisotropy. Its scale invariance overcomes the shortcomings of the traditional statistical sense, which the geometric length changes with the measurement scale. It well reflects the uneven morphology and complex structure of the area to model the rough surface [24]. The expression of two-dimensional fractal geometry W-M model is: In the equation, D is the fractal dimension. The geometric plane is a combination of countless lines, and the figure formed in the whole process is the fractal dimension between 1 and 2. G is the characteristic scale coefficient, and it reflects the magnitude of z(x) and determines the dimensional magnitude as well as the final size of z(x). γ is a parameter representing the spectral density of the surface profile, which determines the self-similarity of the fractal and the spectral density. It is related to the spatial frequency of the profile curve, usually greater than 1. n 1 is the starting frequency corresponding to the lowest frequency in the surface profile curve, which is related to the sampling length L of the profile, n 1 = 1/L. n u is the upper limit of high frequency, n u = 1/2 δ, among δ is the resolution of the contour.
From Eq. (11), it can be seen that surface roughness is a non-stationary random process. As γ is determined, three parameters (D, G, L) determine the shape of the surface. Refer to some parameters of polishing pad in the paper [25][26][27][28], and take G = 1.2 × 10 −11 m, D = 1.61, γ = 1.5, and L = 4.8 × 10 −4 m. By substituting these parameters into Eq. (11), we can get the micro-profile of the polishing pad cut surface as shown in Fig. 6. It can be seen that the figure contains many rough structures with different ranges. According to Sum's research results [29], the measured roughness profile height and frequency are consistent with the calculated values. It is a surface-following method for the modeling of two and more immiscible fluids under a fixed Eulerian grid. In the VOF model, to track each mesh by fluid flow, the concept of fluid volume function, and its control equations are introduced. The solid-liquid two-phase fluids are simulated by solving separate momentum equations and treating the volume fraction of each fluid that passes through the region [30,31]. For the VOF model, the tracking of the interface between phases is performed by solving the continuity equation for the volume fraction. For the q phase, the volume fraction equation is: In this equation, a q is the volume fraction of phase q, v q is the velocity of phase q, ρ q is the density of phase q, and S aq is the source item of a q , with default value 0. m pq is the mass moved from phase p to phase q, and m qp is the mass moved from phase q to phase p. For the initial phase (liquid phase in this paper), the volume fraction expression is as follows: α q values range between the interval [0,1], and according to the values of α q , there are three different conditions: 1) α q = 1, indicating that the cell is completely occupied by solids 2) α q = 0, indicating that the unit is completely occupied by liquid 3) 0 < α q < 1, indicating that the unit is partly liquid and partly solid, and solid-liquid interface exists.
The physical parameters of the mixed fluid in each cell are obtained by calculating the volume fraction of the physical parameters of each phase and then getting the weighted average value. In the mixed fluid computational cell, the density and viscosity are computed as follows.
According to the continuous surface force model proposed by Brackbill [32], the surface tension is explained as a continuous three-dimensional effect acrossing the interface. The effect of surface tension is simulated by the addition of a source term to the momentum equation. The divergence theorem can be expressed as a volume force and added to the momentum equation as a source term; its form is as follows: In Eq. (16), the surface curvature w = ▽·n/|n|, and the interface normal vector n = ▽α q . α ij is the surface tension coefficient between adjacent liquid levels phase i and j, N/m. α i and α j are the volume fractions of phase i and j; ρ i and ρ j are the densities of phase i and j, kg/m 3 ; and w i and w j are the surface curvatures of phase i and j.

Establishment of DPM
One of the most widely used methods for describing the motion of the fluid with particles is the discrete phase model (DPM). In the DPM, the particles are treated as a discrete system. It is first to calculate the continuous phase flow field. Then, the force on each particle is solved with considering the flow variables. As a result, we can obtain the particle velocity and track the trajectory of each particle. In the robotic wet polishing operation, the abrasive particles in the polishing fluid are equivalent to the discrete phase and the liquid is equivalent to the continuous phase. The DPM is used to solve the coupling problem between the abrasive particles and the liquid, which can track the trajectory as well as the force of each abrasive particle. The trajectory of discrete phase abrasive particles is predicted by differential calculation of the force in the Lagrangian coordinate system [33]. The equation in the Cartesian coordinate system is: In this equation, u is the velocity of the fluid, u p is the velocity of the abrasive particles, ρ is the density of the fluid, ρ p is the density of the abrasive particles, and g is the gravitational acceleration. F is an additional force, containing the Basset force, Saffman lifting force, and Magnus lifting force [34].
In Eq. (17), F D (u-u p ) is the traction force of per unit mass in the abrasive particle, for abrasive fluids under sub-scale conventional: In this equation, d p is the mean diameter of abrasive particles and Re is the relative Reynolds number, calculated as: C D is the traction coefficient, which can be obtained from the following equation: Based on the discrete time step, Eq. (17) is solved by integral operation to obtain the particle velocity at each position on the particle track. On the solid boundary, for abrasive particles of larger diameter, their motion is primarily influenced by the inertia and collision between abrasive particles and wall. Abrasive particles of smaller diameter at the boundary are primarily influenced by turbulence.

Establishment of two-dimensional computational domain model for wet polishing process
In the schematic diagram of wet polishing operation presented in Fig. 1, we choose a micro-area which has the Re 2 same velocity at a distance r/2 from the rotation center of the polishing pad, and simplify the two dimensions' computational domain to this micro-area. For the microscopic outline of the polishing pad cut surface as shown in Fig. 6, it is smoothed in MATLAB, and the generated curve is imported into SolidWorks to create the two dimensions' surface. It results in a relatively smooth two-dimensional computational domain shape as shown in Fig. 7. In Fig. 7, the top wall of the computational domain in two dimensions is the polishing plate wall; the bottom wall is the workpiece wall; the left wall acts as an inlet for the polishing fluid; and the right wall serves as the outlet. The two-dimensional equivalent velocity of the polishing wall relative to the fluid is v p ; the two-dimensional equivalent velocity of the workpiece wall relative to the fluid is v g ; the length of the computational domain in two dimensions is L 1 ; the polishing pad surface valley height is h a ; and the surface peak height is h b . The two-dimensional computational domain surface created within SolidWorks is imported into Fluent. The mesh processing is carried out on the 2D computational domain using the mesh module. The expansion layer is added to the polishing pad wall in order to simulate the near-wall flow gradient more completely. In the model setting, the VOF model is used for multi-phase flow to track the interface between polishing fluid and abrasive particles. The primary item is set as water, and the secondary item is set as aluminum oxide. The RNG k-ε model [35] is selected for the viscosity model. It provides an analytic formulation of the flow viscosity considering low Reynolds number flow based on the standard k-ε model, which can handle the flow better with higher strain rate and curved flow lines. It is more accurate for the fluid in the region close to the wall.
The motion of abrasive particles is tracked with the DPM, and the motion of the abrasive particles follows Newton equation as well as the Stokes-Cunningham resistance law. The entrance of the two-dimensional computational domain is defined as the jet source. The physical parameters of the abrasive particles are set. When the boundary conditions were set, the fluid and abrasive particles flow in from the inlet at the same rate. To simulate the effect of heat exchange between the polishing fluid and the polishing pad as well as the workpiece surface, the energy equation is used, and the heat transfer method is set for the up and bottom walls separately. The SIMPLE algorithm was used for the coupling  between pressure and velocity in the solution configuration. The PRESTO and second-order upwind are used to conduct spatial dispersion of pressure and velocity, respectively [36]. Depending on the typical working conditions and capability of robotic wet polishing auxiliary equipment, the main parameters of the solid-liquid coupled two-phase flow model are set as given in Table 2.

Calculation results of solid-liquid two-phase fluid model
The rotating rate of the polishing pad ω is selected to be 100 rpm, 300 rpm, and 600 rpm, respectively. The swing speed of the robot v 2 is chosen as 100 mm/s, 300 mm/s, and 600 mm/s. The step speed v 1 is 10 mm/s. Then, the equivalent velocities of the polishing pad wall surface relative to the polishing fluid in the two-dimensional computational domain v p are 104.4 mm/s, 312.5 mm/s, and 625 mm/s, respectively. The equivalent velocities of the workpiece wall relative to the polishing fluid v g are 100.5 mm/s, 300.17 mm/s, and 600.08 mm/s, respectively. The other selected parameters are listed in Table 2 by default. The flow field simulations under different working conditions are performed separately, and the contours of velocity and temperature distribution of the flow field are obtained as shown in Fig. 8. Based on the velocity cloud diagram of the solid-liquid mixture flow shown in Fig. 8, in the area close to the wall of workpiece and polishing pad, the velocity of the fluid and the wall surface is basically consistent, which is in agreement with lubrication theory [37]. Due to the existence of the rough peak, there is a low-speed area within the rough peak, where the fluid movement is mainly affected by the inertia force. The common polishing pads have the adsorption effect on the liquid. It is easy for this region becoming the agglomeration area of impurity particles and abrasive particles. It makes the surface of the polishing pads becoming dirty, reducing the ability as carriers of polishing fluid, and finally reducing the polishing efficiency in actual robotic wet polishing operations. Agglomerated impurity particles or abrasive particles drop out during polishing, making it easy to leave a large scratch on the workpiece surface, which affects the quality of the workpiece polishing. By comparing the results of a 1 , a 2 , and a 3 in Fig. 8, it can be seen that increasing the rotation rate of the polishing pad ω results in an improvement in the ability to disrupt the liquid. The fluid molecules of the polishing fluid and the debris shed from the workpiece surface are primarily affected by the traction of the fluid, and the low-velocity regions within the rough peak are reduced.
Circulating and renewing the polishing fluid helps to take away excess heat, reduce the deformation, or change in material properties of thin-walled workpiece. Comparing the results of b 1 , b 2 , and b 3 in Fig. 8, with the increase in the rotation rate of the polishing pad ω, the work and heat generated in the wet polishing process increase. The flow rate of the polishing fluid molecules increases under the effect of centrifugal force. Although the renewal of the polishing fluid with cooling effect on the per unit area is accelerated, the temperature on the surface of the polishing pad and the workpiece is still gradually rising. Figure 9 shows the temperature changes around the peaks on the polishing pad surface and at the center of the liquid film as well as on the workpiece surface in the two-dimensional calculation domain. Assuming that the temperature of the renewed polishing fluid T at the inlet of the 2D computational domain is 300 k, due to the heat exchange among the polishing fluid, the polishing pad, and the workpiece, the temperature of the polishing fluid gradually increases when it flows to the outlet of the 2D calculation domain.
Any point at the inlet can be taken as a starting point for abrasive particles. The mean abrasive particle diameter d p is set at 1 µm and 4 µm, while the molecular velocity of polishing fluid v q is set at 0.3 m/s. Figure 10 shows that, because of the presence of turbulence in the flow field, the abrasive particles move with the fluid. Mutual collisions take place between the abrasive particles. Comparing Fig. 10a and b, for the same flow speed v q , the majority of abrasive particles with a mean diameter of 1 µm are suspended in the middle or top position of the 2D computational domain, whereas more abrasive particles with a mean diameter of 4 µm impact on the workpiece surface. Comparing Fig. 10c and d, it can be seen that most of the impacting abrasive particles at the workpiece surface are larger than the mean diameter. Abrasive particles with large diameters impact the workpiece surface more rapidly due to the influence of gravity and high inertia. It is difficult to alter the trend of motion after collision. The fluid traction force is primarily applied to small abrasive particles. So when the flow motion is more violent, it is easier to change one's own trajectory after collision, and impact at the workpiece surface with a slower velocity. These abrasive particles impacting the workpiece surface have high kinetic energy under the action of fluid turbulence, polishing pad pressure P and centrifugal force, which form material removal on the workpiece surface. Therefore, polishing fluids containing abrasive particles with a mean diameter of 4 µm have a higher material removal efficiency than polishing fluids with a mean diameter of 1 µm under the same conditions. On the other hand, when using a polishing fluid with a smaller abrasive particle size for robotic wet polishing applications, the flow rate and rotation rate of the polishing fluid can be appropriately reduced to wait until more abrasive particles can impact the workpiece surface so as to improve the efficiency of material removal.
During the robotic wet polishing process, the abrasive particles mechanically impact the uneven area of the workpiece 1 3 at a certain angle, which cause material removal from the surface of the workpiece. The direction of motion and size of the abrasive particles significantly affect the removal of material from the workpiece surface. The mean abrasive particle diameter of 4 µm is chosen, and the polishing speed is set to be 100 rpm, 300 rpm, and 600 rpm. Figure 11a shows the speed of the abrasive particle in the y direction when it reaches the surface of the workpiece. It can be seen that, upon entering the 2D computational domain, since the abrasive particles do not reach the workpiece surface, the abrasive velocity at the workpiece surface is zero. For the abrasive particles that hit the workpiece surface earlier, the kinetic energy is primarily affected by gravity. The angle between the velocity direction and the workpiece surface is greater, so that the velocity in the y direction is high. The influence of turbulence on the abrasive particle is not obvious at the initial stage where the particles are flowing from the inlet to the outlet in the 2D computational domain. As the angle between velocity and workpiece area decreases, for a short period of time, the particle velocity decreases slightly in the y direction. Turbulence significantly increases the kinetic energy of abrasive particles located in the middle area of the 2D computational domain. Though the angle between velocity and workpiece area becomes smaller, there is a temporary tendency for the velocity in the y direction to increase. The effect of turbulence on the abrasive particles has been stabilized in the rear region of the 2D computational domain. The speed in the y direction shows a decreasing trend as the angle between speed and workpiece area continues to decrease.
The velocity variation of the abrasive particles in the x direction is shown in Fig. 11b. The particles' velocity at the entrance is primarily determined by the velocity of the jet port. During the abrasive particles leading to the outlet, the kinetic energy increases due to turbulence, and the angles of abrasive particles impacting on the workpiece surface become smaller. As a result, the abrasive particles' velocity in the x direction increases gradually. The growth rate of abrasive particles in the x direction slows down when the turbulence gets stable during the middle and end of the 2D computational domain. As can be observed in Fig. 11, increasing the polishing rotation rate ω causes the abrasive particles to have a higher kinetic energy. The abrasive particles impact the workpiece surface at a certain angle, and the velocity in both x and y directions shows a general tendency to increase. As the abrasive particles continuously collide in the process of moving with the fluid, some of them lose kinetic energy, while some of them change the

Simulation of material removal
In order to explore the material removal rate under different working conditions, the following material removal model is used to simulate the material removal under different process parameters.
where R is the amount of erosion wear; m p is the mass flow rate of abrasive particles; C(d p ) is the particle diameter function; α is the incidence angle of abrasive particles colliding with the workpiece surface; f(α) is the impact angle function; v is the velocity of abrasive particles relative to the workpiece; b(v) is the velocity function of abrasive particles; and A p is the surface area of abrasive particles eroding the workpiece. The material removal rate MRR of the workpiece surface is defined as: where A is the total grid area, A i is the area of a single grid, t is the polishing time, and R i is the amount of material removed within a single grid area.
The standard deviation of material removal σ is defined as: where N is the overall number of cases. The 2D computational domain illustrated in Fig. 5 is transferred to a three-dimensional (3D) computational domain with a thickness of 30 µm. The solid-liquid coupled two-phase flow simulation is carried out in Fluent to analyze the erosion contours on the workpiece surface with multiple sets of typical working conditions. From Fig. 12, it can be seen that the mean material removal is high in some typical working conditions while the distribution of material removal is uniform in some other working conditions.
The physical properties of the polishing fluid have a significant influence on material removal. The MRR and the standard deviation of material removal σ are computed considering different parameters. The mean diameter d p of the abrasive particles is set to 1 µm, 4 µm, 7 µm, and 15 µm, respectively, while the volume fraction f p of the abrasive particles is set to 5 wt%, 10 wt%, 15 wt%, and 20 wt%, respectively. The flow rate v q of the polishing fluid at the inlet is set to 0.1 m/s, 0.2 m/s, 0.3 m/s, and 0.4 m/s. As shown in Fig. 13a, assuming the flow rate v q of the polishing fluid is a constant, the MRR starts increasing for a while, and then turns to decreasing when the mean abrasive particle diameter d p and abrasive particle volume fraction f p keep increasing. This is because that abrasive particles are more likely to impact the workpiece surface through gravity and inertia when d p increases, causing more material removal. Besides, with the increase of d p , the kinetic energy of each single abrasive particle is also increased, so the material removal Variation curve of the velocity in the y-direction with position when the abrasive particle reaches the surface of the workpiece (b) Variation curve of the velocity in the y-direction with position when the abrasive particle reaches the surface of the workpiece Fig. 11 Relation curve between abrasive particle velocity and position on the workpiece surface R i of single abrasive particle increases. However, if d p is too large, the abrasive particles will settle to the workpiece surface too quickly. In this case, the main kinetic energy of the abrasive particles becomes the initial kinetic energy. As a result, the effect of turbulence on the abrasive particles kinetic energy is weakened; and the material removal R i of single abrasive particle is weakened. Therefore, the MRR will first increase and then slightly decrease with the increase of d p . On the other hand, when d p and v q of the polishing fluid are constant, with the increase of f p , the number of abrasive particles in contact with the workpiece surface increases, and MRR is improved to a certain extent. However, exorbitant f p increases the collision between abrasive particles, causing a certain amount of lose in their kinetic energy. It will also lead to the viscosity increase of the polishing fluid, which reduces the effect of turbulence on the kinetic energy of the abrasive particles, causing MRR to decrease within a certain range. With the d p and the f p unchanged, the flow rate v q of the polishing fluid increased, the probability of the abrasive particles contacting the workpiece surface decreases. Although the kinetic energy of the abrasive particles increases, the included angle between the velocity direction and the workpiece surface is small; and the polishing normal force is weak. Therefore, with the increase of the v q , the MRR of the workpiece surface decreases slightly. However, it can be seen from Fig. 13b that the increasing of v q can increase the tangential force of polishing, causing reduction of the microscopic surface unevenness on workpiece. It leads to a small decrease in the σ. The σ increases as d p and f p increase. This is because with the increase of the kinetic energy of the abrasive particles, the non-uniformity of R i at different regions of the workpiece surface is amplified. The increase in f p will also improve the agglomeration of abrasive particles, which results in large removal in local areas of the workpiece surface, and finally improves σ.
The kinematic parameters are also important to the material removal of robotic wet polishing. Therefore, the MRR and σ of the wet polishing process are calculated for different kinematic parameters. The rotation rate ω in robotic wet polishing is set to 100 rpm, 300 rpm, 500 rpm, and 700 rpm, respectively, while the robot swing speed v 2 is set to 100 mm/s, 200 mm/s, 400 mm/s, and 600 mm/s, respectively. The polishing pad pressure P is set to 30N, 100 N, 300 N, and 500 N. From Fig. 14a, it can be seen that the increase of ω and P both can improve MRR, while P has a stronger effect on the increase of MRR compared ω. From Fig. 14b, it can be seen that σ goes up while P increases. This is because that the normal polishing force of the abrasive particles increases with a higher P. Therefore, the non-uniformity of the particles impacting on the workpiece surface rises, leading to an increase in σ. However, the increasing rate of MRR slows down when pressure P keeps going up. The improvement of ω primarily increases the tangential polishing kinetic energy of the abrasive particles, causing elimination of micro-unevenness and bulge on the workpiece surface, therefore, reduces σ. From Fig. 14a and b, it can be seen that the change of the robot swing speed v 2 does not have a significant effect on MRR and σ.

Experiments
Experiments are carried out to further verify the influence of the parameters C i and P i on MRR and σ. The experimental setup is shown in Fig. 15. An industrial robot equipped with a force-controlled polishing spindle is used to control the trajectory, pressure P, and rotation rate ω in wet polishing process. The force-controlled polishing spindle mainly consists of a servo motor, a force control device, a floating spindle, and a motion transmission mechanism. The floating spindle can rotate and move axially at the same time; and the servo motor is used to drive the floating spindle to rotate at the expected velocity ω. The force control device is connected to the bottom of the floating spindle, which drives the floating spindle to move axially and provides pressure P on the workpiece surface;  of the workpiece before and after the polishing test. The arithmetic mean deviation of the profile Ra represents the arithmetic mean of the absolute value of the profile deviation within the sampling length. It can reflect the material removal diverse within the measured microscopic area of the workpiece surface. Ten different points with equal intervals are selected around the middle area of the polished workpiece surface to measure Ra; and the experimentally measured standard deviation of material removal σ T is obtained according to Eq. (25).
where MRR T is the measured material removal rate test value; Δm is the change of mass before and after robotic wet polishing; S T is the wet polishing area; t is the polishing time; and ρ T is the workpiece density. The mean abrasive particle diameter d p , volume fraction f p , and polishing fluid flow velocity v q used in the above simulations are defined as a combination C i of the polishing fluid properties, leading to 64 combinations of the polishing fluid properties C 1 -C 64 in total. The robotic wet polishing experiments are carried out on 64 identical workpiece to obtain the experimentally measured material removal rates MRR T1 -MRR T64 and material removal standard deviation σ T1 -σ T64 . The experimental results (MRR T , σ T ) are compared with the theoretical results (MRR, σ) obtained from the proposed solid-liquid coupled twophase flow model. The error distribution between (MRR T , σ T ) and (MRR, σ) are analyzed, and illustrated in Figs. 16a and b. Each sphere in the figure represents an error value.
A larger radius of the sphere indicates a bigger error.
The results show that all the errors are less than 5%. The mean error of the material removal rate is 2.94%; and the mean error of the standard deviation of material removal is 2.27%. Similarly, the values of polishing rotation rate ω, robot swing speed v 2 , and polishing pad pressure P used in the above simulation are defined as a combination P i of the polishing kinematic parameters, leading to 64 sets of combinations P 1 -P 64 in total. The robotic wet polishing experiments are carried out as well on another 64 identical workpieces to obtain the experimentally measured material removal rate MRR T65 -MRR T128 and the material removal standard deviation σ T65 -σ T128 . The experimental results (MRR T , σ T ) are compared with the theoretical results (MRR, σ) obtained from the proposed modeling method. The errors between (MRR T , σ T ) and (MRR, σ) are shown in Fig. 17. From  Fig. 17a, it can be seen that the mean error between the experimental and simulated values of the material removal rate is 4.19%. But when the polishing rotation rate ω is high, there is a relatively large error between the experimental values and the simulation values of MRR in two tests, with an error of 6.7% and 8.1% respectively. In these two tests, the MRR T obtained experimentally are 45.39 nm/min and 46.15 nm/min respectively. They are lower than the MRR obtained in the simulation, which are 48.65 nm/min and 50.22 nm/min respectively. Figure 18 shows the surface state of the polishing pad and workpiece after robot wet polishing, wherein Fig. 18a shows the surface state of the polishing pad in high rotation rate wet polishing. It can be found that high rotation rate ω of the polishing pad causes a large centrifugal force on the polishing fluid and abrasive particles under the polishing pad area. Close to the center of the polishing pad, there is a low efficiency polishing area, in which the linear speeds of the abrasive particles are relatively low while the centrifugal forces on the abrasive particles are relatively high. When the rotation rate ω of the polishing pad rises, the radius of the low efficiency polishing area increases gradually. The abrasive grains in the polishing fluid under this area move away from the center of the polishing pad rapidly, causing a short effective polishing time in this region. As a result, the experimental value of the material removal rate MRR T is lower than the simulated value of MRR when the rotation rate ω is high. The results from Fig. 17b show that the experimentally obtained σ T and the simulated σ have a mean error of 3.77%. However, when P is high, there is a relatively large error between the experimental value σ T and the simulated value σ in two tests, with an error of 6.4% and 7.8%, respectively. The experimentally obtained σ T in these two tests are 52.41 nm and 63.35 nm, which are higher than the simulated values (49.26 nm and 58.77 nm). From Fig. 18b, it can be seen that when other process parameters are consistent and the polishing pressure P is high, there are obvious wet polishing marks locally on the surface of the workpiece, and higher roughness values Ra are measured for this region. Large scratches are normally caused by large particles in  Fig. 19. The diameter distribution of the abrasive particles in the polishing fluid is uniform. Although some of the abrasive particles agglomerate and become a larger size particle, they are separated with the effect of liquid turbulence or external force; and the possibility of generating large scratches on the workpiece surface is low. There are large particles on the polishing pad surface after wet polishing process. They are composed of abrasive particles, material powders removed from the workpiece surface, and scraps dropped from the polishing pad. These particles are embedded in the valley of the polishing pad surface. There is less turbulent effect on this region under the polishing pad area. With a large pressure P, the polishing pad is pressed to be more compact; and the particles get more close to the surface of the polishing pad. Therefore, it is easier to produce large scratches in the robotic wet polishing process.

Conclusions
In this paper, the distribution of the radial velocity and the thickness of the liquid film at different radii is clarified by exploring the flow characteristics of the polishing fluid. Through the construction of a solid-liquid coupled two-phase flow model for the wet polishing process, the interactions between the polishing pad wall, the workpiece wall, the abrasive particles, and the liquid are clarified; and the polishing fluid velocity, temperature change attributes, the transient position, and velocity distribution of abrasive particles under different typical working conditions are analyzed. Based on the material removal rate model considering