Towards the Determination of Machining Allowances and Surface Roughness of 3D-Printed Parts Subjected to Abrasive Flow Machining

Abrasive flow machining (AFM) is considered as one of the best-suited techniques for surface finishing of laser powder bed fused (LPBF) parts. In order to determine the AFM-related allowances to be applied during the design of LPBF parts, a numerical tool allowing to predict the material removal and the surface roughness of these parts as a function of the AFM conditions is developed. This numerical tool is based on the use of a simplified viscoelastic non-Newtonian medium flow model and calibrated using specially designed artifacts containing four planar surfaces with different surface roughnesses to account for the build orientation dependence of the surface finish of LPBF parts. The model calibration allows the determination of the abrasive medium-polished part slip coefficient, the fluid relaxation time and the abrading (Preston) coefficient, as well as of the surface roughness evolution as a function of the material removal. For model validation, LPBF parts printed from the same material as the calibration artifacts, but having a relatively complex tubular geometry, were polished using the same abrasive medium. The average discrepancy between the calculated and experimental material removal and surface roughness values did not exceed 25%, which is deemed acceptable for real-case applications. A practical application of the numerical tool developed was demonstrated using the predicted AFM allowances for the generation of a compensated computer-aided design (CAD) model of the part to be printed.


Introduction
In recent years, additive manufacturing (AM), or 3D printing, became a widely used and rapidly developing technology that simplifies the production of complex parts with reduced weight and improved functionality. At the same time, laser powder bed fusion (LPBF) became one of the most mature technologies for the 3D printing of metal parts. Despite the advantages presented by LPBF, an excessive as-built surface roughness inherent in the process impacts the fluid flow and the heat transfer characteristics of LPBF parts, and alters their mechanical properties. Strano et al. [1] studied the impact of the LPBF build orientation on the surface topography, while Kruth et al. [2] provided a detailed explanation of mechanical limitations associated with the surface finish of LPBF parts. They determined that the bigger the angle between the surface of a printed part and the building platform, the higher the surface roughness. For example, as reported by Urlea and Brailovski [3], if the build angle of a Ti-6Al-4V component increases from 0 to 135 • , the as-printed surface roughness Ra increases from 4 to 23 µm.
Different approaches and techniques were developed to improve the surface roughness of LPBF parts. For example, Yadroitsev and Smurov [4] improved the surface finish by optimizing the laser scanning strategy. From their side, Cherry et al. [5] considered the improvement of the as-built surface roughness by varying the scanning speed. However, for most applications, the as-built surface roughness still requires additional surface finishing. In this context, Alrbaey et al. [6] proposed the use of the surface laser meltingremelting strategy, but this approach is limited to upper surfaces of printed parts. Urlea and Brailovski [7] worked on electrochemical polishing of LPBF parts, but pointed out the difficulty of electropolishing parts with intricate internal geometries, such as internal cavities and channels. To overcome the latter difficulty, chemical polishing appears to be a promising alternative, since it facilitates access to difficult-to-reach surfaces. However, as shown by Pyka et al. [8], chemical polishing alone cannot entirely eliminate the LPBF surface imperfections. To improve the preceding technique, Mohammadian et al. [9] combined the chemical and the abrasive flow polishing techniques to polish 3D printed parts but pointed out the risks of surface contamination, as well as the health and environmental hazards related to the use of chemicals.
Considering the above, forced flow of a chemically passive abrasive medium across the surface of a part to be polished can be considered as an effective and safe surface finishing technology for LPBF parts and particularly those with complex internal channels. As discussed by Rhoades [10], this process called abrasive flow machining (AFM) has been used for decades to deburr, polish or radius surfaces and edges of parts with complex outer and internal shapes. For example, Cheng, et al. [11] studied AFM of gas turbine engine blades, while Ferchow et al. [12] used this technology to polish curved channels of 3D-printed nozzles.
It should be noted that abrasive media used in AFM mainly consist of a liquid polymer carrier and dispersed abrasive particles. Kumar and Hiremath [13] reviewed the AFM media applications and found that the mostly used polymer carriers are polyborosyloxane and silicone rubber, while for abrasives, they are silicon and boron carbides, aluminum oxide, and polycrystalline diamond particles. A typical abrasive particle size lies in the range of 32-1035 µm, as mentioned by Trengove [14], while concentration of the abrasive may approach 80%, after which the medium begins to act inefficiently, as demonstrated by Kar et al. [15]. An interesting method to assess the cutting forces of abrasive flow machining using a two-phase viscoelastic flow approach was proposed by Dong et al. [16] This method combines an analytical model for calculating the cutting factors and an experimental technique to calibrate this model. While this approach provides an interesting insight into the abrading behavior of AFM media and offers a new perspective for the optimization of the media composition, it cannot directly be applied to predict the material removal (MR) during the AFM process and, therefore, to determine the machining allowances of the parts to be polished.
As indicated by Bouland et al. [17], this last endeavor represents the main challenge in the use of AFM for the finishing of 3D-printed parts. Thus, a numerical tool allowing the MR and the surface roughness prediction, depending on the AFM process parameters and the geometric and material attributes of the part to be polished, is required. Different computer fluid dynamic-based AFM MR models were reviewed when preparing this paper. Jain et al. [18] proposed to calculate MR using an analytical model that considers the material hardness, the abrasive grain size and concentration, the normal stresses and the velocity of a Newtonian (linear viscosity) fluid. Wan et al. [19] extended this formulation by adopting a non-Newtonian fluid (power law) model and introducing the MR dependence on the initial and limiting surface roughnesses. Cheng, Shao, Bodenhorst and Jadva [11] also used the non-Newtonian fluid model, but in this case, they reduced the MR model formulation by using experimentally found constants dependent on the media type and the workpiece material and a variable flow shear rate.
Upon close examination, all these models represent the application of a classical theory of metal-cutting processes, with the main challenge consisting in an adequate computational fluid dynamics (CFD) simulation of the AFM media flow across the surface of a part to be polished, and allowing the calculations of the velocity and force fields acting on this surface. However, in all these models, the AFM medium is simulated as a purely viscous fluid, thus neglecting the fact that it should rather be considered as a viscoelastic fluid, as indicated by Uhlmann et al. [20]. However, no publications that account for the viscoelastic behavior of AFM media were found at the time of writing. Moreover, as can be found in Dash and Maity [21], no-slip boundary conditions are generally applied on walls during the CFD analysis of the AFM process, which contradicts the assessment of Uhlmann et al. [22] that slip cannot be neglected while simulating flow of viscoelastic fluids.
Based on the above considerations, the ultimate goal of this study is to develop a numerical model allowing the prediction of AFM allowances and the generation of compensated 3D CAD models of LPBF parts as a function of the selected AFM media, process parameters and required surface finish. From the theory of metal-cutting processes described in [23], the volume of the material removed MR v (m 3 ) may be calculated as: where P c is the cutting force, N; E 1 is the specific energy required for cutting, J/m 3 ; K 1 is the material removal factor (as shown by Brinksmeier et al. [24], K 1 corresponds to the Preston coefficient from the theory of abrasive polishing), m 3 /J; l is the distance over which the cutting force acts, m; v is the cutting velocity, m/s; and t is the time, s. By analogy with Equation (1), the thickness of the material removed during AFM, MR t (m), could be expressed as a function of the cutting stress N (Pa) and the abrasive medium velocity v (m/s) at each segment of the polished body, and the time of polishing, t (s): where A c is the surface of the cutting segment, m 2 . The abrading (cutting) coefficient K 1 , as well as the cutting stress N and velocity v fields, depend on the following characteristics:

•
Polished part hardness and geometric attributes (shape and surface roughness); • Rheological properties of the abrasive medium that depend on its composition (type of the viscous carrier + type and concentration of abrasive particles) and the temperature of use; • AFM operation conditions: inlet and back pressures.

CFD Simulations: Simplified Viscoelastic Model
For a given "abrasive medium/polished part" combination, the MR t prediction requires the knowledge of the abrading coefficient K 1 and the v and N fields generated by a flow of the abrasive medium across the surface of the part to be polished. To calculate the unknown v and N fields, ANSYS Polyflow software was selected. Three CFD models proposed by ANSYS Polyflow were evaluated: the generalized Newtonian model, the differential viscoelastic model and the simplified viscoelastic model. The first two models were excluded from consideration because the generalized Newtonian model could not predict normal stresses and the differential viscoelastic model was too computationally expensive. Moreover, as described by Walters and Webster [25] and Niethammer et al. [26], the latter presents numerical difficulties for CFD modeling with high Weisenberg numbers (We > 1), and under some flow conditions of this study, We can indeed reach a value of 140. Based on the above considerations, the simplified viscoelastic ANSYS Polyflow model was ultimately selected to solve the momentum (3) and the incompressibility (4) equations: where p is the pressure, Pa; T is the total extra-stress tensor, Pa; f is the volume force, N; ρ is the density, kg/m 3 ; a is the acceleration, m/s 2 ; and v is the velocity, m/s. The shear and normal stresses acting on the viscoelastic element are shown in Figure 1: Macosko [28] describes that the first normal stress difference (N 1 ) is necessary to keep the block at a constant thickness x 2 , at this T 11 − T 22 ≥ 0: As shown in Figure 1, three normal stresses act on the element (T 11 , T 22 , T 33 ) and, therefore, the second normal stress difference N 2 is: For an ideal rubber N 2 = 0, whereas for real polymers, N 2 is typically much smaller as compared to N 1 . As a result, for a simple shear flow, the simplified viscoelastic model incorporates N 1 into the force balance and takes T 22 = T 33 = 0, such that the total extrastress tensor T is given by: γ is the shear stress component, Pa; ψ·µ χ is the first normal stress, Pa; η . γ is the shear-rate dependent viscosity, Pa·s; µ . χ is the normal viscosity, Pa·s; . γ is the shear rate, 1/s; . χ is the viscoelastic variable, 1/s; and ψ is the weighing factor. To express both the shear rate-dependent η . γ and normal µ . χ viscosities of an AFM medium, the Carreau-Yasuda model of a non-Newtonian fluid was adopted: where η 0 is the fluid viscosity at zero shear rate, Pa·s; η ∞ is the fluid viscosity at infinite shear rate, Pa·s; λ is natural time, s; and n and a are the power-law indexes. As a result, the first normal stress difference (N 1 ) becomes: where the viscoelastic variable . χ obeys the transport equation involving the characteristic (relaxation) time τ . γ , which is given by: The Equation (11) is such that the solution γ is recovered in simple shear flow. γ of the medium were carried out using a HR-2 discovery hybrid rheometer (TA Instruments, New Castle, DE, USA), with a parallel plate-plate setup (d = 25 mm; t = 1 mm). To control the slip conditions, a sandpaper (grit 120) was attached to the plates, as recommended by the rheometer manufacturer in the training manual [29]. The following protocol was used: (a) Conditioning: 25 • C soak time for 10 s followed by pre-shear performed at a shear stress of 100 Pa for 60 s, followed by final equilibration for 300 s; (b) Frequency oscillation: 25 • C soak time for 60 s followed by testing under a standard controlled stress of 800 Pa with a logarithmic sweep (angular frequency: 628.3-1 rad/s; 10 points per decade); (c) Data acquisition: conditioning and sampling times 3 s; 64 points in waveform; 0.0 µN·m controlled flow torque, and (d) Conditioning at end of test: 25 • C temperature re-set.
Two types of 3D-printed parts were designed to support the model development: (a) V-shape artifact for model calibration ( Figure 2a); S-shape specimen for model validation ( Figure 2b). The calibration artifact was designed to reflect the dependence of the surface roughness of LPBF parts on their build orientation. To this end, the V-shape artifact contains four distinct zones, A, B, C and D, which correspond to the build orientation angles of 90 • , 45 • , 90 • and 135 • , respectively. At the same time, the V-shape artifact possesses a simple flat-plane geometry facilitating measurements of its surface roughness and dimensions. The S-shape specimen was designed to represent one of the most typical 3D parts, namely, tubular channel components for conformal cooling applications. To facilitate inspection, the S-shape specimen ( Figure 2b) is divided into two halves (top and bottom) to be assembled before polishing. To obtain identical initial surface roughness distributions (Ra 0 ) in the top and bottom halves of the specimen before polishing, their semi-channels were oriented upward during printing, thus providing a symmetrical build angle variation of 0 • to 90 • in both of them.
All parts in this study were manufactured using an EOSINT M280 400W Ytterbium fiber laser powder bed fusion system (EOS GmbH, Munich, Germany), an EOS Stain-lessSteelCX powder, and a "Customized CX" parameter set with a 30 µm layer thickness. After printing, based on the EOS [30] recommendations, the parts were annealed at 850 • C for 30 min under argon atmosphere to relieve residual stresses caused by the process, and cut from the building plate. The V-shape artifacts were divided into two groups: the first group was kept in the as-built state and all A, B, C and D surfaces of the second group were mechanically pre-polished down to Ra = 1.0 µm before AFM. The reason for creating the two groups of V-shape artifacts was the need to separate the influence of the overall part geometry on the material removal during the AFM process from that of its surface roughness. As it will be demonstrated herein, the pre-polished artifacts were used for the calibration of the material removal (MR t ) model, whereas the as-built artifacts were used for the extraction of the Ra(MR t ) evolution function.

AFM Setups
Two AFM setups were used in this study, one dedicated to the AFM processing of V-Shape artifacts (calibration setup, Figure 3a), and the second, to the AFM processing of S-shape specimens (validation setup, Figure 3b). Note that the first setup was built as a laboratory testing rig, while the second setup represents a commercially available industrial-size piece of equipment. The AFM calibration setup ( Figure 2a) is shown in Figure 3a and described in details in [17]. Briefly, under the action of the two moving pistons (1), an AFM medium flows reciprocally from the two feeding chambers (2) into the fixture (3) holding the V-shape artifact (4). The V-Shape pre-polished and as-built artifacts were polished under the same AFM conditions (flow rate 31.18 × 10 −5 m 3 /s; duration up to 300 cycles). The polishing process was discretized as follows: 302, 454, 756, 1512, 3024, 6047, 7559 and 9071 s (10, 15, 25, 50, 100, 200, 250 and 300 sequential passes). After each polishing sequence, the V-shape artifact was removed and evaluated (weighing, MR t profile and Ra measurements). During the V-shape calibration trials, RAW data (p inlet ) were captured with 1 s intervals using a custom LabView-based system. Note that since the pistons (1) (Figure 3a) move back and forth under the action of a single actuator (moving assembly), this setup does not allow for the controlled application of back pressure during polishing.
The AFM validation setup (Figure 2b) is shown in Figure 3b and reproduces an industrial piece of AFM equipment, Vector 150 (Extrude Hone LLC, Irwin, PA, USA). In this case, the force was applied only by one piston at a time (p inlet = 1200 psi), such that the second piston moved under the pressure exerted by the polishing medium. Although no specific back pressure was applied to the second piston, an effective positive p back existed due to an unspecified AFM system resistance. The polishing process was discretized as follows: 2966, 2211, 1787, 5272, 6398, 11,287, and 10,760 s (1, 1, 1, 4, 5, 7 and 7 sequential passes). After each polishing sequence, the S-shape validation specimen was removed and evaluated (weighing and Ra measurements). During the S-shape trials, RAW data (p inlet , p back , flow rate Q) were captured with 0.5s intervals using pressure sensors (3100 30CPS (Gems, Plainville, CT, USA), 0-3000 psi) and linear encoders (flow rate measurements) connected to a dedicated data acquisition system.

Measuring Equipment and Protocols
Before and after polishing, the V-and S-shape parts were weighed using electronic scales SECURA324-1S (Sartorius, Goettingen, Germany) and AP-210 Analytical Plus (Ohaus, Parsippany, NJ, USA), both 0.1 mg. Next, the Ra and MR measurements on the V-shape artifacts were carried out after each polishing series using a VR-5200 3D optical microscope (Keyence, Osaka, Japan) (12× magnification, high resolution, fine mode, ±2.5 µm) and Keyence VR-5000 series software. On the other hand, the Ra measurements on the S-shape specimens after each polishing series were carried out using a Surfcom 1500SD3 surface roughness tester (Accretech, Tokyo, Japan) with an DM43825 stylus (Accretech, Tokyo, Japan) (tip radius 2 µm), while the 3D optical microscope MR and Ra measurements on these specimens were carried out for their initial and final states only.
MR measurements, V-shape: As an example, Figure 4 illustrates the MR measurements realized on a 90 • -oriented surface A of the V-shape artifact (Figure 4a). The same routine was repeated for the resting 45 • -, 90 • -and 135 • -oriented surfaces. For these measurements, the initial and polished (final in this case) states were first compared using a reference plane passing by two zones located at the left and right edges of the artifact, which were not subjected to polishing (Figure 4a). A "volume and area" module of the 3D microscope software was then used (see Figure 4b) to determine the void (V void , mm 3 ) and solid (V s , mm 3 ) volumes, the first corresponding to the void located below the reference plane, while the second corresponded to the material located above the reference plane. For a given number of passes (n), the corresponding material removed volume MR v,n was determined according to Equation (12): To determine the MR t value at each point of the surface, the corresponding initial and final AFM point clouds were exported from the 3D microscope into the CATIA software environment and aligned together using a 3D CAD model of the part. The MR t calculations were performed in the MATLAB environment using a customized point-by-point calculation routine based on the "findNearestNeighbors" function.
MR measurements, S-shape: For the MR t measurements on the S-shape specimens, the corresponding initial and final (AFM) point clouds were divided on the "base"/"channel" regions ( Figure 5a) and exported from the 3D microscope software. Because of the complexity of the profile measured, the alignment routine was realized using the "Align" module in the MeshLab software environment. The initial state was selected as the reference and the base regions were used for alignment. Then, the matched channel STL point clouds ( Figure 5b) were exported for point-by-point MR t calculations in the MATLAB environment. For the model validation, the experimentally found V-and S-shape MR t fields were assigned to the corresponding CFD mesh vertices.  Ra measurements, V-shape: To establish the Ra(MR t ) function, the same reference planes as for the MR v measurement were used (Figure 4a). However, for this study, 8 × 8 mm segments with a relatively uniform MR t evolution on each of the surfaces (MR t,n ) were selected, as shown in Figure 6a. Next, the "Multiline roughness" module of the 3D microscope with 30 equally spaced (~250 µm) lines was used to calculate the average Ra (Figure 6b). The MR t,n was found as defined by Equation (13): where A segm = 64 mm 2 : area of the segment studied. Ra measurements, S-shape: In addition to the Ra measurements for building angles 0 • and 40 • realized after each polishing series using a surface roughness tester (see Figure 7, region of interest ROI), these measurements were repeated using a 3D microscope as shown in Figure 8a (straight channel ROI), by positioning the specimens at an angle of 45 • . The "Cylinder" surface shape correction was then applied in the 3D microscope software to determine the Ra corresponding to a build angle (α) varying from 0 to 90 • (Figure 8b).

CFD Simulation Setups
The V-and S-shape CFD simulation setups are presented in Figure 9a,b and the corresponding CFD parameters are presented in Table 1. Each CFD setup represents a symmetrical system divided into two volumes: polishing volume and global volume, with the former corresponding to the zone containing a part subjected to polishing, either V-or Sshape, and the latter corresponding to the remaining zone of an entire AFM setup. For each system, the boundary conditions (BC) on walls were also divided into two corresponding BC setup conditions, with and without a part to be polished. For the CFD analysis, the media rheology [η . γ , τ], the inlet flow rate (Q), the back pressure (p back ) and the wall slip coefficient (k slip ), characterizing friction between the abrasive medium and the polished part, were used as input parameters, whereas the simulated inlet pressure (p inlet ) was used as a control parameter.

CFD Simulation Setups
The V-and S-shape CFD simulation setups are presented in Figure 9a,b and the corresponding CFD parameters are presented in Table 1. Each CFD setup represents a symmetrical system divided into two volumes: polishing volume and global volume, with the former corresponding to the zone containing a part subjected to polishing, either V-or Sshape, and the latter corresponding to the remaining zone of an entire AFM setup. For each system, the boundary conditions (BC) on walls were also divided into two corresponding BC setup conditions, with and without a part to be polished. For the CFD analysis, the media rheology [ , ], the inlet flow rate ( ), the back pressure ( ) and the wall slip coefficient ( ), characterizing friction between the abrasive medium and the polished part, were used as input parameters, whereas the simulated inlet pressure ( ) was used as a control parameter.  During the CFD analysis, simulations were performed by swapping the inlet/outlet boundary conditions (BC). The justification for this swapping was two-fold: (a) each of the modeling cases represents the actual back and forth movements during the AFM process and (b) the simulated normal stress fields (N 1 ) of the simplified viscoelastic model are sensitive to the surfaces, either inlet or outlet, to which the BC are applied. Next, the average (back + forth)/2 solutions were obtained. To reduce the computational cost, "symmetry" BCs (boundary conditions) were applied during the simulations of 3D volumes, which resulted in a 50% reduction in the actual (Q real ) inlet flow rate (Q CFD = Q real /2). The following assumptions were made to facilitate the simulations: γ was taken at the cross-over of the G and G graphs and then calibrated using the MR results obtained with the pre-polished V-shape artifacts.  Figure A1.) Outlet p back = 0 psi Wall Generalized Navier's slip: f s = −k slip ·v s where, f s : shear force k slip : material parameter ** v s : tangential velocity at wall * Simulation flow rate Q CFD was reduced in half compared to the actual value due to the applied "symmetry" BC. ** τ, k slip : independent variables.

Model Calibration and Validation
Following the CFD sensitivity study on V-shape artifacts, the wall slip coefficient (k slip ) and the material relaxation time (τ) were found to have the most impact on the CFD simulation results (v, N 1 ) with k slip influencing primarily the p inlet and the v and N 1 fields and τ influencing primarily the N 1 field. Therefore, for the MR model calibration, k slip and τ were chosen as independent variables, and the MR t solution was established in the form: The removed material thickness (MR t ) and time (t) of the AFM process are considered as control parameters, since they could be determined directly from experiments. In contrast, K 1 as well as k slip and τ are unknown parameters requiring calibration. Thus, the MR model development is divided into two phases (see Figure 10):

1.
The model calibration phase using two types of V-shape artifacts: (a) The pre-polished V-shape artifact serves for the calculations of the stress and velocity fields of the polishing medium, and it is realized in two steps: -First step: k slip calibration. The entire AFM setup containing the V-shape artifact ( Figure 9) is considered to be fully polished and having an identical and constant k slip = k pol coefficient, irrespective of the polishing stage. The k slip value is adjusted to equalize the simulated and the experimentallymeasured inlet pressures: p inlet,CFD = p inlet,exp , and then applied as the BC on walls. -Second step: τ and K 1 calibration. The relaxation time (τ) and the material abrading coefficient (K 1 ) are adjusted to reach the best fit between the numerically predicted (MR t,mod ) and the experimentally measured (MR t,exp ) material removal values. To assess the degree of fitness, the maximum coefficient of determination (R 2 ) corresponding to the proportion of variance between the dependent and independent variables is found according to Chicco et al. [31] as: (b) The as-built V-shape artifact is used to establish the dependence of the experimentally measured surface roughness on the numerically calculated material removal, Ra(MR t , Ra 0 ), where Ra 0 is the initial (as-built) wall roughness. 2.
The model validation phase using S-shape specimens: At this stage, the k slip and τ values and the calibrated MR rate (N 1 , v) and Ra(MR t , Ra 0 ) models are used to calculate the material removal and the surface roughness at each point of the S-shape specimens during their polishing, and the results obtained are compared with their experimentally obtained equivalents to conclude on the validity of the proposed modeling approach.

Rheology of the LV-60B Abrasive Medium
The results of the complex viscosity (η * ), storage modulus (G ) and loss modulus (G ) measurements are given in Figure 11. As described by Cox and Merz [32], the shear rate correlates with the angular frequency ( . γ ∝ ω) whereas the shear viscosity correlates with the complex viscosity (η ∝ η * ). Analysis of the η * (ω) dependence demonstrates a shear-thinning behavior of the LV-60B medium (when the shear rate increases, the complex viscosity decreases). The best-fit correlation with the Carreau-Yasuda model η . γ (8) was obtained using the ANSYS [33] Polymat module and is presented in Table 2. Another important parameter characterizing viscoelastic fluids is the relaxation time (τ), which corresponds, at a first approximation, to the G and G moduli cross-over (see Table 2):

V-Shape: Weight and Material Removal Evolutions
The results of the MR weight measurements are presented in Figure 12. It can be noted that for both the pre-polished and as-built artifacts, the MR(t) functions are almost linear, and that the amounts of the material removed are close, while being slightly less significant for the pre-polished artifact. It can also be noted that at the beginning of the process, the material removal rates are unstable, especially for the as-built artifact, and that they then decrease, stabilize and become relatively constant for longer AFM sequences. The surface observations and the MR(t) field evolution during polishing of the pre-polished and as-built artifacts are given in Figure 13, while average numerical values of the removed material on surfaces A, B, C and D are plotted in Figure 14.

Calibration of the CFD Model Using Pre-Polished V-Shape Artifacts First
Step: Calibration of the k pol Value During the AFM experiment, inlet pressure measurements were constantly taken and an average value of p inlet,AFM = 1.46 × 10 6 Pa was calculated. Next, it was assumed that in the polished state, both the AFM setup and the pre-polished V-shape artifact possess the same k slip . Thus, a series of CFD simulations with k slip varying from 0.5 × 10 5 to 5.0 × 10 5 [N·s/m] on the V-shape artifact were executed. The inlet pressures of the CFD solutions (p inlet,CFD ) were then compared with the inlet pressure measured during the AFM experiments (p inlet,AFM ), and k pol = 2 × 10 5 [N·s/m] was found to be the one satisfying the p inlet,CFD = p inlet,AFM = 1.46 × 10 6 [Pa] condition. Second Step: Calibration of the τ and K 1 Values From Figure 14a, we concluded that the MR evolves linearly as a function of time. Therefore, for the pre-polished artifact, the experimental material removal rate, MR t,rate−exp , must be equal to the numerical material removal rate, MR t,rate−mod (v, N 1 ). A series of the CFD simulations were executed with k pol = 2 × 10 5 [N·s/m] and a relaxation time τ varying from 0.00625 to 0.12100 s. Next, for each CFD solution, a specific set of the τ, and K 1 constants (see Equation (15)) that maximizes the R 2 coefficient (15), when comparing the experimental (MR t,rate−exp ) and model (MR t,rate−mod ) values, was found. Table 3 and Figure 15 present the results of the optimization problem calculated after 950 AFM passes.
The optimized solution with a maximum R 2 =0.4759 was found for τ = 0.025 [s] and K 1 = 0.86 × 10 −12 m 3 /J , which leads to the following model: The MR t,rate ( N 1 , v) function is then used to calculate the MR t,model fields for each of the four surfaces of the V-shape artifact, compare them to the MR t,exp fields, and calculate the discrepancy between the simulation and measurements results. Note that Figure 15a presents the experimentally measured MR t fields, Figure 15b presents their CFD equivalents calculated using Equations (2), (9) and (10), and Figure 15c illustrates the Error calculated using Equation (18) (note that the prerequisite v, . χ, µ . χ and N 1 fields can be seen in the Appendix A, Figure A2).  Note that the experimental results correlate well with the ones obtained in Bouland, Urlea, Beaubier, Samoilenko and Brailovski [17] with characteristic "river bed" zones at the specimen tips. Next, comparison between the experimental and modeling results shows similarity in the MR t trends for surfaces A, C and D. However, the calculated surface B does not manifest distinct MR t zones seen from the experiment. Moreover, significant discrepancies between the calculated and experimental values are observed in zones with drastic variations in (a) flow directions (zones with minimum MR) and (b) cross-sections (zones with maximum MR), this reflecting the limitations of the model to simulate non-stationary processes.
Next, the volume-based MR experimental (3D microscope) and numerical results are presented in Figure 16. In terms of local trends, the experimental and modeled MRs range similarly from greater to lower as D, C, A, and B. A good convergence (maximum error of 3.2%) is observed for surfaces A and B, while greater discrepancy is observed for surfaces C and D (maximum error of 15.8%). Globally, as shown in Figure 17, the predicted MR values are situated below the material removed measurements obtained using the scales and the 3D microscope, with a maximum error not exceeding 17%.

Calibration of the Ra(MR t ) Function Using As-Built V-Shape Artifacts
Before starting AFM, the initial wall roughness (Ra 0 ) as a function of the build orientation (α) was evaluated for the V-shape artifact and plotted in Figure 18a. Then, an Ra(MR t ) function was extracted from the experimental measurements of the material removal and surface roughness evolution, and they are shown in Figure 18b. Note that RAW data presenting the Ra(t) and MR t (t) evolutions as functions of the polishing time are shown in the Appendix A ( Figure A3). Surfaces A and C with similar Ra(MR t ) distributions are merged into a single curve. Next, the Ra = f (Ra 0 , MR t ) function is established and a best-fit power function Ra(Ra 0 , MR t ) was then determined to make the Ra values approaching those in the fully polished state Ra pol at the end of the AFM sequence: where a 1 = −1.570; b 1 = 0.500; c 1 = 2.848 µm; d 1 = 0.333; Ra pol = 0.500 µm.  From the analysis of Figure 18b it can be concluded that the Ra(MR t ) distribution depends on the initial surface roughness Ra 0 and varies non-linearly. Surfaces A and C (90 deg) with the same initial Ra 0 possess similar Ra(MR t ) evolutions. All the surfaces follow constantly decreasing roughness improvement rates (dRa/dMR t ) and reach their asymptotic limits at Ra pol ≈ 0.45-0.55 µm.

S-Shape Specimen: Validation of the MR Model
The detailed evolutions of the flow rate and mass removal during the S-shape AFM trial are shown in Figure 19. CFD simulations of an entire AFM setup were run with the calibrated relaxation time τ = 0.025 s, and the slip coefficient k slip varying accordingly to the test sequence (see Appendix A, Table A2). To adjust the slip coefficient for this specific setup, the next steps were taken: (1) We started from the final "Test 7" results, by assuming that the last case corresponded to the completely polished S-shape state. The entire AFM system was considered polished as well. By applying the experimental flow rate at the inlet (Q = 3 × 10 −6 m 3 /s), the back pressure at the outlet (p back = 0 psi), and the calibrated slip coefficient (k pol = 2 × 10 5 [N·s/m]), a solution for the inlet pressure was found (p inlet,CFD = 4.32 × 10 6 Pa ≈ 600 psi). The simulated p inlet,CFD was approximately x2 lower as compared to the actual AFM process (p inlet,AFM = 1200 psi). This difference was attributed to an additional back pressure resulting from the AFM system resistance. (2) For the 1 to 6 cases, for the chamber/reducer/fixture, we kept the same k slip = 2 × 10 5 [N·s/m], while adjusting k slip only for the S-shape specimen, in order to maintain an inlet pressure of p inlet = 4.32 × 10 6 Pa. From Figure 20a,b, a non-uniformity of the MR t distribution is observed (higher/lower MR zones). To compare the experimental (MR t,exp ) and calculated (MR t,model ) values on a point-to-point basis, the Error field was calculated (Figure 20c) (note that these calculations were carried out after 26 AFM passes). It can be observed that the developed model predicts the experimental MR t trends with higher MR t at entrances and larger radii and lower MR t at smaller radii of the channel. For the point cloud considered, the model underpredicted results by~14% with an R 2 of 0.4063, where the biggest error was observed at larger radii and entrances, where drastic variations in (a) flow direction and (b) cross-section occur (similar to the V-shape). The scales, 3D microscope and MR model results were compared globally in terms of the material removed during the AFM process ( Figure 21). The discrepancy between the scales and 3D microscope measurements was~4.4%. Globally, the MR model underpredicted the 3D microscope measurements by~6.2%, which is in line with the discrepancy found during the MR t field analysis. Before starting AFM, an initial wall roughness (Ra 0 ) was evaluated for the S-shape component as a function of the build orientation (α) and is plotted in Figure 22a. Then, the corresponding initial roughness distribution Ra 0 (α) was found as follows: where a 0 = 58.81, b 0 = −0.01311, c 0 = −53.23, d 0 = −0.06065. Based on the calibrated Ra(Ra 0 , MR t ) evolution function, the S-shape Ra 0 (α) and modelled MR t (t) values were used to predict the Ra(t) evolution during the AFM process. It can be seen in Figure 22b that for the 0 • build angle, the MR model significantly overpredicted the experimental values (by up to 300%). However, for the 40 • build angle, the predicted Ra correlated well with the experimental values, and the error did not exceed 23% when compared with the average Ra measurements. Based on the Figures 21 and 22 results, model overpredicts the experiment, but follows the same trends.

Discussion
From the present research, it was concluded that the MR model calibration protocol requires the use of calibration artifacts in two states: pre-polished and as-built. The use of the pre-polished artifacts eliminates the influence of as-built surface roughnesses on the AFM process and allows the calibration of the following parameters of the material removal model: (a) slip coefficient (k slip ), (b) relaxation time (τ) and (c) abrading coefficient (K 1 ). The use of the as-built calibration artifacts allows measuring the surface roughness evolution as a function of the material removal.
Slip coefficient (k slip ). In the literature dealing with CFD modeling of the AFM process, the authors apply no-slip conditions at walls based on the classical principles of fluid dynamics. However, the results of this study demonstrate that k slip is a cornerstone parameter responsible for the velocity field distribution, especially for curved surfaces and channels, and must be carefully calibrated.
Relaxation time (τ) and abrading coefficient (K 1 ). The relaxation time is another cornerstone parameter responsible for the distribution of normal stresses on the surfaces subjected to AFM. Whereas k slip was calibrated using the experimentally determined inlet pressures, τ and K 1 were calibrated simultaneously such that, when comparing the experimentally obtained MR results with those obtained from modeling, a maximum R 2 of 0.4759 and an average error of 35.83% were obtained. Despite this discrepancy between the calculated and measured material removal values observed in zones with drastic variations in flow direction and velocity, a value of the abrading (Preston) coefficient (K 1 = 0.86 × 10 −12 m 3 /J) found in this study corresponds well to the literature data. For example, for the grinding of a hardened steel part using a cast iron lapping tool and diamond suspension, Speich and Börret [34] found a K 1 of 0.76 × 10 −12 m 3 /J, while for the polishing of a tool steel part using a water-glycol compound containing monocrystalline diamond particles, Dambon et al. [35] found K 1 in the (0.33 − 1.00) × 10 −12 m 3 /J range.
The calibrated k slip , τ, and K 1 values were used in the validation study on the S-shape specimens. Results demonstrated that for a specific AFM medium/part combination, there exists a critical roughness (Ra cr ) after which the slip coefficient remains relatively constant ( Figure 23): For the LV-60B medium (EH) and a stainless steel (EOS) part (S-shape), Ra cr ∼ = 10 µm. The results show that~83% of the polishing time corresponded to the calibrated k slip = 2 × 10 5 [N·s/m]. Thus, from a practical standpoint, such a calibrated slip coefficient may be taken as a first approximation during the CFD simulations of the process. However, a significant discrepancy was observed when comparing the experimental and CFD-simulated inlet pressures (p inlet,exp = 1200 psi vs. p inlet,CFD = 600 psi) during AFM of S-shape specimens. When retaining p inlet,CFD = 600 psi as a control parameter resulted in a relatively low average error <14% between the modelling and experimental MR results (both thickness and weight), while using p inlet,CFD = 1200 psi doubled this error. Based on these findings, two groups of possible reasons for this inconsistency were put forward, with the first related to the experimental conditions of the study, while the second pertained to the AFM conditions. (a) Experiment: the use of different setups for calibration and validation, which contributed to creating inequivalent polishing conditions in terms of the inlet and back pressures applied. (b) Modeling: the use of a simplified viscoelastic model and identical laws for the normal and shear viscosities dependences, as well as the use of constant values of the relaxation time and slip coefficient during the entire process duration.
Finally, based on the results of this study, a compensated 3D model (CAD comp ) can be generated by adding the calculated AFM machining allowances (MR t ) to the initial 3D model (CAD initial ). In the example shown in Figure 24, a MATLAB routine generates a compensated surface (STL comp ) based on the initial point cloud (STL initial ) and the predicted MR t field. To this end, each initial vertice of STL initial is displaced by the MR t value along the direction corresponding to an average of the surrounding face normals (averaging is required since one vertice may be attributed to different faces with different normals). Next, a CAD comp is generated using the CAD initial and the previously obtained STL comp using CATIA V5 and Autodesk Inventor through the trivial CAD modeling routines. Using the predicted MR t field, the non-uniform AFM machining allowances are added to the initial CAD and a 3D printing-ready compensated CAD model is generated (Appendix A, Figure A4). This compensated CAD takes into consideration the AFM-related non-uniformity in the material removal, such as higher allowances at the entrances and the larger radii of the channel (Figure 20b).

Conclusions
(1) The developed MR model based on the simplified viscoelastic model (ANSYS Polyflow software) and the calibration methodology using the V-shape calibration artifacts shows an average discrepancy with the experimental results not exceeding 25%, which is deemed acceptable for real-case applications; (2) The slip coefficient (k slip ) and the viscoelastic relaxation time (τ) are two parameters that greatly influence the MR field distribution and require special attention in their determination. It is proposed to calibrate the k slip and τ values using pre-polished V-shape calibration artifacts; (3) A strong dependence of k slip on Ra was demonstrated (both as-built V-shape and Sshape). From the CFD analysis of the S-shape specimens, a critical value of roughness value (Ra cr ) was determined, such that for Ra < Ra cr , k slip could be considered relatively constant; (4) To predict the velocity and normal stress fields, it is recommended to study an entire AFM system by simultaneously controlling the flow rate and the inlet/back pressures. With this approach, the k slip adjustments can be achieved using real operational RAW data. Funding: This research received no external funding.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.
Appendix A