Non-uniform material removal in robotic compliant grinding for flexible free-form surfaces

Abstract

In a convolution material removal process, taking grinding free-form surfaces as an example, the workpiece’s complex shape may lead to dynamic tool–workpiece contact state, and the curved tool path results in an uneven dwell time distribution. These factors contribute to non-uniform material removal (NMR), causing over-grinding or under-grinding in localized areas. This work aims to model NMR accurately and propose a method to enhance material removal uniformity. First, a dynamic tool–workpiece contact model integrating the workpiece’s complex shape, contact force, and the mechanical properties of the tool and the workpiece is proposed by introducing the measured workpiece point cloud. Second, path geodesic curvature is employed to calculate the dwell time distribution. Third, a material removal model that combines the dynamic tool–workpiece contact and the uneven dwell time distribution is introduced. Then, the tool influence function is optimized by adjusting the tool orientation to improve material removal uniformity. Finally, the proposed material removal model and optimization method are validated through experiments, with results showing a remarkable improvement in material removal uniformity using this approach.

Graphical abstract

Keywords

robotic grinding / material removal / tool influence function

1 Introduction

Flexible free-form surfaces refer to free-form surfaces made of flexible materials, which play a key role in the formation of aircraft skins [1,2]. Traditionally, these surfaces are machined by manual grinding, which is inefficient and makes it difficult to meet the high precision requirements [3,4]. With the development of industrial robots, fully automated grinding is expected to improve the efficiency and accuracy of free-form surfaces, which have received much attention [5,6]. The material removal in robotic grinding is a convolution process. In these cases, surface evolution can be understood as the convolution of a tool influence function (TIF) and the dwell time at each path point [7]. The TIF represents the material removed per unit time, influenced by the contact state, and the dwell time is reciprocal of the feed rate [8].

In our previous work [9], we focused on the flexible planar workpiece and proposed a dual flexible contact material removal model to address the issue of inaccurate material removal prediction caused by the flexibility of the workpiece. However, when a grinding tool moves along a curved path on a free-form surface with constant contact force, unexpected grinding on the workpiece surface, such as over-grinding or under-grinding in localized areas, was observed. This issue is referred to as a non-uniform material removal (NMR) problem, which can be attributed to two main factors: (1) dynamic pressure distribution in the contact region, resulting from the workpiece’s complex geometric shape, and (2) uneven distribution of dwell time. Although the dwell time at each path point is usually considered a scalar [10,11], the authors believe that the tool size must be considered because it results in varying dwell times for each microelement of the tool surface within the contact region. For example, when the tool path is curved, the workpiece surface along the positive direction of the geodesic curvature vector experiences a longer grinding time than the negative direction. Therefore, to improve the uniformity of material removal, a NMR model considering the free-form surface shape and the curved tool path is necessary.

Accurately modeling the tool–workpiece contact is crucial to establishing the material removal model [12], which includes the contact geometry, force, and pressure distribution [13]. The Hertz model is a widely used contact theory [14,15] that assumes that two contacting elastic bodies form a localized elliptical contact region with symmetric pressure distribution [16]. In cases such as irregular contact regions in disk grinding or large tool deformations, Hertz contact theory assumptions are violated, making it ineffective. Wang et al. [17] discovered that in robotic disk grinding, the contact region is enclosed by two circular arcs instead of forming a complete ellipse. Therefore, the differential geometric methods are used to model the contact geometry and pressure distribution. In ball-end grinding, Kim and Kim [18] found that the tool–workpiece contact pressure conforms to a Gaussian distribution. Recently, finite element techniques have also been adopted to calculate contact geometry and pressure. Wang et al. [19] proposed a physical simulator aiming to explore the impact of varying contact wheel hardness and contact forces on material removal in belt grinding. Zhang et al. [20] proposed a novel model combining support vector regression with finite element analysis data to enhance the accuracy of contact pressure distribution estimation. In robotic compliant disk grinding for free-form surfaces, Hertz contact theory is no longer applicable because the assumption of an elliptical contact shape is violated. The differential geometry method hardly captures the varying contact geometry when grinding complex free-form surfaces, which may affect the accuracy of the contact pressure distribution model. The finite element method demands large computational resources. Therefore, developing an efficient and accurate contact geometry model that dynamically adjusts to the curved toolpath poses a considerable challenge. After studying the contact model, many researchers have investigated the material removal mechanisms for quantifying the material removal depth. These studies are mainly classified into three categories: analytical methods [21,22], empirical wear equations [23,24], and data-driven methods [25,26]. In this work, we model material removal on the basis of the empirical wear equation because it greatly simplifies the material removal model compared with the other two methods, making it more suitable for subsequent optimization.

In terms of dwell time, a scalar is sufficient to predict the material removal map [17]. The tool’s influence on the workpiece surface correlates directly with the geodesic curvature. Workpiece points located in the positive direction of the geodesic curvature vector will be influenced by more path points, whereas those in the opposite direction will be affected by fewer. However, predicting the complete material removal map takes substantial time, which is not always necessary. Sometimes, modeling the material removal profile (hereinafter referred to as “profile”) at a single path point is more important [27]. In this case, proposing the dwell time of a distribution to compute the profile at a single path point would save much time.

The material removal model can be used to predict the machined surface of the workpiece, facilitating the evaluation of NMR. Furthermore, addressing the NMR issue is crucial. Typically, researchers improve the uniformity of material removal by optimizing the tool path. Some scholars are devoted to generating equidistant tool paths that covers the entire workpiece [28,29]. Notably, these works do not consider the grinding mechanism, which may result in the NMR problem. Some scholars have improved material removal uniformity by optimizing the spacing of tool paths. Zhang et al. [30] proposed the physically uniform coverage. The contact model is utilized to compute the variation in grinding width. Then, they adjusted the spacing between paths, guaranteeing that the grinding ribbons overlap uniformly along neighboring grinding paths. Similarly, Wen et al. [31] estimated the contact shape with a given normal force by differential geometry. Then, they optimized the tool path to achieve minimal overlap considering variations in contact shape as the tool moves on the surface. The above work improves the uniformity of material removal between multiple paths. Additionally, material removal uniformity within individual paths is also crucial, and the most effective way to improve it is by optimizing the TIF. The TIF is influenced by key factors including tool orientation, contact force, and tool rotation speed, with tool orientation primarily determining the shape of the TIF. In many works, the tool orientation is determined by the normal of the workpiece surface, which is not conducive to ensuring uniform material removal [32]. Thus, Gao et al. [33] concentrated on enhancing the pressure distribution uniformity in the contact region by optimizing the tool orientation, contributing to the workpiece’s surface quality and the grinding tool’s longevity. Their model is insufficient to describe the workpiece shape’s effect on dynamic tool–workpiece contact accurately.

Based on the analysis above, the shape of the free-form surface can lead to dynamic tool–workpiece contact, changing contact pressure distribution. In addition, the geodesic curvature of the path leads to an uneven distribution of dwell time. Both phenomena lead to NMR within an individual tool path, and a few studies have fully elucidated this mechanism. Therefore, this work focuses on the following two aspects:

● NMR model. An accurate material removal model reveals the mechanism for NMR and guides the subsequent optimization.

● TIF optimization. Optimizing the TIF by adjusting the tool orientation to compensate for the NMR caused by dynamic contact geometry and uneven dwell time distribution.

This paper is organized as follows: The precise model of the material removal is presented in Section 2. Section 3 deduced the objective function and optimization methods. The proposed models and optimization methods are validated in Section 4. The conclusions are given in Section 5.

2 Modeling

Accurate modeling of the NMR process helps quantify the material removal uniformity. In this section, we model the two main factors affecting the material removal uniformity: the dynamic contact state and the uneven dwell time distribution. Then, the material removal map, which describes the complete workpiece surface after machining, and the material removal profile at a single path point is derived.

2.1 Dynamic tool–workpiece contact model considering workpiece complex surface shape

1) Modeling of contact geometry

Workpiece shape is commonly represented by a point cloud. The coordinates of an individual point on the workpiece are {\mathit{w}}_i=[x_i,y_i,z_i,1{]}^{\mathrm{\top }}, i=1,\dots ,n_w. These points constitute the matrix \mathit{W}=[{\mathit{w}}_1,{\mathit{w}}_2,\dots ,{\mathit{w}}_{n_w}]. The tool–workpiece contact geometry is related to the tool pose on the workpiece surface; thus, a tool path is needed to determine the initial tool position. The tool path consists of a series of ordered points: {\mathit{q}}_k=[x_k,y_k,z_k{]}^{\mathrm{\top }}, k=1,\dots ,n_q. Three reasonable assumptions are made here to model the contact geometry at path point {\mathit{q}}_k. (1) The contact geometry is only related to the normal force, and the influence of tangential force on contact geometry is ignored. (2) The deformation is elastic, with no plastic deformation occurring. (3) The material properties of the tool and the workpiece are considered isotropic. On the basis of these assumptions, the tool and workpiece are modeled as a series of springs aligned along the z-axis of \{B_k\}. This spring-based model effectively represents the relationship between axial deformation and contact pressure, providing a simplified yet accurate depiction of the contact behavior. At this point, the tool–workpiece contact geometry can be interpreted as the z-coordinate of the workpiece point cloud in the frame \{B_k\}:

(1)u_{k,i}=\left|{}^{B_k}{z}_i\right|,i\in I_k,

where \left|\cdot \right| represents the absolute value operation, u_{k,i} is the contact depth at workpiece point {\mathit{w}}_i when the tool moves to path point {\mathit{q}}_k, {}^{B_k}{z}_i is the z-coordinate of the workpiece point cloud in the frame B_k, I_k is the set of index i of all activated workpiece surface points (defined as points within the contact region). Given that the tool is a cylinder, the activated workpiece points can be identified as those lying inside the tool:

(2)I_k=\{i\in \left\{1,\dots ,n_w\right\}|{\left({}^{B_k}{x}_i-r\right)}^2+{\left({}^{B_k}{y}_i\right)}^2⩽r^2,{}^{B_k}{z}_i⩾0\},

where r denotes the tool radius.

The workpiece point cloud \mathit{W} is fixed in world frame \{C\}. The tool moves within frame \{C\}, with the frame \{B_k\} representing the tool’s pose. When modeling the contact geometry, the coordinates of the workpiece point cloud in frame \{B_k\} must be calculated: {}^{B_k}\mathit{W}={}_{A_k}^{B_k}\mathit{T}{}_C^{A_k}\mathit{T}{}^C\mathit{W}. As illustrated in Fig. 1, frame \{A_k\} is established to represent the initial tool position at the path point {\mathit{q}}_k. In this frame, the z-axis aligns with the surface normal direction {\mathit{n}}_k\in {\mathbb{R}}^3. The x-axis lies in the plane formed by {\mathit{n}}_k and the vector {\mathit{\tau }}_k={\mathit{q}}_k-{\mathit{q}}_{k-1}. It is perpendicular to {\mathit{n}}_k and indicates the direction of tool motion. The y-axis is then determined in accordance with the other two axes. Hence, the transformation matrix {}_C^{A_k}\mathit{T}\in \mathrm{S}\mathrm{E}(3) from the world frame \{C\} to the tool frame \{A_k\} is denoted as

(3){}_C^{A_k}\mathit{T}=\left[\begin{array}{cccc}{\displaystyle \frac{{\mathit{n}}_k\times {\mathit{\tau }}_k\times {\mathit{n}}_k}{{\Vert {\mathit{n}}_k\times {\mathit{\tau }}_k\times {\mathit{n}}_k\Vert }_2}}& {\displaystyle \frac{{\mathit{n}}_k\times {\mathit{\tau }}_k}{{\Vert {\mathit{n}}_k\times {\mathit{\tau }}_k\Vert }_2}}& {\displaystyle \frac{{\mathit{n}}_k}{{\Vert {\mathit{n}}_k\Vert }_2}}& {\mathit{q}}_k\\ 0& 0& 0& 1\end{array}\right],

where {\Vert \cdot \Vert }_2 means the Euclidean norm. Frame \{B_k\} is the tool frame when grinding. The tool axis is usually obtained by tilting the z-axis of frame \{A_k\} to improve the removal of chips and prevent zero linear velocity in the tool center area. The associated rotation matrix is denoted as {}_{A_k}^{B_k}\mathit{R}. Additionally, grinding is a force-controlled process. When the tool grinds the workpiece with a given normal force F, it undergoes elastic deformation [34,35]. The maximum contact depth in the current contact region, i.e., z-direction displacement between frame \{A\} and frame \{B\}, is denoted as h_k. Therefore, the transformation matrix {}_{A_k}^{B_k}\mathit{T} from \{A_k\} to \{B_k\} is

(4){}_{A_k}^{B_k}\mathit{T}=\left[\begin{array}{cc}{}_{A_k}^{B_k}\mathit{R}& {[0,0,-h_k]}^{\mathrm{\top }}\\ \mathbf{0}& 1\end{array}\right].

2) Pressure distribution with a given contact force

The nonlinear dual flexible contact force model proposed in our previous work [9] is used to describe the relationship between contact depth and contact pressure:

(5)p_{k,i}=u_{k,i}^{\alpha }\hat{E},

where p_{k,i} signifies the contact pressure, \hat{E} denotes the equivalent contact coefficient, and \alpha represents the nonlinear power index. The force at each point is the product of the pressure and the area, and then the normal contact force F_k on the tool is the sum of the forces at all points within the contact region:

(6)F_k=\sum _{i\in I_k}{p}_{k,i}{s}_{k,i},

where s_{k,i} indicates the area associated with point w_i.

The workpiece point cloud does not contain area s_{k,i}. In our case, the point cloud is processed through voxel grid downsampling with a voxel size of l. Following this, in frame \{B_k\}, s_{k,i}=l^2\mathrm{c}\mathrm{o}\mathrm{s}{\vartheta }_{k,i}/\phantom{\mathrm{c}\mathrm{o}\mathrm{s}{\vartheta }_{k,i}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i, where {\theta }_i=\arccos \left([0,0,1]{\mathit{n}}_i\right) indicates the angle between the normal vector of workpiece point {\mathit{w}}_i and the z-axis of the world frame \{C\}, and {\vartheta }_{k,i}=\arccos [\frac{1}{2}\mathrm{t}\mathrm{r}({}_{A_k}^{B_k}\mathit{T}{}_C^{A_k}\mathit{T})-1] represents the angle between the z-axis of the world frame \{C\} and tool frame \{B_k\}. The maximum contact depth with a given contact force F_{\mathrm{c}\mathrm{m}{\mathrm{d}}_k} is obtained by solving the following optimization problem:

(7)h_k^{\ast }=\underset{h_k}{\arg min}{\Vert F_{\mathrm{c}\mathrm{m}{\mathrm{d}}_k}-F_k(h_k)\Vert }_2.

In this work, Eq. (7) is solved using the bisection method. Then, by substituting h_k^{\ast } into Eqs. (4) and (5), the contact geometry and contact pressure are derived.

2.2 Material removal model based on tool path geodesic curvature

1) Modeling of material removal map

Within the contact region at path point {\mathit{q}}_k, we assume that material removal is proportional to the contact pressure and sliding velocity, ignoring the effects of friction and temperature on material removal. Therefore, the material removal depth can be determined using the Preston equation [36]:

(8)d_{k,i}=K{p}_{k,i}{v}_{k,i}{t}_k,

where K is the wear coefficient; d_{k,i} is the material removal depth; p_{k,i} is the contact pressure; t_k is the dwell time scalar determined by path point spacing and planned feed rate f_k, t_k={\Vert {\mathit{q}}_k-{\mathit{q}}_{k-1}\Vert }_2/\phantom{{\Vert {\mathit{q}}_k-{\mathit{q}}_{k-1}\Vert }_2{f}_k}{f}_k; and v_{k,i} is the sliding velocity determined by the spindle speed \omega, given by v_{k,i}=2\pi \omega [{\left({}^{B_k}{x}_i-r\right)}^2+{\left({}^{B_k}{y}_i\right)}^2{]}^{1/2}. The subscripts of (\cdot {)}_{k,i} denote the contribution to workpiece point i when the tool moves to the path point k.

Material removal in robotic disk grinding can be considered a convolution process. The tool grinds a workpiece surface point at many different path points. The material removal matrix {}^{B_k}{\mathit{M}}_k is proposed to represent the material removed at each path point:

(9)\begin{array}{c}{}^{B_k}{\mathit{M}}_k=[{\mathit{m}}_{k,1},{\mathit{m}}_{k,2},\dots ,{\mathit{m}}_{k,n}],\\ {\mathit{m}}_{k,i}=\{\begin{array}{ll}{[0,0,0,1]}^{\mathrm{\top }},& i\notin I_k,\\ [0.3em]{[0,0,d_{k,i},1]}^{\mathrm{\top }},& i\in I_k.\end{array}\end{array}

Given that these material removal matrices are calculated in frame \{B_k\}, they must be transformed to the world frame \{C\}. Then, the material removal map, i.e., the finished workpiece, is obtained by subtracting all material removal matrices from the initial workpiece point cloud:

(10){}^C{\mathit{W}}_{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{s}\mathrm{h}\mathrm{e}\mathrm{d}}={}^C\mathit{W}-\sum _{k=1}^m{}^C{\mathit{M}}_k.

The colormap can be represented by the sum of the material removal matrices \sum _{k=1}^m{}^C{\mathit{M}}_k.

2) Material removal profile based on tool path geodesic curvature

The material removal map provides a comprehensive visualization of the entire finished workpiece on the point cloud. Additionally, the profile {\gamma }_k(y) at each path point is important because it helps evaluate the local uniformity of material removal.

Typically, as shown in Fig. 2, a profile requires calculating many material removal matrices which contribute to {\gamma }_k(y). Then, the material removal map is drawn on the basis of numerical convolution, and the {\gamma }_k(y) is obtained by interpolation. The computation of multiple material removal matrices requires considerable time. To this end, a weight factor t_{w_{k,i}} based on geodesic curvature is proposed to improve computational efficiency. The use of multiple removal matrices is replaced by the weighted sum of a single matrix to generate the material removal map.

As shown in Fig. 3, T_{{\mathit{q}}_k} is the tangent plane of the workpiece at point {\mathit{q}}_k. The motion of the tool in T_{{\mathit{q}}_k} can be interpreted as a rotation around the center of geodesic curvature {\mathit{r}}_k. The dwell time of each point is related to its distance from the center of the geodesic curvature circle. Therefore, the weighting factor is

(11)t_{w_{k,i}}=t_k\frac{{\rho }_k}{{\Vert {}^{A_k}{\mathit{w}}_i-{\mathit{r}}_k\Vert }_2},

where {}^{A_k}{\mathit{w}}_i is the coordinate of workpiece point {\mathit{w}}_i in frame \{A_k\}, and {\rho }_k is the geodesic curvature radius. {\rho }_k and {\mathit{r}}_k are determined by the geodesic curvature vector {\mathit{\kappa }}_g:

(12){\rho }_k=\frac{1}{{\Vert {\mathit{\kappa }}_{gk}\Vert }_2},{\mathit{r}}_k=\frac{{\rho }_k{\mathit{\kappa }}_{gk}}{{\Vert {\mathit{\kappa }}_{gk}\Vert }_2}=\frac{{\mathit{\kappa }}_{gk}}{{\Vert {\mathit{\kappa }}_{gk}\Vert }_2^2}.

Given that the curvature of the path point {\mathit{q}}_k is {\mathit{\kappa }}_{qk}, the geodesic curvature vector {\mathit{\kappa }}_{gk} is denoted as

(13){\mathit{\kappa }}_{gk}={\mathit{\kappa }}_{qk}{\mathit{n}}_k^{\mathrm{\top }}{\mathit{n}}_k-{\mathit{n}}_k.

After obtaining the weighting factors, the weighted material removal depth is

(14)d_{w_{k,i}}=d_{k,i}{t}_{w_{k,i}}.

To obtain the profile, d_{w_{k,i}} is re-interpolated along the x-axis and y-axis in the frame \{B_k\}, denoted as d_{w_{k,x,y}}. This step is realized in MATLAB^® by interpolation, with a spacing of l. Then, the profile {\gamma }_k(y) is obtained by summing d_{w_{k,x,y}} along the x-axis:

(15){\gamma }_k(y)=l\sum _{x=0}^{2r}{d}_{w_{k,x,y}}.

3 TIF optimization

An optimal TIF can improve material removal uniformity by changing the pressure distribution, and an optimized pressure distribution can compensate for NMR resulting from the dynamic tool–workpiece contact state and uneven dwell time. Once the workpiece and tool are determined, the tool orientation is the primary factor determining the TIF; therefore, it is selected as the optimization variable.

3.1 Tool orientation description

The tool orientation in \{B_k\} can be conveniently represented by a rotation vector {\mathit{\varphi }}_k=[{\varphi }_{k,1},{\varphi }_{k,2},0{]}^{\mathrm{\top }}\in {\mathbb{R}}^3. Here, the third term of {\mathit{\varphi }}_k is 0 because the rotation of the tool around itself does not affect the pressure distribution. A linear skew-symmetric operator, denoted as (\cdot {)}^{\wedge }, is the mapping from the vector space to Lie algebras \mathfrak{s}\mathfrak{o}(3), namely,

(16){\mathit{\varphi }}_k^{\wedge }=\left[\begin{array}{ccc}0& 0& {\varphi }_{k,2}\\ 0& 0& -{\varphi }_{k,1}\\ -{\varphi }_{k,2}& {\varphi }_{k,1}& 0\end{array}\right].

Each matrix {\mathit{R}}_k\in \mathrm{S}\mathrm{O}(3) and its Lie algebra {\mathit{\varphi }}_k^{\wedge } have exponential and logarithm mapping: {\mathit{R}}_k=\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{\varphi }}_k^{\wedge }), {\mathit{\varphi }}_k^{\wedge }=\mathrm{l}\mathrm{o}\mathrm{g}({\mathit{R}}_k).

3.2 Analytical expression of NMR problem

The modeling work in Section 3 is suitable for calculating the global material removal map but is inefficient for optimizing tool orientation. As shown in Fig. 4, all the points in the possible contact region, i.e. ({}^{A_k}{x}_{k,i}-r{)}^2+({}^{A_k}{y}_{k,i}{)}^2⩽r^2, are selected in frame \{A_k\} to reduce computation time. In addition, these randomly distributed points are rasterized to improve computational efficiency. The rasterized points have n levels for the x and y coordinates, denoted as x_j^{\mathrm{l}\mathrm{e}\mathrm{v}} and y_j^{\mathrm{l}\mathrm{e}\mathrm{v}}, j=1,\dots ,n. The sub-region {\mathit{W}}_k=[{\mathit{w}}_{k,1},{\mathit{w}}_{k,2},\dots ,{\mathit{w}}_{k,n^2}]\in {\mathbb{R}}^{3\times n^2} is formed by these rasterized points {\mathit{w}}_{k,i}=[x_{k,i},y_{k,i},z_{k,i}{]}^{\mathrm{\top }}, i=1,\dots ,n^2.

At path point {\mathit{q}}_k, the corresponding sub-region {\mathit{W}}_k is adopted to deduce an objective function to evaluate the removal uniformity problem. This objective function describes the sum of the differences in the weighted material removal depths on both sides of the tool path:

(17)f({\mathit{\varphi }}_k)={\Vert \left({\mathit{d}}_k({\mathit{\varphi }}_k)\mathit{\Lambda }\right)\odot \mathit{e}\Vert }_2,

where \odot is the Hadamard product, {\mathit{d}}_k is the weighted material removal (given that all subsequent material removal depths are weighted, {\mathit{d}}_{w_k} is abbreviated as {\mathit{d}}_k). \mathit{e} is the weight, and its element e_i is determined by the y-coordinate in frame \{B_k\} of the point e_i=\frac{1}{\sigma \sqrt{2\pi }}\exp [\frac{{(r-y_i)}^2}{2{\sigma }^2}]. \mathit{\Lambda }\in {\mathbb{R}}^{n^2\times n} is an operator used to calculate the difference in material removal depth:

(18){\lambda }_{i,j}=\{\begin{array}{ll}1,& \mathrm{i}\mathrm{f}{y}_{k,i}=y_j^{\mathrm{l}\mathrm{e}\mathrm{v}},\\ -1,& \mathrm{i}\mathrm{f}{y}_{k,i}=-y_j^{\mathrm{l}\mathrm{e}\mathrm{v}},\\ 0,& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\end{array}

In sub-region {\mathit{W}}_k, according to Section 2, the weighted material removal can be derived as

(19)\begin{array}{c}{\mathit{d}}_k({\mathit{\varphi }}_k)={\mathit{u}}_k({\mathit{\varphi }}_k{)}^{\alpha }\odot {\mathit{b}}_k,\\ {\mathit{b}}_k={\mathit{v}}_k\odot {\mathit{t}}_{w_k}K\hat{E}.\end{array}

Here, (\cdot {)}^{\alpha } is the element-wise power. The z-coordinate of the workpiece point cloud when in contact is given as

(20){\mathit{z}}_k({\mathit{\varphi }}_k)=[001]({\mathit{R}}_k({\mathit{\varphi }}_k){\mathit{W}}_k).

An analytical expression of {\mathit{d}}_k helps improve optimization efficiency. Given that the speed of analytical differentiation is much faster than numerical differentiation, an approximation of the modified sigmoid function is employed in place of the Boolean evaluation in Eq. (2):

(21){\mathit{u}}_k({\mathit{\varphi }}_k)=\frac{{\mathit{z}}_k({\mathit{\varphi }}_k)}{1+\exp [-2\beta {\mathit{z}}_k({\mathit{\varphi }}_k)]},

where \beta is a tunable parameter used to adjust the curvature of the function. A larger value of \beta reduces the approximation error. In our work, we set \beta = 1000.

3.3 Optimization algorithm

To avoid sharp fluctuation of the axis orientation, we defined a dynamic penalty function at iteration t:

(22)g_s^t({\mathit{\varphi }}_k^t)=\Vert \ln [\exp {({\mathit{\varphi }}_k^{t\wedge })}^{\mathrm{\top }}{\mathit{R}}_L]{\Vert }_2,

where {\mathit{R}}_L={\mathit{R}}_{k-1}^{t-1}\exp [\frac{1}{2}\ln ({\mathit{R}}_{k-1}^{t-1\mathrm{\top }}{\mathit{R}}_{k+1}^{t-1})] represents the interpolation at the midpoint between the rotation matrices {\mathit{R}}_{k-1}^{t-1} and {\mathit{R}}_{k+1}^{t-1} at iteration t-1 [37]. {\Vert \ln ({\mathit{R}}_1^{\mathrm{\top }}{\mathit{R}}_2)\Vert }_2 is the angular distance between {\mathit{R}}_1 and {\mathit{R}}_2 in \mathrm{S}\mathrm{O}(3) [38].

During grinding, excessive tilt angles result in excessive radial force on the spindle, whereas a too small tilt angle causes unintentional interference between the tool and the workpiece. Hence, u{b}_1,l{b}_1,u{b}_2,l{b}_2 represent the upper bounds and lower bounds of {\varphi }_1 and {\varphi }_2, respectively. In this work, u{b}_1,l{b}_1,u{b}_2, and l{b}_2 are set to −3°, 3°, 2°, and 20°, respectively. These values were determined on the basis of the mechanical limits of the grinding equipment and the geometric characteristics of our workpiece.

On the basis of the analysis above, the objective function is denoted as

(23)\underset{{\mathit{\varphi }}_k}{\arg min}\sum _{k=2}^{n_q-1}\left(f({\mathit{\varphi }}_k)+\mu g_s({\mathit{\varphi }}_k)\right),

(24)\mathrm{s}.\mathrm{t}.\{\begin{array}{c}l{b}_1⩽{\varphi }_{k,1}⩽u{b}_1,\\ l{b}_2⩽{\varphi }_{k,2}⩽u{b}_2,\end{array}

where \mu is the weight of the smoothness penalty. A small \mu means achieving material removal uniformity to a greater extent, but it may result in fluctuations in tool orientation. Conversely, a large \mu means more smoothness but less uniformity. In our work, a value of \mu =10 was found to achieve a remarkable improvement in path smoothness while maintaining acceptable material removal uniformity.

The iterative optimization process on a manifold \mathcal{M} involves addressing perturbations within a vector space and then projecting them to the manifold \mathcal{M} until convergence. In \mathrm{S}\mathrm{O}(3), the perturbation is defined as the vector space \mathit{\psi }\in {\mathbb{R}}^3 of its Lie algebra \mathfrak{s}\mathfrak{o}(3) [39]. In this case, the perturbation {\mathit{\psi }}_k is obtained by the negative gradient direction. According to the Baker-Campbell-Hausdorff formula, the update step is denoted as {\mathit{R}}_k^{t+1}=\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{\psi }}_k^{\wedge }){\mathit{R}}_k^t. One iteration is completed by updating all rotation matrices of path points, except for the first and last points, then iteration is continued until convergence. The tool–workpiece contact is coupled with the tool rotation. The contact depth h_0 must be calculated before each differentiation to update the tool–workpiece contact under force control. The pseudocode and flowchart of the algorithm are provided in Algorithm 1 (Table 1) and Fig. 5, respectively.

3.4 Computational cost analysis

The input size of our algorithm is the number of path points n_q. For each path point, the objective function f({\mathit{\varphi }}_k) and smoothness penalty g_s({\mathit{\varphi }}_k) are computed. Traversing all path points constitutes one iteration, and the algorithm runs for a total of K_{\mathrm{m}\mathrm{a}\mathrm{x}} iterations. The time complexity of the algorithm is \mathrm{O}(n_q\cdot K_{\mathrm{m}\mathrm{a}\mathrm{x}}). We set the number of iterations to K_{\mathrm{m}\mathrm{a}\mathrm{x}}=50 and conducted several tests with input sizes of n_q=500,1000,1500,2000 to demonstrate the algorithm complexity. Tests were performed on a workstation with an Intel^® Core^TM i7-13700F CPU and NVIDIA^® GeForce RTX^TM 4060 Ti GPU. The runtime results are summarized in Table 2, where we define the computational efficiency metric: \mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}={\displaystyle \frac{n_q\cdot K_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{n}\mathrm{d}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}. The ratio remains nearly constant during the tests, indicating that the efficiency of the algorithm does not degrade as the input size increases. Additionally, the loss function consistently converged to a similar range in all trials, demonstrating robustness to input scale.

4 Verification and analysis

Grinding experiments are conducted to validate the proposed model’s accuracy and the optimization methods’ effectiveness. In Fig. 6, an industrial robot with a flexible force-controlled device is used for grinding. The workpieces are measured using a line laser 3D measurement equipment with an error margin of less than 10 μm. The material removal map is obtained by calculating the distance between two point clouds.

The workpiece is a 5 mm-thick ethylene propylene diene monomer adhered to an aluminum alloy free-form surface. Its Young’s modulus is 0.186 MPa, and Poisson’s ratio is 0.2. The sandpaper is 75 mm in diameter and 180# in grit size, and the abrasive is made of aluminum oxide. The new sandpaper is adopted to grind an aluminum workpiece for 1 min before regular use to prevent the sandpaper’s wear from affecting the grinding result; hence, it can reach a stable wear stage. In all experiments, the spindle speed is 1000 r/min, the contact force is 10 N, and the feed rate is 10 mm/s.

4.1 Verification of the dynamic tool–workpiece contact model

Grinding experiments on a free-form surface are conducted to verify the suggested model of dynamic contact. The varying shape of the workpiece when the tool moves along the path leads to the dynamic contact geometry, which in turn alters the distribution of contact pressure. A straight tool path is chosen to avoid the influence of uneven dwell time distribution on material removal. In the experiments, the generation of each path point follows the same procedure, where the position is determined by the workpiece point, and the tool axis is obtained by tilting the normal vector of the workpiece point to a certain angle. In addition, the dwell time, contact force, and spindle speed of each contact point remain consistent. In other words, all factors except for the surface shape of the workpiece are consistent.

The experiments are divided into two parts. The first part, shown in Figs. 7(a)–7(c), is grinding with discrete path points. It aims to verify the accuracy of our model in predicting the tool–workpiece contact state at each path point. It shows that grinding at different path points using identical parameters results in varying material removal depths and regions. The second part, shown in Figs. 7(d)–7(f), is grinding with continuous paths, which is used to validate the proposed material removal model. The modeling material removal map also agrees with the experiments and reflects the varying contact geometry and material removal depth when grinding free-form surfaces. The material removal map is symmetric about the y=0 plane because the workpiece and the tool path are also symmetric. This demonstrates that the influence of complex workpiece geometry on material removal follows an underlying law, which our model has accurately captured.

4.2 Verification of the material removal profile model

The uneven dwell time distribution is validated by grinding experiments with different geodesic curvatures. A plane workpiece is chosen to avoid the effects of dynamic tool–workpiece contact on material removal.

Figures 8(a) and 8(b) show the modeled and experimental material removal map, respectively. The geodesic curvature of the paths is 0.02, 0.01, 0.003, and 0 mm⁻¹ (from left to right, respectively). Figure 8(c) illustrates the results of the existing (based on numerical convolution) and proposed (based on geodesic curvature) models for the profile. The material removal uniformity increases as the geodesic curvature decreases. The RMSE values between the proposed method and the existing method with experimental results are 0.0915 and 0.0950 mm, respectively. Both models can clearly predict the profile well. As the path curvature decreases, the red dashed line (indicating the position of the maximum material removal profile) and the black solid line (representing the location of the tool path) get closer, indicating that the uniformity of material removal is improving.

To illustrate the efficiency advantage, we calculated the elapsed time for the two methods of profile calculation. The existing method based on numerical convolution requires drawing a complete material removal map and then obtaining the profile through interpolation. The time consumed by this method is related to the number of path points. In our example, with a tool diameter of 75 mm and a path point spacing of 0.5 mm, calculating a contour requires the computation of 150 material removal matrices and the corresponding material removal map, taking a total of 3.339 s. The proposed method only requires calculating the material removal matrix for a single path point, then using the geodesic curvature of the current path point to calculate the weight, and summing to obtain the material removal profile. It takes 0.042 s, reducing the time by 98.7%.

4.3 Verification of TIF optimization

To validate the effectiveness of the proposed tool orientation optimization method, a spiral tool path with changing curvature was planned on the workpiece mentioned in Section 4.1. This experiment simultaneously considers the effects of dynamic tool–workpiece contact and uneven dwell time distribution on material removal, showing generality. Then, the tool orientation is optimized using the method described in Section 3, and machining experiments are conducted using the tool paths before and after optimization.

Figure 9(a) depicts the optimized tool tilt angle. The initial values of {\varphi }_1 and {\varphi }_2 are 0° and 10°, respectively. After optimization, the tilt angle dynamically adapts during the grinding process to maintain uniform material removal. Figure 9(b) displays the objective function value. Before optimization, the objective function exhibited drastic variations, with a maximum value of 1.2. By contrast, with the optimization method proposed in this study, the objective function is controlled under 0.25, and it reveals a notable decrease (from 0.39 to 0.14, 64% on average). Within the 0–500 interval, a slight decrease in the objective function is observed at certain path points, attributed to the trade-off between global and local optima. Despite this localized deviation, the proposed method achieves a substantial improvement in overall material removal uniformity.

Figure 10 shows the grinding results. The proposed model can accurately predict material removal in free-form surface disk polishing. Moreover, many areas with NMR in Fig. 10(b) have been notably improved in Fig. 10(e). The proposed tool orientation optimization method remarkably enhanced the uniformity of material removal.

To demonstrate the advantages of our method more specifically, we provide the profiles of six positions on the path. As shown in Fig. 11, the blue curves represent the profiles of material removal before optimization. These curves reveal the inconsistent removal depth on either side of each path point and considerable variations in the removal profiles across different path points. After optimization, each profile (orange curves) displays a consistent shape, peaking in the center and tapering on both sides, thereby remarkably improving the uniformity of material removal. The figure also displays the value of the objective function f, demonstrating a considerable reduction after optimization.

5 Conclusions

The paper presents a precise modeling of NMR due to the complex workpiece shape and the curved tool path during robotic disk grinding of free-form surfaces and proposes a tool orientation optimization method, greatly improving the material removal uniformity. The main contributions are as follows:

● A dynamic tool–workpiece contact model integrating the workpiece’s complex shape, contact force, and the mechanical properties of the tool and the workpiece is proposed by introducing the measured workpiece point cloud. It can accurately predict the dynamic contact geometry and contact pressure distribution during free-form surface grinding.

● A material removal profile model based on geodesic curvature is introduced. Compared with traditional methods reliant on numerical convolution, the proposed model can considerably reduce computational time without compromising prediction accuracy.

● An analytical metric for evaluating NMR is introduced. The TIF is optimized by adjusting the tool orientation to improve the material removal uniformity by minimizing this metric.

This work focuses on optimizing the uniformity of material removal within individual paths. The width of the grinding path is not consistent, which poses a challenge to path planning. In the future, we will address this issue by optimizing the spacing between paths to achieve deterministic grinding of free-form surfaces.

6 Conflict of Interest

The authors declare no competing financial interests or personal relationships that may have influenced the work reported in this paper.

References

[1] Xin C B , Gu Y Z , Li M , Luo J , Li Y X , Zhang Z G . Experimental and numerical study on the effect of rubber mold configuration on the compaction of composite angle laminates during autoclave processing. Composites Part A: Applied Science and Manufacturing, 2011, 42(10): 1353–1360

[2] Wang G , Li W L , Jiang C , Zhu D H , Li Z W , Xu W , Zhao H . Trajectory planning and optimization for robotic machining based on measured point cloud. IEEE Transactions on Robotics, 2022, 38(3): 1621–1637

[3] Ma K , Han L , Sun X , Liang C , Zhang S , Shi Y , Wang X . A path planning method of robotic belt grinding for workpieces with complex surfaces. IEEE/ASME Transactions on Mechatronics, 2020, 25(2): 728–738

[4] Ochoa H , Cortesão R . Impedance control architecture for robotic-assisted mold polishing based on human demonstration. IEEE Transactions on Industrial Electronics, 2022, 69(4): 3822–3830

[5] Zhu W L , Beaucamp A . Compliant grinding and polishing: A review. International Journal of Machine Tools & Manufacture, 2020, 158: 103634

[6] Bari P , Kilic Z M , Law M , Wahi P . Rapid stability analysis of serrated end mills using graphical-frequency domain methods. International Journal of Machine Tools & Manufacture, 2021, 171: 103805

[7] Yamato S, Sencer B, Beaucamp A. Tool path planning and feed scheduling for time-dependent processes. In: Tolio T (ed). CIRP Novel Topics in Production Engineering: Volume 1. Lecture Notes in Mechanical Engineering. Cham: Springer, 2024, 185–231

[8] Wang R , Li Z L , Fan Z , Zhang X Q , Ren M J , Zhu L M . Surrogate-model-based dwell time optimization for atmospheric pressure plasma jet finishing. International Journal of Mechanical Sciences, 2024, 281: 109547

[9] Xu B Z , Zhang X J , Yang Z Y , Wang J Y , Yan S J , Ding H . Dual flexible contact material removal model for robotic disk grinding. Journal of Manufacturing Processes, 2024, 124: 867–876

[10] Yamato S , Sencer B , Beaucamp A . Concurrent process and feedrate scheduling with convoluted basis functions and its application to fluid jet polishing. International Journal of Machine Tools & Manufacture, 2024, 197: 104135

[11] Han Y J , Zhu W L , Zhang L , Beaucamp A . Region adaptive scheduling for time-dependent processes with optimal use of machine dynamics. International Journal of Machine Tools & Manufacture, 2020, 156(April): 103589

[12] Han Y J , Zhang H Y , Yu M H , Yang J Z , Qian L M . Simulation model optimization for bonnet polishing considering consistent contact area response. Frontiers of Mechanical Engineering, 2024, 19(4): 27

[13] Zhou H Y , Zhao H , Li X F , Xu Z R , Ding H . Accurate modeling of material removal depth in convolutional process grinding for complex surfaces. International Journal of Mechanical Sciences, 2024, 267: 109005

[14] Qu S , Wang Z X , Zhang C , Ma Z L , Zhang T Q , Chen H , Wang Z , Yu T B , Zhao J . Material removal profile prediction and experimental validation for obliquely axial ultrasonic vibration-assisted polishing of K9 optical glass. Ceramics International, 2021, 47(23): 33106–33119

[15] Li L F , Ren X K , Feng H J , Chen H B , Chen X Q . A novel material removal rate model based on single grain force for robotic belt grinding. Journal of Manufacturing Processes, 2021, 68: 1–12

[16] Hertz H. Ueber die Berührung fester elastischer Körper. In: Crelle A L, Borchardt W C, Schellbach (eds). Band 92. Journal für die reine und angewandte Mathematik. Berlin: De Gruyter, 1882, 156–171 (in German)

[17] Wang Q H , Liang Y J , Xu C Y , Li J R , Zhou X F . Generation of material removal map for freeform surface polishing with tilted polishing disk. The International Journal of Advanced Manufacturing Technology, 2019, 102(9–12): 4213–4226

[18] Kim D W , Kim S W . Static tool influence function for fabrication simulation of hexagonal mirror segments for extremely large telescopes. Optics Express, 2005, 13(3): 910–917

[19] Wang C , Wu Y R , Liao H Z , Deng C , Luo J , Huang Y . Influence of contact force and rubber wheel hardness on material removal in abrasive belt grinding investigated by physical simulator. Precision Engineering, 2022, 78: 70–78

[20] Zhang X , Kuhlenkötter B , Kneupner K . An efficient method for solving the Signorini problem in the simulation of free-form surfaces produced by belt grinding. International Journal of Machine Tools & Manufacture, 2005, 45(6): 641–648

[21] Yang Z Y , Chu Y , Xu X H , Huang H J , Zhu D H , Yan S J , Ding H . Prediction and analysis of material removal characteristics for robotic belt grinding based on single spherical abrasive grain model. International Journal of Mechanical Sciences, 2021, 190: 106005

[22] Zhu W L , Yang Y , Li H N , Axinte D , Beaucamp A . Theoretical and experimental investigation of material removal mechanism in compliant shape adaptive grinding process. International Journal of Machine Tools & Manufacture, 2019, 142: 76–97

[23] Tsai M J , Huang J F , Kao W L . Robotic polishing of precision molds with uniform material removal control. International Journal of Machine Tools & Manufacture, 2009, 49(11): 885–895

[24] Ke X L , Yu Y H , Li K S , Wang T Y , Zhong B , Wang Z Z , Kong L B , Guo J , Huang L , Idir M , Liu C , Wang C J . Review on robot-assisted polishing: Status and future trends. Robotics and Computer-Integrated Manufacturing, 2023, 80: 102482

[25] Song Y X , Lv H B , Yang Z H . An adaptive modeling method for a robot belt grinding process. IEEE/ASME Transactions on Mechatronics, 2012, 17(2): 309–317

[26] Pandiyan V , Murugan P , Tjahjowidodo T , Caesarendra W , Manyar O M , Then D J H . In-process virtual verification of weld seam removal in robotic abrasive belt grinding process using deep learning. Robotics and Computer-Integrated Manufacturing, 2019, 57: 477–487

[27] Yang Z Y , Xu X H , Li J , Zhu D H , Yan S J , Ge S S , Ding H . Knowledge-wrapping method for prediction and evaluation of material removal behavior in robotic belt grinding. Mechanical Systems and Signal Processing, 2024, 208: 110914

[28] Amersdorfer M , Meurer T . Equidistant tool path and cartesian trajectory planning for robotic machining of curved freeform surfaces. IEEE Transactions on Automation Science and Engineering, 2022, 19(4): 3311–3323

[29] Negi V S , Wang T Y , Garg H , Pullen W C , Ke X L , Kumar R R S , Choi H , Tiwari U K , Karar V , Kim D . Random adaptive tool path for zonal optics fabrication. Optics Express, 2022, 30(16): 29295–29309

[30] Zhang L , Han Y J , Fan C , Tang Y , Song X P . Polishing path planning for physically uniform overlap of polishing ribbons on freeform surface. The International Journal of Advanced Manufacturing Technology, 2017, 92(9): 4525–4541

[31] Wen Y L , Jaeger D J , Pagilla P R . Uniform coverage tool path generation for robotic surface finishing of curved surfaces. IEEE Robotics and Automation Letters, 2022, 7(2): 4931–4938

[32] Huang Z , Song R , Wan C B , Wei P X , Wang H Y . Trajectory planning of abrasive belt grinding for aero-engine blade profile. The International Journal of Advanced Manufacturing Technology, 2019, 102(1): 605–614

[33] Gao Y , Lin X J , Shi Y Y . Path generation for flexible flank polishing on freeform surfaces with uniform over-cuts. IEEE/ASME Transactions on Mechatronics, 2022, 27(3): 1796–1807

[34] Inoue T , Hirai S . Elastic model of deformable fingertip for soft-fingered manipulation. IEEE Transactions on Robotics, 2006, 22(6): 1273–1279

[35] Xiao M B , Ding Y , Fang Z J , Yang G L . Contact force modeling and analysis for robotic tilted-disc polishing of freeform workpieces. Precision Engineering, 2020, 66: 188–200

[36] Fan C , Zhao J , Zhang L , Zhou W S , Sun L N . Local material removal model considering the tool posture in deterministic polishing. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 2016, 230(15): 2660–2675

[37] Barfoot T D. State Estimation for Robotics. Cambridge: Cambridge University Press, 2017

[38] Hartley R , Trumpf J , Dai Y C , Li H D . Rotation averaging. International Journal of Computer Vision, 2013, 103(3): 267–305

[39] Wang Z W , Yan S J , Wu L , Zhang X J , Chen B J . Robust point clouds registration with point-to-point lp distance constraints in large-scale metrology. ISPRS Journal of Photogrammetry and Remote Sensing, 2022, 189: 23–35

RIGHTS & PERMISSIONS

The Author(s). This article is published with open access at link.springer.com and journal.hep.com.cn