Numerical Simulation and Veriﬁcation of Laser-Polishing Free Surface of S136D Die Steel

As a novel polishing technology, polishing by laser beam radiation can be used to improve the sample surface ﬁnish without causing material losses. In order to study the effect of laser polishing on the surface morphology of S136D die steel, an L 16 (4 4 ) orthogonal experiment was designed to describe the variation trend of surface roughness with energy density. The two-dimensional transient model of laser polishing was established to simulate the evolution process of material surface morphology during laser polishing by combining numerical simulation with the experiment. The model uses a moving laser heat source to study the effects of capillary pressure and thermocapillary pressure in the laser polishing process. The experimental results show that the minimum roughness can be reduced to 0.764 µ m, and the error between the actual molten pool depth and the simulated molten pool depth is 5.3%.


Introduction
Polishing is a material removal process and is generally used to produce a smoother surface topography by mechanical, chemical, and electrochemical action [1][2][3]. Conventional polishing not only constitutes materials loss but also has the following limitations as well: low polishing rate, difficulty in automation, high consumption cost, waste treatment issue. Recent studies suggest that laser polishing (LP), as a novel technology for material surface polishing, has greater advantages than conventional polishing. The LP process can be proceeded without any material losses and available to achieve freeform surface polishing [3]. It is a surface finishing technique characterized by rapid heating and cooling in the process. During LP operations, the material surface can be heated to the melting point in a short period of time without material loss. Its efficiency is eight times bigger than that of the traditional polishing and 30 times bigger than that of the manual polishing of skilled workers under the same prerequisite [4][5][6]. As the heat is accumulated on the metallic surface to a certain amount a micro peak on the material surface will first reach the melting point and start melting, and the melted material flows from the peak to the valley due to the tension [7][8][9].
At the beginning, the problem of finite element simulation focused on the solid-liquid change, and how to distinguish and track the solid-liquid boundary was a major difficulty. T.A. Mai [10] used a self-developed finite difference heat conduction model to numerically simulate the effects of temperature changes, pool morphology changes, and potential phase changes on the liquid-solid boundary during rapid melting and solidification. The Stefan boundary equation was introduced into the model during the simulation process to track the actual position of the moving solid/liquid boundary in the mesh element where the phase transition occurred.
In the process of LP, the force driving the melt movement in the remelting region is mainly in the form of capillary pressure and thermocapillary pressure. D.A. Willis [11] calculated the transient velocity field, temperature field, and free surface motion by the numerical method. The possible mechanism of droplet formation was analyzed by using the experimental and numerical results. It is found that under the condition of low pulse energy, with the increase of pulse energy, the radial flow becomes fast enough, and the inertial effect dominates, which leads to the separation of droplets from the molten pool. In the subsequent study, it is concluded that the change of surface morphology and the formation of droplets are due to the flow caused by surface tension, and the recoil pressure effect does not play a significant role in the development of surface morphology. Pasandideh Fard [12] studied the phase transition, hydrodynamics, and heat transfer of molten metal dripping on a flat substrate by using a three-dimensional model, and they further studied the effect of surface tension on fluid motion. Hong Shen [13] proposed a forming mechanism of pulse laser melting stainless steel. A two-dimensional axisymmetric finite element model considering the relationship between surface tension and temperature was established by COMSOL 5.4a finite element analysis software. The melt flow and free surface deformation were analyzed by using normal force and shear force. It was proven that the thermal capillary force and capillary force had an effect on the melt flow and the corresponding free surface deformation. Chi Zhang [14] established a twodimensional axisymmetrical numerical model coupled with heat transfer and fluid flow through COMSOL, analyzed the role of capillary pressure and thermalcapillary pressure in the LP process, and simulated the specific process of molten pool flow. Shirzad Mohajerani [15] established a three-dimensional unsteady heat transfer CFD (Computational Fluid Dynamics) model in ANSYS Fluent, a finite element analysis software, to numerically analyze the continuous wave (CW) laser polishing of H13 tool steel. The temperature distribution of the whole workpiece was obtained by calculating the width and depth of the molten pool with the finite volume model. The dependence of heat capacity, thermal conductivity, and density on temperature was considered in the established model. The appropriate value of the heat absorption coefficient of the workpiece was determined by the calibration experiment.
Perry [16][17][18] proposed a critical frequency f cr , which describes the cut-off point of the surface spatial frequency content. The amplitude above the cut-off point is expected to decrease significantly. The research shows that the maximum melt duration t m-max can be used to determine the minimum critical frequency f cr . If it exceeds f cr , the amplitude of surface and its surface spatial frequency content should decrease significantly. Madhu Vadali [19] used the transient two-dimensional axisymmetric heat transfer model with temperature dependence to predict the depth and duration of the melt, and they deduced the minimum critical frequency by determining the maximum melt duration required by the critical frequency. Then, inverse Fourier transform was carried out on the predicted spatial frequency content of the polished surface to obtain the predicted surface height data by no more than 15% error compared with the actual value.
Chao Ma [20] established a two-dimensional axisymmetric transient model coupled with heat transfer and fluid flow by using the finite element method, which not only predicted the evolution of the free-deformed surface profile but also studied the influence of laser pulse duration on the molten pool flow. In the research process, it was found that longer pulses would produce more obvious fluid flow. Furthermore, the melting and re-solidification of a molten pool are simulated, which provides more real information for the transient evolution and surface morphology change of the molten pool.
As can be observed from previous studies, simulation results of laser polishing are mainly concentrated on stationary heat source, compared with the stationary light source, and the simulation results by moving light source are more accurate and representative. The purpose of this paper is to study the capillary force and thermocapillary force with time and the evolution of the free surface under the moving heat source by the method of combining numerical simulation and experiment. In order to verify the numerical model, the molten pool morphology was measured and compared with the simulated molten pool morphology.

Experiment Materials
The S136D die steel is used in the experiment. The basic chemical composition is shown in Table 1. The S136D steel has good corrosive and abrasive resistance and polishing performance, due to the high chromium and low sulfur. Prior to the experiment, the sample is quenched and tempered, cut into a 120 × 115 × 5mm 3 steel plate by wire, and cleaned before polishing. The initial surface roughness Ra = 5.358 µm.

Experiment Device and Testing Equipment
Two lasers were used in the experiment, which are continuous fiber laser. The power range is 150-1000 W, the focal length is 720 mm, the zoom range is 60 mm, and the laser spot diameter is 0.3 mm at the focal point. The experiment device and principle used in the experiment are shown in Figure 1a,b. The whole laser equipment is composed of laser emitter (MFSC-1000W, Chuangxin Laser, Shenzhen, China), Scan header, 2-Axis CNC rotation table, Process chamber, and other devices. Since the laser emitter is inside the whole equipment, it is not shown in Figure 1a. The laser beam was transmitted by the laser emitter. The diameter and divergence angle of the laser beam were adjusted by the beam expander. Inside the laser galvanometer, there were mirrors to adjust the laser beam position and to ensure the laser irradiation on the workpiece surface. The white light interferometer (BRUKER WYKO Contour GT-K, Bruker, Billerica, MA, USA) and metallographic microscope were used to observe the sample surface topography and measure the surface roughness Ra. The elemental content of the initial surface was measured by EDS (AZtec X-MaxN 80, Oxford Instruments, Abingdon, UK). The polished sample was cut into small cubes by wire. To avoid the heat influence on the test results during wire cutting, the roughest sandpaper was used to grind a 1 mm cross-section of the cube and polish it to a mirror finish. The molten pool shape of the polished cross-section was observed by SEM (ZEISS Gemini 300, Gina, Germany).

Experimental Methods
Prior to the polishing experiment, the initial surface roughness of the sample is 5.358 µm measured by a white light interferometer. The samples were first polished by continuous laser with different parameters in an area of 10 mm × 10 mm. After polishing, the effect of various laser factors on surface roughness was to find the best set of parameters. The laser beam in the experiment moved in the x-axial direction at first to cover the whole polishing area and began to move along the y-axis until the same area was covered again.
In order to explore the influence of energy density on surface roughness and spatial frequency, an L 16 (4 4 ) orthogonal experimental table is established, as shown in Table 2, in which P is laser power (W), V is scanning speed (mm/s), D is scanning distance (mm), FO is focus offset distance (mm), and ED is energy density (J/mm 2 ).  Figure 2 is a scatter plot of the relationship between surface roughness and energy density. The surface roughness (Ra) decreased first and then increased with the increase of ED. The polishing effect of ED is better when it is about 15-25 J/mm 2 . When ED is below 20 J/mm 2 , the relationship between ED and Ra is not randomly distributed. In this range, Ra decreases gradually with the increase of ED. There are two main reasons for the error. Firstly, it is due to the difference in the material surface microstructure. Neither of the polishing areas feature the identical microstructure. The reflection and the contact area are different with laser irradiation on the material surface.

Results and Discussion
Secondly, it is due to the material temperature. The material selected in the experiment is a S136D die steel plate of 120 × 115 × 5mm 3 , which means hundreds of areas can be polished. Therefore, the residual heat generated during the first polishing exists when polishing the second area. Furthermore, the energy density before and after polishing differs, leading to unidentical residual heat and initial temperature.
In the whole data, with the increase of ED, Ra first decreases and then increases, because at low energy density, the laser energy is not enough to melt the material surface completely, resulting in insufficient melting of the material surface and the initial surface cannot be smoothed. When the energy density is too high, the surface of the material melts too deep, which leads to the cooling of the molten material before the oscillation stops, and the surface of the material regenerates the uneven profile, forming the SOM phenomenon, so it can not achieve good polishing effect [21]. When the energy density is 18 J/mm 2 , the polishing effect is better, and the roughness can be reduced to 0.764 µm.
The measurement area of the white light interferometer is set at 2.5 mm × 1.9 mm; Figure 3a shows the unpolished three-dimensional topography image; Ra is 5.358 µm. The value of the peak Rp = 32.975 µm, the valley Rv = -26.151 µm, and the maximum distance between the peak and valley Rt = 59.127 µm. Figure 3a shows the three-dimensional image of the sample 11 surface topography. It can be noted that the surface roughness Ra = 0.764 µm, peak value Rp = 6.374 µm, the valley value Rv = −5.448 µm, and the maximum distance between the peak and valley Rt = 11.822 µm. Ra reduced by 85.74% compared with the initial surface. Spatial Fourier analysis transforms a spatial signal into a frequency signal to get shape, frequency, and amplitude. The spectrum analysis method is used to analyze the surface profile and power spectral density (PSD) curve before and after LP is obtained [22]. Figure 4a,c are the spatial frequency amplitudes of the initial surface in X and Y directions. The maximum frequency amplitudes in the X and Y directions are 3.38 and 9.6, respectively. Since the surface after WEDM will produce stripes, the amplitude parallel to the stripe spatial frequency (X direction) will be smaller than that perpendicular to the stripe spatial frequency (Y direction).  After laser polishing, the spatial frequency of samples 1-16 decreases more or less in both the X and Y directions, and it decreases first and then increases with the change of ED. When the ED is about 18 J/mm 2 , a better amplitude reduction effect can be obtained.

Model Assumptions
In order to simplify the calculation and reduce the calculation time, the following assumptions are made for the model: (1) The material is regarded as isotropic, and its thermophysical parameters are only related to temperature. (2) The material was tested in high-purity nitrogen without considering the process of surface oxidation.

Governing Equations
The transient temperature field is controlled by Fourier's law described by the following equation where ρ is the density, t is the time, k is the thermal conductivity, and u is the velocity field in the molten pool calculated according to the Navier-Stokes equation, as shown below: where p is the pressure, µ is the dynamic viscosity, F v is the volume force, including buoyancy force F b and gravity F g [13].
where ρ 0 is the reference density, g is the gravity constant, β is the thermal expansion coefficient, and T re f is the reference temperature. In the numerical, ρ 0 = ρ l and T re f = T l . For incompressible fluid, the continuity equation is simplified as: The continuous laser melting and remelting process of S136D die steel is studied in this paper. The model is solved in a single domain including solid and liquid phases. Considering the effect of latent heat in the melting process, the equivalent heat capacity C * p is used to replace the specific heat capacity where T is temperature, C p is specific heat, L m the is latent heat of melting, and f L is the liquid fraction, which is expressed as follows: where T s and T l are the solidus and liquidus temperatures of the material.

Moving Mesh
To simulate the deformation of the free surface, the Arbitrary Lagrangian-Eulerian (ALE) approach is adopted. In this method, the displacement of boundary nodes is decided by the fluid dynamics. While domain nodes are subject to an Eulerian description, the boundary nodes follow a Lagrangian description. This method can not only deal with the flow with large distortion, but it also accurately describe the internal motion of the fluid. Equation (7) describes the coupling of ALE with fluid flow equations. u muse ·n = u·n (7) where u muse is the mesh velocity (m · s −1 ), and u is the material velocity computed from the Navier-Stokes Equation (2).

Boundary Conditions
Different initial surfaces affect the polishing results. The definition of the initial surface in the simulation is very important. It is not accurate to define the initial surface of the model only by a simple Ra value. In this model, a white light interferometer is used to measure the initial surface profile of the sample, Fourier filtering is performed, and then the filtered curve is imported into COMSOL. A two-dimensional axis symmetry (cylindrical coordinate r, z) numerical model of surface melting that is transient under the influence of surface tension and the Marangoni effect is established by using a COMSOL multi-physical field. The geometric shape and boundary division of the calculation area are shown in Figure 5.

Heat transfer:
The labels of all the boundaries in the calculation domain are shown in Figure 2. Boundary 1 is the axis of symmetry, boundary 2 is under the boundary conditions of generalized heat flux and convective heat flux, and the surface affects both boundary 2 and boundary 3 on the environmental radiation, as described in Equation (8): where h is the natural convection coefficient, ε is the emissivity, σ is the Stefan-Boltzmann constant. T a is the ambient temperature, I is the flat top light source, It is expressed by Equation (9), x 1 is the length from the center of the spot, and f (x 1 ) is the range of laser action. When the laser diameter exceeds the range, the laser energy is regarded as 0, which is expressed by Equations (10) and (11): P is the laser power, r is the spot diameter, α is the absorptivity of the material, and v speed is the laser scanning speed. Boundary 4 is defined as adiabatic and is represented by Equation (12): Fluid flow: Boundary 2 is defined as a free-form deformable surface, on which the capillary pressure (surface tension) is a normal force, expressed by Equation (13), and the thermal capillary pressure (Maragoni effect) exerts a tangential force along the surface, as expressed by Equation (14). σ n = κγ · n (13) where κ is the curvature of the surface profile, γ is the surface tension, ∇ S T is the temperature gradient along the tangential direction of the surface, n and t are the normal and tangential unit vectors of local surfaces, respectively. The tangential force can be applied to boundary 2 by the inherent multi-physics interface of COMSOL, while the normal force is implemented as a weak contribution to the boundary. The detailed theoretical derivation of the weak form is shown as follows.
The first item S γ(∇ u) · dS can be written in COMSOL weak form format, and the detailed implementation is covered in [23]. The second item ∂S γ( u · n)dI, which is called a contour integral, is equal to zero.
The specific boundary conditions are shown in Table 3.

Properties and Parameters
The material used is S136D die steel. Since S136D is a new type of plastic mold steel, many thermophysical parameters can not be found, so most of the thermophysical parameters are quoted in the same 304 stainless steel [13,24,25], and the specific parameters are shown in Table 4. Most of the parameters are related to temperature. For the sake of accuracy, the transition from solidus to liquidus is regarded as a process within a certain temperature range. The melting temperature is assumed to be the average of the solidus and liquidus temperatures. In terms of fluid flow, the effective viscosity method is used to consider both the solid and liquid phases, which allows the liquid phase to flow but prevents any deformation of the solid part. The actual dynamic viscosity at the liquidus temperature is 5 mPa·s. Since there is no fluid flow, the viscosity of the solid phase is basically infinite. In this study, the viscosity of the solid phase was set at least 10 5 Pa·s, which is several times that of the liquid phase. In the melting zone, the dynamic viscosity decreased steadily from 10 5 Pa·s at T s to 5 mPa·s. Meanwhile, the constructed viscosity is a "smooth" step function, which not only satisfies the required difference between liquid and solid states but also achieves numerical convergence.

Mesh and Configurations
The direct solver used for simulation is PARDISO [26]. In the process of molten pool simulation, the influence of high-quality mesh on calculation results is very important. The triangular mesh is easier to converge, so triangular mesh is adopted in the whole model. Theoretically, the smaller the mesh size, the longer the calculation time of fusion, and the higher the accuracy of numerical simulation results. Since boundary 2 is a free-form deformation surface, when the material reaches the melting temperature, the boundary will deform, so it needs to use finer mesh for precise calculation. The maximum mesh element near boundary 2 is 1 µm, and the minimum mesh element is 0.02 µm. The final number of cells is 57,231, and the number of grid vertices is 29,300. The specific mesh division is shown in Figure 6 and Table 5.

Simulation Results and Discussion
In order to study the evolution of molten pool on the material surface caused by laser, the laser parameters used in the simulation are 180 W, 50 mm/s (parameter 11 in orthogonal experiment), and the spot diameter is 0.3 mm without considering the effect of defocus.

Evolution of Molten Pool Surface Profile
When the laser acts on the surface of the material, the material begins to absorb heat and the temperature rises. When the heating time reaches 1.3 ms, the material begins to melt and the temperature reaches 1550 K. As time goes on, the temperature increases, and the molten pool increases. Figure 7a is the temperature curve at 3, 6, 9, 12, and 15 ms. The black dotted line is the melting temperature of the material. Since the laser heat source of the model is at the speed of 50 mm/s, only the region of about 150 µm reaches the melting point above 3 ms, and the farther away from the laser heat source, the lower the temperature. With the moving of the laser heat source, the area above the melting point temperature increases gradually, and the maximum temperature moves with the laser heat source. The physical quantity used to characterize the top surface of the molten pool is spatial curvature, which is defined as follows: κ = − r ∂ 2 ϕ/∂r 2 + (∂ϕ/∂r) 1 + ϕ 2 r r + 1 + (∂ϕ/∂r) 2 3/2 (16) where ∂ϕ/∂r is the vertical displacement of the top surface. In order to better explain the polishing process, Figure 7b reflects the curvature curves at 0, 3, 6, 9, 12, and 15 ms. the larger the spatial curvature, the greater the undulation of the surface profile. At 3 and 6 ms, the curvature in the melting zone is smaller than the initial surface profile in the same position under the action of surface tension, gravity, and other forces. At 9, 12, and 15 ms, the space curvature of the molten pool is lower than the initial space curvature. With the movement of the laser heat source, the space curvature of the front molten material is higher than that of the melting material after re-solidification. This is because the force acting on the surface of the material begins to solidify before equilibrium during cooling, which leads to the formation of an uneven surface, The spatial curvature of the polished surface is much lower than that of the initial surface, which indicates that the polished surface is smoother. Due to the different peak temperature of the heating center, the surface tension distribution profile is also different. Figure 8 shows the surface profile change and the flow velocity of the molten material during the laser polishing process; Figure 8a-e represent the surface profile and flow state of molten material at 3, 6, 9, 12, and 15 ms, respectively. At 3 ms, the molten pool is initially formed, and the flow velocity of the molten material is small. Relatively speaking, due to the influence of gravity and surface tension, the velocity is larger at the place with larger curvature, and the left side of the molten pool has been circulating in this time period. Since the molten pool is small and the flow velocity is small, the effect of gravity and buoyancy is not obvious.
At 6 ms, the molten pool has begun to take shape, and the molten material flow velocity increases from 20 × 10 −3 m/s at 3 ms to 0.12 m/s. Due to the temperature difference between the center and the edge of the laser heat source, the temperature at the center is higher than that at the edge, and the surface tension will decrease with the increase of temperature, forming a temperature gradient. Therefore, in the range of the material covered by the laser heat source, there is a gradually increasing surface tension from the center to the edge, forming a surface tension gradient, and the molten material will flow from the position with low surface tension to the position with high surface tension, Furthermore, the surface tension at the bottom is higher than that at the top, so the molten material flows from the right side along the bottom to the left. Due to the buoyancy, the molten material flows on the net, and then due to the lowest temperature on the left side and the maximum surface tension, the molten material continues to flow to the left side until the solid-liquid boundary; then, it flows to the bottom of the molten pool under the influence of gravity and then flows to the right side, thus forming convection.
In the range of laser action, the initial surface is gradually smooth. At 6 ms, it can be observed that the right contour gradually melts. However, when the molten pool reaches 9 ms, it can be observed that the previous contour is not completely flat. This is because at the beginning of laser operation, the molten pool is small, the action time is short, and the molten pool maintenance time is not enough to make the surface completely smooth.
At 9, 12, and 15 ms, the molten pool is stable, and the flow velocity of the molten material increases gradually. It can be seen from Figure 8c that the left and right sides form a reflux, respectively, and the center area of the laser heat source is in the middle. The temperature is the highest and the surface tension is the lowest, so the molten material will flow to both sides; then, due to the influence of gravity, it will flow downward, then to the center, and finally, due to the buoyancy, it will flow upward to form convection. Compared with the initial surface, the material surface after laser polishing is much smoother.
A comparison between the values of capillary pressure |κγ| and thermocapillary pressure |(∂γ/∂T) · (∂T/∂r)| along the r direction reveals where the capillary and thermocapillary forces dominate, respectively, at a specific time and development of capillary and thermocapillary regimes, as shown in Figure 9. The graph reflects the changes of capillary force and thermocapillary force at 3, 6, 9, 12, and 15 ms. At 3 ms, when the laser starts to run, the capillary force in the molten pool is dominant, so the capillary force first determines the deformation of the free surface and effectively smoothes a large number of surfaces. Under the influence of temperature gradient, the Marangoni effect has been produced, but the influence of the thermal capillary force on the surface of the molten pool is less than that of the capillary force. With the increase of temperature, the melting zone expands, the dynamic viscosity decreases, and the temperature gradient at the top of the molten pool becomes obvious. At 6 ms, the Marangoni effect at the edge of the molten pool is more significant, so the capillary force is dominant in the center of the molten pool, and the thermal capillary force is dominant at the edge of the molten pool. However, compared with the right side, the thermal capillary force on the left side is more obvious, because the curvature on the right side is larger, which leads to larger surface tension. Moreover, the laser runs from left to right, and the molten pool on the right side is smaller than that on the left side, which makes the Marangoni effect less obvious than that on the right side; that is, the thermal surface tension is less obvious than that on the right side. At 9 ms, the capillary force is the main force, at 12 ms and 15 ms, the thermal capillary force is the main force on both sides, and the capillary force is the main force in the middle.

Experimental Validation
At 15 ms, the molten pool is completely stable, so this position is selected as the final shape of the molten pool, as shown in Figure 10. The simulated molten pool width is 353.3 µm (Figure 10a), and the actual molten pool width is 408.6 µm (Figure 10b). Compared with the actual size, the simulated error is 13.5%. The laser spot diameter is 300 µm. When the laser acts on the surface of the material, there will be energy transfer, so the actual molten pool width will be greater than 300 µm. The simulated molten pool depth is 67.5 µm, the actual molten pool depth is 71.3 µm, and the error is 5.3%. Since the actual molten pool is scanned by a laser line on the sample surface, its cross-section is obtained perpendicular to the scanning direction. The simulated molten pool is parallel to the laser path. Therefore, the simulated molten pool width has a larger error than the actual molten pool width.

Conclusions
The effect of laser energy density on the surface morphology was investigated. When the ED was 18 J/mm 2 , the polishing effect was the best, and the Ra decreased from 5.358 to 0.764 µm. When the ED is less than 15 J/mm 2 , the energy absorbed by the material is not enough to completely melt the rough surface, resulting in poor polishing effect. When the energy density is greater than 25 J/mm 2 , the forces of gravity, buoyancy, and surface tension inside the molten pool have not reached equilibrium, and the molten pool has cooled and solidified, resulting in the production of uneven surfaces.
The influence mechanism of capillary pressure and thermocapillary pressure on the surface morphology evolution of a molten pool is simulated. The simulation results show that when the surface curvature of the molten pool is large, the capillary pressure plays a leading role, and the material at the top of the molten pool flows into the groove. After 6 ms, the molten pool gradually stabilized, and the thermocapillary pressure caused by the temperature difference played a leading role at the edge of the molten pool, while the capillary pressure still played a leading role in the center of the molten pool.
The experiment results show that the actual molten pool width is 408.6 µm. The simulated molten pool width is 353.3 µm. Compared with the actual size, the simulated error is 13.5%. The depth of molten pool is 71.3 µm in the actual situation, and the simulated depth is 67.5 µm, with an error of 5.3%.