metals

Article
Spatial Path Planning for Robotic Milling of Automotive
Casting Components Based on Optimal Machining Posture
Hao Wu 1,2 , Yudi Wang 1,2 , Xiaoxu Wei 1,2, * and Dahu Zhu 1,2

1 Hubei Key Laboratory of Advanced Technology for Automotive Components,

Wuhan University of Technology, Wuhan 430070, China; wuhao20228@[Link] (H.W.);
wyudi2022@[Link] (Y.W.); dhzhu@[Link] (D.Z.)
2 Hubei Collaborative Innovation Center for Automotive Components Technology,
Wuhan University of Technology, Wuhan 430070, China
* Correspondence: wxx2014@[Link]; Tel.: +86-27-87858200

Abstract: The robotic milling of automotive casting components can effectively reduce human
participation in the production process and enhance production efficiency and quality, but the
premise addresses the reasonable planning of machining paths. To address major challenges, this
paper proposes a spatial path planning method for the robotic milling of casting flash and burrs on
an automotive engine flywheel shell based on the optimal machining posture. Firstly, an improved
stereolithography slicing algorithm in arbitrary tangent plane direction is put forward, which solves
the problem that the existing stereolithography slicing algorithm cannot accurately extract the
contour of complex components. Secondly, the contour path curve fitting of the slicing points of
the flywheel shell is realized based on the B-spline curve. Next, a machining posture evaluation
function is established based on the robot’s stiffness performance, and the optimal machining posture
is solved and verified with simulation according to the machining posture evaluation function and
posture interpolation. Finally, the experiments indicate that the proposed method can significantly
Citation: Wu, H.; Wang, Y.; Wei, X.; enhance the machining quality, with an average allowance height of 0.33 mm, and reduce the
Zhu, D. Spatial Path Planning for machining time to 9 min, compared with the conventional manual operation, both of which satisfy
Robotic Milling of Automotive the machining requirements.
Casting Components Based on
Optimal Machining Posture. Metals Keywords: robotic milling; path planning; automotive components; slicing algorithm; robot’s
2022, 12, 1271. [Link] stiffness performance; machining posture
10.3390/met12081271

Academic Editor: Sergey N.

Grigoriev
1. Introduction
Received: 5 July 2022
Accepted: 26 July 2022
Robotic machining, as an advanced manufacturing technology [1] that conforms to
Published: 28 July 2022
the national situation [2,3], receives wide concerns in recent years due to its advantages in
improving the level of automation, as well as manufacturing efficiency and quality. This
Publisher’s Note: MDPI stays neutral
type of automated technology has been extensively used in the fields of aerospace and rail
with regard to jurisdictional claims in
transit but poses major challenges in the automotive industry. For example, removing flash
published maps and institutional affil-
and burrs on automotive casting components such as engine flywheel shells still depends
iations.
on manual operation in automotive enterprises, exhibiting the typical problems of low
efficiency, poor consistency, and high labor intensity. This situation, however, is being
changed by the robotic machining mode by virtue of flexibility, reconfigurability, and cost-
Copyright: © 2022 by the authors.
effectiveness. The premise for the robotic milling of such complex casting components lies in
Licensee MDPI, Basel, Switzerland. spatial path planning, including the machining path generation [4,5], the stereolithography
This article is an open access article (STL) slicing algorithm, and the optimal posture solution.
distributed under the terms and For the complex automotive component of flywheel shells, reasonably planning the
conditions of the Creative Commons spatial path is a primary task in robotic machining processes such as milling. Currently,
Attribution (CC BY) license (https:// many methods are put forward to generate machining paths based on complex surface
[Link]/licenses/by/ features, for instance, the section plane method [6], the parameter method [7], and the equal
4.0/). residual height method [8]. Due to the complex structure of flywheel shells, machining

Metals 2022, 12, 1271. [Link] [Link]

Metals 2022, 12, 1271 2 of 17

paths are distributed on multiple free-form surfaces. As the section plane method [9,10]
intercepts the machining surface by setting a set of equidistant section planes, it is suitable
for the task requirements of this paper.
The section plane method usually takes a triangular mesh model as the intercepting
object. STL is considered one of the most widely used file formats for triangular mesh
models. The STL slicing algorithm based on the section plane method [11] obtains the
intersection point by intersecting the section plane and the triangular patch, which can
ensure the accurate extraction of path fitting points [12,13]. However, most STL slicing
algorithms perform hierarchical slicing in the height direction and can only obtain the
section profile of a certain height, which is difficult for the generation of contour paths for
complex components. Aiming at the requirement of the contour path extraction of complex
components such as flywheel shells, this paper proposes an arbitrary tangent-plane-based
STL slicing algorithm to obtain slicing points and generates the task path by fitting the
slicing points with a B-spline curve. This method has strong versatility and can also extract
the machining path for other complex components.
In terms of robotic machining posture optimization, the majority of the existing
methods eliminate the redundant degrees of freedom in the robotic machining process by
adding constraints or optimization indicators. Based on the dexterity index, Xiao et al. [14]
and Zhu et al. [15] established a multi-objective optimization model for the elimination of
redundant degrees of freedom of machining robots and integrated the joint avoidance limit
index and the robot dexterity index to complete robot pose optimization. By virtue of the
stiffness performance index, Xiong et al. [16] proposed a new frame invariant performance
index; Chen et al. [17] used the inverse of the volume of compliance ellipsoids as the stiffness
performance index; Sun et al. [18] proposed the measured stiffness performance evaluation
metrics that explain the difference in stiffness in two opposite directions. In addition,
the constraints that comprehensively consider multiple specific indicators have also been
proposed. For example, Chen et al. [19] proposed a method to achieve posture optimization
by controlling the functional redundancy of the robot, which comprehensively considers
the deformation caused by the spindle weight and the deformation caused by the cutting
force. Then, the cutting trajectory of the robot is optimized with the kinematic performance
index as the optimization goal. For the robotic milling system in this paper, the robot should
avoid singularity and joint overrun and ensure the maximum stiffness performance.
The studies above indicate that the existing STL slicing algorithms can only obtain the
cross-sectional contour of a certain height, which cannot meet the extraction requirements
for the complex and special-shaped flywheel shell. Therefore, this paper aims to develop an
improved STL slicing algorithm for reasonable and effective robotic milling path planning.
To achieve this objective, in Section 2, an improved STL slicing algorithm in an arbitrary
direction is proposed to extract the slicing points of a flywheel shell, and the machining
curve fitting is completed through a B-spline curve from the slicing points. Meanwhile, in
Section 3, a weighted posture evaluation function is established to optimize the robotic
posture by considering robot singularity, joint avoidance limit, and stiffness performance.
Through the path point posture interpolation, the minimum value of the posture evaluation
function of each path point is obtained to realize the robotic milling posture optimization.
Finally, in Section 4, the spatial path planning of the robotic milling of the engine flywheel
shell is completed, and the superiority of the proposed method is verified with experiments.

2. STL-Based Path Planning for Robotic Milling of Flywheel Shells

2.1. Contour Path Segment Based on Geometric Features of the Flywheel Shell
STL is a file format for describing three-dimensional graphics [20], which represents
the three-dimensional graphics as several triangular patches. Before acquiring the STL file
of the flywheel shell, it is necessary to determine the contour curve segment according to its
features and task requirements. As shown in Figure 1, the flywheel shell’s flash and burrs
are mainly distributed on the upper edges, and the geometric features of the task path are
 [Link]-Based
STL-BasedPath
PathPlanning
Planningfor
forRobotic
RoboticMilling
Millingof ofFlywheel
FlywheelShells
Shells
2.1.
2.1. Contour Path Segment Based on Geometric Features of the FlywheelShell
Contour Path Segment Based on Geometric Features of the Flywheel Shell
STL
STLisisaafile
fileformat
formatfor fordescribing
describingthree-dimensional
three-dimensionalgraphics
graphics[20],
[20],which
whichrepresents
represents
the
the three-dimensional graphics as several triangular patches. Before acquiringthe
three-dimensional graphics as several triangular patches. Before acquiring theSTL
STL file
file
Metals 2022, 12, 1271 3 of 17
of
of the flywheel shell, it is necessary to determine the contour curve segment accordingto
the flywheel shell, it is necessary to determine the contour curve segment according to
its features and task requirements. As shown in Figure 1, the flywheel
its features and task requirements. As shown in Figure 1, the flywheel shell’s flash and shell’s flash and
burrs
burrsare aremainly
mainlydistributed
distributedon onthetheupper
upperedges,
edges,andandthe
thegeometric
geometricfeatures
featuresofofthe
thetask
task
characterized
path are by straight
characterized lines,
by curves,
straight and
lines, circles.
curves, Accurately
and [Link] the
Accurately
path are characterized by straight lines, curves, and circles. Accurately extracting contour
extractingcurve
the
the
of the flywheel
contour curve shell
of the is critical shell
flywheel to taskis path generation.
critical to task path generation.
contour curve of the flywheel shell is critical to task path generation.

(a)
(a) (b)
(b)

Figure
[Link]
Figure The
The automotive
automotive engine
engine
automotive flywheel
flywheel
engine shell
shell with
flywheel with
with flash/burrs
shellflash/burrs (a) and (a)
flash/burrs the and
and the
the corresponding
(a) corresponding geometric
corresponding
geometric
(b). model
geometric
model model(b).
(b).

The contour
Thecontour curve
contourcurve isisextracted
curveis extracted based
extractedbased
based on
on the
the
on section
section
the section plane
plane method,
method,
plane method,and and
and the
the task
task
the path
path
task is
path
is divided
divided into
into horizontal
horizontal plane
plane and
and non-horizontal
non-horizontal plane
plane areas
areas according
according to
to
is divided into horizontal plane and non-horizontal plane areas according to the geometric the geometric
features
featuresofofthe [Link]
thecontour. InFigure
In Figure2,2,both
Figure boththe
both thegreen
greentask
taskpaths
pathsand
andyellow
yellowcircular
circulartask
task
paths are classified into horizontal plane paths, while the red task paths are classifiedinto
paths are classified into horizontal plane paths, while the red task paths are classified into
non-horizontal
non-horizontalplane
non-horizontal planepaths
plane pathswith
paths withvertical
verticalplanes
planesandandinclined
inclinedplanes.
[Link]
canbe
can beseen
seenthat
that
the
thered
redtask
taskpaths
pathsdivide
dividethe thegreen
greenones
onesinto
intoseven [Link]
sevensections. totalof
of10
10horizontal
horizontalplane
plane
paths exist with different heights, and their actual positions are shown
paths exist with different heights, and their actual positions are shown inFigure
heights, and their actual positions are shown in
in Figure
Figure3.3.3.

Figure
Figure2. The
Thetask
taskpath
pathdistributions
distributionsof the
theflywheel shell.
Figure [Link] task path distributions ofofthe flywheelshell.
flywheel shell.

Table 1 further shows the plane heights of these 10 horizontal plane paths. According
to the height values, the plane equation can be obtained, so the section plane of the task
paths on the horizontal plane is given by z = −h.

Table 1. The path heights on the horizontal plane.

Path Segment Height h (mm) Path Segment Height h (mm)

1 42.24 6 42.19
2 11.47 7 53.00
3 42.79 8 15.02
4 3.99 9 18.81
5 13.94 10 15.48
Metals 2022, 12, x FOR PEER REVIEW 4 o
Metals 2022, 12, 1271 4 of 17

Figure 3. The actual positions of the task path segments on the horizontal plane.

Table 1 further shows the plane heights of these 10 horizontal plane pa

to the height values, the plane equation can be obtained, so the section pla
paths on the horizontal plane is given by z = −h.

Table 1. The path heights on the horizontal plane.

Path Segment Height h (mm) Path Segment Heig

1 42.24 6
2 11.47 7
3 42.79 8
4 3.99 9
5 13.94 10

Figure 3.
By
Figure 3. The actual
contrast,
The actual positions
for of
positions theofnon-horizontal
the the path
task task segments
path segments
onplane,on the
thehorizontal
the horizontalplane plane.
plane. equation needs to
by Bythree points
contrast, that
for the are not collinear.
non-horizontal Theequation
coordinates of determined
the three poin
Table 1 further shows the planeplane, the plane
heights of these needs to be
10 horizontal plane paths. Accord
collinear
by three pointscanthat
bearetaken from the
not collinear. Thegeometric
coordinatesmodel andpoints
of the three correspond to the ac
that are not
to the height values, the plane equation can be obtained, so the section plane of the t
collinear can be taken
For aonstraight-line from the
segment, geometric model
the two and correspond to the actual positions.
paths the horizontal
For a straight-line segment,plane
the twoisends
given z ends
by line
of the = −h. of the line segment and the cent
segment and the center of the circle are
are selected
selected to constructtotheconstruct
section plane, the section
as shown plane,
in Figure 4. Theas showncoordinates
three-point in Figure of 4. Th
coordinates
Table
each 1. The
red taskpath ofheights
path each on
segment red
arethetask path
Tablesegment
horizontal
listed in plane.
2. are listed in Table 2.
Path Segment Height h (mm) Path Segment Height h (mm)
1 42.24 6 42.19
2 11.47 7 53.00
3 42.79 8 15.02
4 3.99 9 18.81
5 13.94 10 15.48

By contrast, for the non-horizontal plane, the plane equation needs to be determi
by three points that are not collinear. The coordinates of the three points that are
collinear can be taken from the geometric model and correspond to the actual positio
For a straight-line segment, the two ends of the line segment and the center of the ci
are selected to construct the section plane, as shown in Figure 4. The three-p
coordinates of each red task path segment are listed in Table 2.
Figure
Figure 4. The
4. The section
section plane constructed
plane constructed by the
by the three points three
that points
are not that are not collinear.
collinear.
 Table 2. The three-point coordinates of non-horizontal path segments.

Path
P1 P2 P3
Metals 2022, 12, 1271 Segment 5 of 17
11 [285.43 232.34 −42.24] [286.66 236.66 −11.47] [0 0 0]
12 [−90.04 291.47 −11.47] [−90.00 292.77 −36.23] [0 0 0]
Table132. The three-point
[−90.00 292.77 −36.23]
coordinates [−99.32 290.04
of non-horizontal path−37.86]
segments. [−98.01 261.15 −38.10]
14 [−98.01 261.15 −38.10] [−90.09 248.40 −49.44] [−90.50 236.58 −56.69]
Path P P2 P3
15
Segment [−90.50 236.58 1−56.69] [−126.40 248.00 −76.57] [−176.21 265.52 −42.79]
16 11 232.34 −
[−211.39 263.10 −42.79] [−207.61 261.80 −101.79]
[285.43 42.24] [286.66 236.66 − 11.47] [0 0 0]
[0 0 0]
12 [−90.04 291.47 −11.47] [−90.00 292.77 −36.23] [0 0 0]
17 13 [−331.62 40.66
[−90.00 −42.79]
292.77 −36.23] [−332.55 39.90−37.86]
[−99.32 290.04 −3.99] [−98.01 [0 0 0]−38.10]
261.15
14 [−98.01 261.15 −38.10] [−90.09 248.40 −49.44] [−90.50 236.58 −56.69]
18 15 [−176.69 −176.69 −3.99]
[−90.50 236.58 −56.69] [−183.14 −183.14 −68.25]
[−126.40 248.00 −76.57] [−173.85 −173.84 −121.67]
[−176.21 265.52 −42.79]
16 [−211.39 263.10 −42.79] [−207.61 261.80 −101.79] [0 0 0]
19 17 [176.44[−−176.44 −13.94] [183.11
331.62 40.66 −42.79]
−183.11 −67.87] [173.87 −173.63
[−332.55 39.90 −3.99] [0 0 0]
−123.11]
20 18 [−176.69
[204.84 −120.39−176.69 −3.99] [205.16
−13.94] [−183.14−121.93
−183.14 − 68.25]
−34.06] [−173.85 −[0 0 0]−121.67]
173.84
19 [176.44 −176.44 −13.94] [183.11 −183.11 −67.87] [173.87 −173.63 −123.11]
21 20 [270.55 −36.74
[204.84 −42.19]
−120.39 −13.94] [272.46
[205.16 − −34.26
121.93 −−51.19]
34.06] [245.06 −24.94
[0 0 0] −111.18]
21 [270.55 −36.74 −42.19] [272.46 −34.26 −51.19] [245.06 −24.94 −111.18]
22 22 [313.87 19.71
[313.87 −52.98]
19.71 −52.98] [327.15
[327.1543.81
43.81 −−42.24]
42.24] [0
[0 000]0]

Given
Given the the coordinates of three
coordinates of three points
pointsthat
thatare
arenot
notcollinear,
collinear,PP 1(x1, y1, z1), P2(x2, y2, z2),
1 (x1 , y1 , z1 ), P2 (x2 , y2 , z2 ),
P
P3 (z3 , y3 , z3 ), the plane equation passing through the three points isgiven
3 (z3 , y 3 , z3 ), the plane equation passing through the three points is givenby: by:
AX + BY + CZ + D = 0 (1)
AX + BY + CZ + D = 0 (1)
where the coefficients A, B, C, and D of the plane equation are represented by the
where the coefficients
coordinates of the threeA, B, C, and
points P1,DPof the plane
2, and equation are represented by the coordinates
P3, respectively:
of the three points P1 , P2 , and P3 , respectively:
1 y1 z1 x1 1 z1 x1 y1 1 x1 y1 z1
A1 = y11 yz2 1 z2 , B =x1x2 1 1 z1z2 , C = x21 y21 11 , D = x2 x1 y 2 y1 z2 z1
A = 1 y12 yz3 2 ,zB3 = x2x3 1 1 z2z3, C = x32 y32 11 , D =x3 x2 y 3 y2 z3 z2
1 y3 z3 x3 1 z3 x3 y3 1 x3 y3 z3
In addition, due to the small number of triangular patches in the STL file exported
by 3D In modeling
addition, duesoftware
to thesuch
smallas SolidWorks
number (Solidworks2020,
of triangular SolidWorks,
patches in the Concord,
STL file exported by
MA, USA), the number of fitting points needs to be increased to
3D modeling software such as SolidWorks (Solidworks2020, SolidWorks, Concord, MA, ensure the accuracy of
the path
USA), thepoint fitting
number ofcurve.
fittingTherefore,
points needs meshto subdivision
be increasedis to
performed
ensure the to accuracy
ensure thatof the
the
mesh density
path point is incurve.
fitting the same order of
Therefore, meshmagnitude
subdivisionas the minimum distance
is performed to ensureofthat1 mm from
the mesh
the pathispoint
density in theinsame
robot offline
order programming
of magnitude as thesoftware
minimum such as RobotStudio.
distance of 1 mm from Finally, the
the path
pointfile
STL in of
robot offline programming
the flywheel shell model is software
[Link]
RobotStudio.
subdivisionFinally, the STL
is shown file of
in Figure
the flywheel shell model is exported. The mesh subdivision is shown in Figure 5.
5.

(a) (b)

Figure
Figure 5.
5. Comparison
Comparison of
of before
before (a)
(a) and
and after (b) mesh subdivision of the flywheel shell STL file.

2.2. STL Slicing Algorithm in an Arbitrary Direction Based on a Complex Parting Surface
The STL slicing algorithm is a commonly used printing method in the field of 3D
printing. However, most of the current STL slicing algorithms perform layered slicing in
the height direction and can only obtain the cross-sectional contour of a certain height.
Thus, it is difficult to generate the contour path of complex components such as flywheel
shells. Therefore, the STL slicing algorithm needs improvement. In this paper, the method
 2.2. STL Slicing Algorithm in an Arbitrary Direction Based on a Complex Parting Surface
The STL slicing algorithm is a commonly used printing method in the field of 3D
printing. However, most of the current STL slicing algorithms perform layered slicing in
the height direction and can only obtain the cross-sectional contour of a certain height.
Metals 2022, 12, 1271 Thus, it is difficult to generate the contour path of complex components such as flywheel 6 of 17
shells. Therefore, the STL slicing algorithm needs improvement. In this paper, the method
of generating the tangent plane in an arbitrary direction is proposed accordingly, and the
constraint of eliminating
of generating the tangent the planeredundant inner direction
in an arbitrary and outer contours accordingly,
is proposed is established.
andThe
detailed steps are as follows:
the constraint of eliminating the redundant inner and outer contours is established. The
Stepsteps
detailed 1. Determine the tangent plane equation in an arbitrary direction.
are as follows:
The
Steptangent planethe
1. Determine of the taskplane
tangent path equation
of the flywheel shell has
in an arbitrary great randomness, and
direction.
The tangent
the tangent plane plane
equationof the taskbepath
must of the flywheel
determined first. It shell has great
is known that randomness,
the three pointsandthat
the tangent plane equation must be determined first. It is known that
are not collinear can determine a plane, and the three points that are not collinear are the three points
that are notselected
arbitrarily collinearincan determine
each a plane,
path of the and the
flywheel three
shell. Forpoints that are not collinear
the straight-line segmentare path,
arbitrarily selected in each path of the flywheel shell. For the straight-line
the two endpoints of the path and the center of the flywheel housing circle are selected, segment path,
the the
and twosection
endpoints of the
plane path and
equation the centeras
is expressed of the flywheel housing circle are selected,
and the section plane equation is expressed as
( Z = h horizontal plane
Z=
AX + BY + CZ + D = 0 non-horizontal
h horizontal plane plane (2)
 (2)
AX + BY + CZ + D = 0 non-horizontal plane
Step 2. Determine whether the triangular patch intersects with the tangent plane.
In [Link]
Step Determine plane equation,
whether the normalpatch
the triangular vector N(A, B,with
intersects C), anythe point
tangent onplane.
the tangent
plane,Inand the tangent
the three plane equation,
vertices of thethetriangular
normal vector
patchN(A, B, C), any
constitute pointvectors.
three on the tangent
The inner
plane, and the three vertices of the triangular patch constitute
product of the three vectors with the normal vector N of the tangent plane can obtain three vectors. The inner three
product of the three vectors with the normal vector N
values. If a value exists that equals 0 among the three values, it means that the three
of the tangent plane can obtain tangent
values.
plane If a value
passes through exists
thethat equals
vertices in 0the
among the three
triangular values,
patch, it means
and no that the
intersection tangent
occurs. If the
plane passes through the vertices in the triangular patch, and no intersection occurs. If the
three values are all positive or all negative, it means that the three vertices of the triangular
three values are all positive or all negative, it means that the three vertices of the triangular
patch are on the same side of the plane; that is, the triangular patch does not intersect with
patch are on the same side of the plane; that is, the triangular patch does not intersect with
the tangent plane. Even if the three values are not all positive or all negative, it means that
the tangent plane. Even if the three values are not all positive or all negative, it means that
the
thethree three vertices
vertices of the triangular
of the triangularpatchpatcharearelocated
locatedonon both
both sides
sides of of
thethe plane;
plane; thatthat
is, theis, the
triangular patch intersects with
triangular patch intersects with the tangent [Link] tangent plane.
As
As shown
shown in in Figure
Figure6,6,thethethree
three vertices
vertices of the
of the triangular
triangular patchpatch are defined
are defined as P4 (xas4 , P
y44 ,(x4,
yz4,4 ),z4P
),5P(x
5 (x 5, y5, z5), and P6 (x6, y6, z6), and the three vertices constitute the vectors P1P4, P1P5,
5 , y5 , z5 ), and P6 (x6 , y6 , z6 ), and the three vertices constitute the vectors P1 P4 , P1 P5 ,
and
and P1 P6 withthe
P 1P 6 with thepoint
point PP11 on
on the tangent plane;thus,
tangent plane; thus,thethejudgment
judgment formula
formula is written
is written as as
−
→( P1 P4 N )  0;
( P1 P4 · N ) ∧ 0;
−
→( P P N )  0; (3)
( P1 P51· N5 ) ∧ 0; (3)
−
→( P P6 N )  0
( P1 P61 · N )∧0

Figure 6. Intersection of the triangular patch with the section plane.

Figure 6. Intersection of the triangular patch with the section plane.
Step 3. Find the intersection point of the triangular patch with the tangent plane.
According to the distribution of the vertices of the triangular patch on the tangent plane
in the previous step, the triangular patch that intersects with the tangent plane must have
three vertices distributed on both sides of the tangent plane. This means that one vertex
is located on one side of the tangent plane, and the other two vertices are located on the
 Metals 2022, 12, 1271 7 of 17

other side of the tangent plane. By drawing a straight line from one vertex to the other two
vertices, the straight-line equation is established according to the vertex coordinates and
then is combined with the tangent plane equation to obtain the intersection coordinates.
Assuming that P1 P4 and P1 P5 are greater than zero, P1 P6 is less than zero, and the
intersections of the triangular patch with the section plane are P4 P6 , P5 P6 , then the straight-
line equation is shown below:
X − X4 Y −Y4 Z − Z4
X6 − X4 = Y6 −Y4 = Z6 − Z4
X − X5 Y −Y5 Z − Z5 (4)
X6 − X5 = Y6 −Y5 = Z6 − Z5

Combined with the section plane Equation (2), the intersection points are acquired,
respectively.
Step 4. Eliminate the inner contours.
By traversing all the triangular patches, the obtained intersections constitute the cross-
sectional contour between the tangent plane and the flywheel shell, which includes the task
path segments, as well as the possible useless inner and outer contours. Therefore, in order
to eliminate the inner contours, two measures are taken. Firstly, the flywheel shell model is
filled as a solid body to eliminate the inner contour, and secondly, for the contours outside
the scope of the task path segment, the desired slicing points are filtered out according to
the coordinates of the starting and ending points of the task segment, as well as the given
definition domain of X and Y. In this way, the cross-sectional contour of the specified area
model and the tangent plane can be solved.

2.3. Path Fitting Based on the B-Spline Curve Slicing Point

By using the STL slicing algorithm, the slicing points of the flywheel shell contour
can be acquired, and the path curve is fitted accordingly. As the fitting curve needs
to pass through all the slicing points, it is suitable for the spline curve in the field of
computational geometry. The common splines include Bezier spline, B spline, and NURBS
spline. Compared with the simple Bezier spline and complicated NURBS spline, the B
spline can meet the fitting requirements and is, therefore, selected to fit the task path curve
in this paper.
Before the curve fitting, the slicing points need to be sorted. Since the set of path points
of each segment is ordered, it is only necessary to consider the connection order between seg-
ments. The path sequence is finally determined as 1→11→2→12→13→14→15→3→17→4
→5→20→6→7→22. The rest of the paths have no order requirement and can be fitted
individually. In addition, denser slicing points are obtained through the above mesh subdi-
Metals 2022, 12, x FOR PEER REVIEW 8 of 1
vision. As shown in Figure 7, the slicing point spacing is significantly reduced by increasing
the initial 1400 slicing points to 12,000.

[Link]
Figure Comparison of the
of the slicing
slicing point
point sets before
sets before (a)after
(a) and and(b)
after (b) preprocessing
preprocessing of the flywhee
of the flywheel
shell STL file.
shell STL file.

After completing the above operation, the 15 segment paths are fitted in sequence
The equation of the B-spline curve is written as follows [21]:
n
p(u) =  di Ni , k (u) (5
i =0

Note that the B-spline curve is composed of the cumulative sum of multiple B-splin
Metals 2022, 12, 1271 8 of 17

After completing the above operation, the 15 segment paths are fitted in sequence.
The equation of the B-spline curve is written as follows [21]:
n
p(u) = ∑ di Ni,k (u) (5)
i =0

Note that the B-spline curve is composed of the cumulative sum of multiple B-spline
basis functions Ni,k (u), and the shape of the curve is controlled by the control point di .
The second subscript k of the B-spline basis function Ni,k (u) represents the degree of the
curve, and i is the serial number from 0 to n. The argument u of the basis function is the
parameter in the node vector U, and the parameters in the node vector U are arranged in
non-decreasing order, from 0 to n + k + 1. The fitting principle is as follows:
Step 1. Determine the node vector.
Assuming that m + 1 data points qi (I = 0, 1, . . . , m) are obtained, a cubic general
non-uniform B-spline curve should be constructed. First, it should be determined whether
the curve to be constructed is an open curve or a closed curve; that is, whether the data
points at the beginning and end are the same. If they are, it is a closed curve; otherwise, it
is an open curve. For a closed curve, if C2 is not required to be continuous at the overlap,
it is equivalent to an open curve. The data point at the closure is used as the beginning
and ending point of the curve, and the other data points are used as intermediate segment
points so that the curve has a total of m segments. The number of spline control points is
equal to the number of data points plus the degree of the curve minus one. Assuming that
there are n + 1 control points dj (j = 0, 1, . . . , n), it has n = m + 2. Since the degree of the
curve is 3 times, and the repeatability of the domain endpoint is 4, the standard domain
u ∈ [u3 ,un + 1 ] = [0,1]. Thus, the endpoints with a repeatability of 4 at the beginning and
ending are given by
u0 = u1 = u2 = u3 = 0
u n +1 = u n +2 = u n +3 = u n +4 = 1
For the standard definition domain of m + 1 data point qi , the normative accumulation
chord length method is used to calculate the chord length of m + 1 data points connected in
turn. Assuming that the m-th chord length is lm , the total chord length L is expressed as
m
L= ∑ lm (6)
1

According to the parameter sequence ui (i = 0, 1, . . . , m) of the node vector, it has

i

 ∑ li
1
_ui = L , i = 1, 2, . . . , m (7)
_u0 = 0


The correspondence between the standard domain and the node vector is u3 + i = ui
(i = 0, 1, . . . , m). For the C2 continuous, closed curve in this paper, the standard definition
domain is the same as above. The endpoints at both ends of the domain are set to

u0 = un−2 − 1, u1 = un−1 − 1, u2 = un − 1, un+2 = 1 + u4 , un+3 = 1 + u5 , un+4 = 1 + u6

Step 2. Inversely calculate the control vertex.

The fitted curve equation can be written as

n i
p(u) = ∑ d j Nj,3 (u) = ∑ d j Nj,3 (u), u ∈ [ui , ui+1 ] ⊂ [u3 , un+1 ] (8)
j =0 j = i −3
Metals 2022, 12, 1271 9 of 17

By substituting the node vector parameter u ∈ [ui , ui+1 ] ⊂ [u3 , un+1 ] of the standard
definition domain into the node vector determined by Step 1, the following interpolation
conditions need to be satisfied:

i
 p(ui ) = ∑ d j Nj,3 (ui ) = qi−3 , i = 3, 4, . . . , n



j = i −3
n (9)
p ( u ) = ∑ d N ( u ) = q

n +1 j j,3 n+3 m



j = n −3

Note that Equation (9) contains n−1 equations in total. Since q0 = qm , one equation is
repeated; thus, there are n−2 after subtraction. Based on the definition of a closed curve,
the three control points at the beginning and end should be the same, that is,

d n −2 = d 0 , d n −1 = d 1 , d n = d 2

The required n + 1 control points are reduced to n−2, and the number of equations is
equal to the number of unknowns, which is solved and represented by a matrix as follows:

N1,3 (u3 ) N2,3 (u3 ) N0,3 (u3 ) d1 q0

    
 N1,3 (u4 ) N2,3 (u4 ) N3,3 (u4 )    d2  q1
.. ..
    
 .. .. .. 
=
  
(10)

 . . . 
   . 
 .
 Nn−4,3 (un ) Nn−3,3 (un ) Nn−2,3 (un )   d n −3   q n −4 
Nn−1,3 (un+1 ) Nn−3,3 (un+1 ) Nn−2,3 (un+1 ) dn−2 q n −3

The matrix equation is further rewritten as

b1 c1 a1 d1 e1
    
 a2 b2 c2  d2
   e2
.. ..
    
 .. .. .. 
=
  
(11)

 . . . 
  . 
 .
 a n −3 bn − 3 c n −3 d n −3
    e n −3 
c n −2 a n −2 bn − 2 d n − 2 e n −2

where
( ∆ i +2)2
ai = ∆ i + ∆ i +1 + ∆ i +2
∆ i +2 ( ∆ i + ∆ i +1 ) ∆ i +1 ( ∆ i +2 + ∆ i +3 )
bi = ∆ i + ∆ i +1 + ∆ i +2 + ∆ i +1 + ∆ i +2 + ∆ i +3
( ∆ i +1 )2
ci = ∆ i +1 + ∆ i +2 + ∆ i +3
ei = ( ∆ i +1 + ∆ i +2 ) q i −1
∆ i = u i +1 − u i
i = 1, 2, . . . , n − 2
The control points can be acquired by solving the linear equations. By substituting the
control points of the task path into Equation (5), the fitting curve is obtained as shown in
Figure 8.
Based on these steps, the independent path segments are drawn in turn to obtain
the complete task path. Finally, the complete task path is imported into the robot offline
programming software (ABB RobotStudio) and matched with the flywheel shell model. As
shown in Figure 9, the two fit well. Specifically, the distance values between the centers of
the three circles on the trajectory are measured as 98.61 mm, 238.33 mm, and 185.75 mm,
and then compared with the model, indicating the accuracy of the task path.
 Metals 2022, 12, x FOR PEER REVIEW 10 of 17

The control points can be acquired by solving the linear equations. By substituting
Metals 2022, 12, 1271 10 of 17
the control points of the task path into Equation (5), the fitting curve is obtained as shown
in Figure 8.

Figure 8. The fitting curve of the flywheel shell task path.

Based on these steps, the independent path segments are drawn in turn to obtain the
complete task path. Finally, the complete task path is imported into the robot offline
programming software (ABB RobotStudio) and matched with the flywheel shell model.
As shown in Figure 9, the two fit well. Specifically, the distance values between the centers
of the three circles on the trajectory are measured as 98.61 mm, 238.33 mm, and 185.75
[Link]
Thefitting
fittingcurve
curveof
ofthe
the flywheelshell
shelltask
taskpath.
path.
Figure
mm, and then compared withflywheel
the model, indicating the accuracy of the task path.
Based on these steps, the independent path segments are drawn in turn to obtain the
(a)
complete task path. Finally, the complete
(b) task path is imported into the robot offline
programming software (ABB RobotStudio) and matched with the flywheel shell model.
As shown in Figure 9, the two fit well. Specifically, the distance values between the centers
of the three circles on the trajectory are measured as 98.61 mm, 238.33 mm, and 185.75
mm, and then compared with the model, indicating the accuracy of the task path.

(a) (b)
Figure 9.
Figure 9. Comparison of the
Comparison of the task
task path
path (a)
(a) with
with the
the standard
standard model
model (b).
(b).

3. Optimal Machining
3. Optimal Machining Posture
Posture Solution
Solution Based
Based on
on the
the Robot’s
Robot’s Stiffness
Stiffness Performance
Performance
3.1. Machining Posture Evaluation Function Based on the Robot’s Stiffness Performance
3.1. Machining Posture Evaluation Function Based on the Robot’s Stiffness Performance
The
The robotic
robotic machining
machining pathpath includes
includes the
the position
position and and orientation
orientation (pose)
(pose) information
information
of
of the path point, in which the pose information is the frame at the path point.
the path point, in which the pose information is the frame at the path point. In
In robotic
robotic
milling, the tool is vertically downward by default, and the frame can only be rotated
milling,
Figure [Link] tool is vertically
Comparison of the taskdownward
path (a) withbythedefault,
standardand modelthe(b).
frame can only be rotated
around the tool axis to adjust the posture. Therefore, aiming at the redundancy of the
around the tool axis to adjust the posture. Therefore, aiming at the redundancy of the
robot’s degrees of freedom, a path posture optimization function is proposed to solve the
robot’s degrees
3. Optimal of freedom,
Machining a path
Posture postureBased
Solution optimization functionStiffness
on the Robot’s is proposed to solve the
Performance
optimal posture of the endpoints of the path segment. Then, the intermediate points of the
optimal posture
3.1. Machining of the endpoints
PosturetoEvaluation of the path Based on the Robot’s Stiffness Performance of the
segment. Then, the intermediate points
path are interpolated complete Function
the path planning.
path are interpolated to complete the path planning.
The to
Due robotic machining
the large pathand
tangential includes
normal theforces
position androbotic
in the orientation (pose)
milling information
process, better
Due to the large tangential and normal forces in the robotic milling process, better
stiffness performance
of the path is required
point, in which in these
the pose two directions.
information Therefore,
is the frame thepath
at the directional
point. Instiffness
robotic
stiffness performance is required in these two directions. Therefore, the directional
index K the
milling, is used
toolasisfollows [22]:
vertically downward by default, and the frame can only be rotated
stiffnessefindex Kef is used as follows [22]:
around the tool axis to adjust the posture. Therefore, 1 aiming at the redundancy of the
Ke f = T (12)
e f C f1d e f
robot’s degrees of freedom, a path posture optimization function is proposed to solve the
Kef = T (12)
optimal posture of the endpoints of the path segment. e C e f Then, the intermediate points of the
where f is the magnitude of the force vectorf atfdthe end, and ef is the direction of the
path are interpolated to complete the path planning.
force
wherevector.
f is the magnitude of the force
Due to the large tangential and vector
normalatforces
the end, androbotic
in the ef is themilling
direction of the better
process, force
Set the directional stiffness index KE as
vector.
stiffness performance is required in these two directions. Therefore, the directional
stiffness index Kef is used as follows [22]: 1
KE = q (13)
2 +1K 2
KX
Kef = Y
(12)
eTf C fd e f
where KX is the normal vector and KY is the tangential vector, and the directional stiffness
whereisfobtained
index is the magnitude
from the of the force vector
compliance [Link] the end, and ef is the direction of the force
vector.     
d C f d Cmd f
= (14)
δ C f δ Cmδ m
Metals 2022, 12, 1271 11 of 17

 
c11 c12 c13
Cfd = c21
 c22 c23  (15)
c31 c32 c33
The directions of each coordinate axis at the end of the robot are [1,0,0], [0,1,0], and
[0,0,1], respectively, and the stiffness in the direction of the end coordinate axis is obtained.

1 1 1
KX = ;K = ; KZ = (16)
c11 Y c22 c33

The directional stiffness index KE is normalized, and the stiffness posture evaluation
function Erigidity is established, as shown in Equation (17).

KE − KEmin
Erigidity (θ ) = (17)
KEmax − KEmin

An accessibility index is proposed for judging which axis configuration is the closest to
the mechanical origin of the robot when multiple solutions and multiple axis configurations
exist for a certain position. Therefore, the accessibility posture evaluation function Ereach is
defined as shown in Equation (18).

6
|θi −θ0i | θi −θimin
∑ ; 0< < 0.9

Ereach (θ ) = (|θi |max −θ0i ) θimax −θimin (18)
i =1
∞; others


where θ 0i is the i-th joint angle when the robot is at zero position, θ i max is the maximum
value among all solutions of the i-th joint, and θ imin is the minimum value among all
solutions of the i-th joint. This formula is finally normalized.
For the singularity index, there are three common singularities in a six-axis serial
robot, namely the wrist joint singularity, the shoulder joint singularity, and the elbow joint
singularity. Since the overall posture of the robot in the robotic milling experiment of the
flywheel shell is downward, both the shoulder joint and elbow joint singularities are not
considered. For the wrist joint singularity, changing the second and third axes can make
the robot reach most of the positions. Note that the joint angle 5 may be zero during milling
operation, which needs to be avoided.
The range of robot joint 5 is −135◦ ~135◦ , and 90% of the range is taken as effective. The
absolute value of the joint 5 should be better than 120◦ . The singularity posture evaluation
function Esingular is shown in Equation (19).

||θ5 | − 120|
Esingular (θ ) = (19)
120
The above posture evaluation functions affect the robotic milling performance, among
which the robot’s stiffness has the greatest impact on the stability of the milling process.
According to the investigation by Huo and Baron [23], the weight of singularity is greater
than that of accessibility. After several attempts in this study, the rigidity weight is set to
0.5, followed by singularity to 0.3 and accessibility to 0.2. The overall posture evaluation
function is finally rewritten as follows:
n o
E = min 0.5∆Erigidity + 0.3∆Esingular + 0.2∆Ereach (20)

Note that the smaller E value indicates a better machining posture.

3.2. Optimal Machining Posture Solution and Simulation Verification

In the process of robotic milling of the flywheel shell, the direction of the milling cutter
is always vertically downward, and its tool coordinate system can be arbitrarily rotated
 Note that the smaller E value indicates a better machining posture.

3.2. Optimal Machining Posture Solution and Simulation Verification

In the process of robotic milling of the flywheel shell, the direction of the millin
Metals 2022, 12, 1271 cutter is always vertically downward, and its tool coordinate system can12be arbitraril
of 17
rotated around the Z axis. Therefore, from the coordinates of the endpoints of the pat
segment, an infinite number of joint angle combinations can be obtained.
aroundWhenthe Z axis. Therefore,
the posture from the coordinates
information of the two of endpoints
the endpoints
of of the path
a curve is segment,
known, the an postur
infinite number of joint angle combinations can be obtained.
of the intermediate point can be obtained by interpolating the posture of the tw
When the
endpoints. Forposture
the two information
endpoints ofMtheand
twoNendpoints of a curve
of any curve, is known,
the initial frametheofposture
the two point
of the intermediate point can be obtained by interpolating the posture of the two endpoints.
is firstly set, and the initial posture is taken into the posture evaluation function. Then, th
For the two endpoints M and N of any curve, the initial frame of the two points is firstly set,
initial frame is rotated around the Z axis. In this way, the evaluation value of all the frame
and the initial posture is taken into the posture evaluation function. Then, the initial frame
isaround
rotated the Z axis
around isaxis.
the Z calculated. Finally,
In this way, the frame
the evaluation withofthe
value smallest
all the framesevaluation
around the value i
selected as the optimal posture, and the rotation angles Ang and
Z axis is calculated. Finally, the frame with the smallest evaluation value is selected
M Ang N are as
recorded.
the
optimal Toposture,
interpolate
and thetherotation
posture of the
angles Angintermediate
M and AngN are path points, the interpolation poin
recorded.
angleTo is
interpolate
determined the byposture of the intermediate
the distance between the path points,
path theThe
points. interpolation
total lengthpoint
of the pat
angle is determined by the distance between the path points. The total
segment MN is set as the sum of the path segments connected by the intermediate pat length of the path
segment
points andMNis is denoted
set as theas sum of theThe
Length. pathdistance
segments connected by the
corresponding intermediate
to each path point i
intermediate
points and is denoted as Length. The distance corresponding to each intermediate point
the sum of all previous line segments and is denoted as Si. The i-th intermediate poin
is the sum of all previous line segments and is denoted as Si . The i-th intermediate point
angle Angi between the path segment MN is given by
angle Ang between the path segment MN is given by
i
Si
Sii =
Ang  ( AngN − AngM ) + AngM (21
Angi = × ( Ang N − Ang
Length M ) + Ang M (21)
Length
For the known coordinate frames of all path points, only the endpoints are verified
For the known coordinate frames of all path points, only the endpoints are verified.
In order to avoid the problem with intermediate points, the posture of all intermediat
In order to avoid the problem with intermediate points, the posture of all intermediate
pointsis is
points inversely
inversely solved
solved to check
to check whetherwhether
they arethey are reachable,
reachable, whether
whether collisions collisions
occur, and occu
whether they are singular. If the above problems exist, the point is set as the endpoint, theendpoin
and whether they are singular. If the above problems exist, the point is set as the
the curve
curve is divided
is divided into twointo two and
sections, sections, andsteps
the above the are
above steps Based
repeated. are repeated. Based on th
on the optimal
optimal
posture posture
solution, thesolution,
feasibilitythe feasibility
of the algorithmof the algorithm
is verified is verified
in the simulation in theThe
software. simulatio
simulation
software. results in Figure 10results
The simulation show that the robot
in Figure 10can run the
show thatpath
thesegments
robot can without
run the pat
collision
segments according
withouttocollision
the obtained optimalto
according posture paths. optimal posture paths.
the obtained

Figure10.10.
Figure Simulation
Simulation of the
of the flywheel
flywheel shell shell path segments.
path segments.

[Link]
Experimental Verification
Verification
4.1. Experimental Setup
4.1. Experimental Setup
Figure 11 shows the experimental platform for the robotic milling of the automotive
engineFigure 11 shell
flywheel shows theand
flash experimental platform forequipment
burrs. The experimental the roboticmainly
milling of the automotiv
included an
engine flywheel shell flash and burrs. The experimental equipment mainly
industrial robot (ABB IRB6700-200/2.60) and an automatic tool-changing electric spindle included a
(NR4040-AQC).
industrial robot The workpiece
(ABB to be machined
IRB6700-200/2.60) andwas
ana automatic
flywheel shell with dimensions
tool-changing of spindl
electric
600 mm × 560 mm × 156 mm in length, width, and height, respectively. It was
(NR4040-AQC). The workpiece to be machined was a flywheel shell with dimensions o made of
aluminum alloy 50021, and the corresponding chemical composition is listed in Table 3. The
flywheel shell flash/burr to be removed was 0~5 mm in width and 1~2 mm in thickness.
A four-edge tungsten steel milling cutter with a diameter of 10 mm, a length of 75 mm,
a helix angle of 35◦ , and a maximum cutting hardness of HRC48 was used. The spindle
Metals 2022, 12, x FOR PEER REVIEW 13 of 17

600 mm × 560 mm × 156 mm in length, width, and height, respectively. It was made of
aluminum alloy 50021, and the corresponding chemical composition is listed in Table 3.
Metals 2022, 12, 1271 13 of 17
The flywheel shell flash/burr to be removed was 0~5 mm in width and 1~2 mm in
thickness. A four-edge tungsten steel milling cutter with a diameter of 10 mm, a length of
75 mm, a helix angle of 35°, and a maximum cutting hardness of HRC48 was used. The
speed and
spindle the robot’s
speed and thefeed speed
robot’s were
feed set were
speed as 8000
setrpm and rpm
as 8000 5 mm/s,
and 5respectively, and the
mm/s, respectively,
transition speed was set as 50 mm/s for efficiency.
and the transition speed was set as 50 mm/s for efficiency.

Electric spindle with

automatic tool changing

Milling
cutter

Industrial
Robot

Flywheel shell

Figure11.
Figure [Link]
Experimentalplatform
platform for
forrobotic
roboticmilling
millingof
ofthe
theautomotive
automotiveengine
engineflywheel
flywheelshell
shellflash
flash
andburrs.
and burrs.

[Link]
Table Chemicalcomposition
compositionofofthe
thealuminum
aluminumalloy
alloy50021.
50021.

Aluminum
Aluminum Ingredient
IngredientContent (%)
Content (%)
Alloy
Alloy Si Si FeFe CuCu ZnZn Mg
Mg Mn
Mn Pb
Pb Sn
Sn Ti + Zr
Ti + Zr
50021 4.0–6.0 ≤1 2.0–40 ≤1 ≤0.5 ≤0.5 ≤1 ≤1 ≤1
50021 4.0–6.0 ≤1 2.0–40 ≤1 ≤0.5 ≤0.5 ≤1 ≤1 ≤1

As
Asshown
showninin
Figure 12, 12,
Figure the path segments
the path for planning
segments the optimal
for planning machining
the optimal posture
machining
are the plane segment paths in red, which accounts for 90% of the overall
posture are the plane segment paths in red, which accounts for 90% of the overallmachining path.
Thus, the red segment paths were used for experimental verification. Before performing
machining path. Thus, the red segment paths were used for experimental verification. the
robotic milling operation, tool calibration and workpiece calibration had
Before performing the robotic milling operation, tool calibration and workpiece to be completed.
The tool coordinate system was calibrated with the four-point method, while the workpiece
calibration had to be completed. The tool coordinate system was calibrated with the four-
coordinate system was calibrated with the three-point method [24]. Additionally, in order
point method, while the workpiece coordinate system was calibrated with the three-point
to ensure machining safety, it is necessary to verify the robotic milling path of the flywheel
method [24]. Additionally, in order to ensure machining safety, it is necessary to verify
shell. The entire machining path height was enhanced by 3 mm with a slow running
Metals 2022, 12, x FOR PEER REVIEW speed
14 of 17
the robotic milling path of the flywheel shell. The entire machining path height was
to observe the collision between the milling cutter and the flywheel shell.
enhanced by 3 mm with a slow running speed to observe the collision between the milling
cutter and the flywheel shell.

Figure 12.
Figure 12. Robotic
Robotic milling
millingpath
pathofofthe
theflywheel shell
flywheel in in
shell experiment.
experiment.

4.2. Experimental Results and Analysis

According to the above experimental steps, the first straight-line segment, the third
curved segment, and the tenth circular segment were selected as typical path segments
for analysis. It can be seen in Figure 13 that the flywheel shell flash and burrs were greatly
removed by the robotic milling in comparison with the manual operation.
Metals 2022, 12, 1271 14 of 17

Figure 12. Robotic milling path of the flywheel shell in experiment.

4.2. Experimental Results and Analysis
4.2. Experimental Results and Analysis
According to the above experimental steps, the first straight-line segment, the third
According to the above experimental steps, the first straight-line segment, the third
curved segment, and the tenth circular segment were selected as typical path segments for
curved segment, and the tenth circular segment were selected as typical path segments
analysis.
for analysis. It canbebeseen
It can seen in Figure
in Figure 13the
13 that that the flywheel
flywheel shell flashshell flashwere
and burrs andgreatly
burrs were greatly
removed
removed bybythethe robotic
robotic milling
milling in comparison
in comparison with the with
manual the manual operation.
operation.

Path segment one Path segment three Path segment ten

Raw

Robotically
processed

Manually
processed

Figure
Figure13.13.
Comparison of theof
Comparison machining effects ofeffects
the machining flywheel
of shell flash/burrs
flywheel shellby robotic milling
flash/burrs and
by robotic milling and
manual operation.
manual operation.
The machining effects were further quantified by measuring the allowance height of
The machining
the flash/burrs effects were
on the straight-line (Str.), further quantified
curved (Cur.), by measuring
and circular the with
(Cir.) segments allowance height
of the flash/burrs on the straight-line (Str.), curved (Cur.),
a dial indicator. Figure 14 shows the measurement process as follows: (1) The dial and circular (Cir.) segments
with a dial
indicator indicator. Figure
was horizontally fixed on 14
theshows the measurement
lifting platform, and the needle process as touched
vertically follows: (1) The dial
the machined
indicator areahorizontally
was of the flywheel shell;on
fixed (2)the
the height
liftingofplatform,
the lifting and
platform
the was adjusted
needle vertically touched
tothe
measure the allowance height; (3) a total of 10 points were selected
machined area of the flywheel shell; (2) the height of the lifting platform in each path for was adjusted
measurement, and the average value was taken to characterize the machining quality of
to measure the allowance height; (3) a total of 10 points were selected in each path for
each path. As listed in Table 4, the average allowance height after robotic milling was 0.33
measurement,
Metals 2022, 12, x FOR PEER REVIEW and the average value was taken to characterize the machining 15 of 17 quality of
each path. As listed in Table 4, the average allowance height after robotic milling was
0.33 mm, with high consistency, which was about half of the value of 0.65 mm by the
conventional
mm, with highmanual operation.
consistency, which Additionally,
was about half inofterms of the
the value of machining
0.65 mm byefficiency,
the it took
conventional manual operation. Additionally, in terms of the machining efficiency,
about 9 min for the robot to process 90% of the flywheel shell flash/burrs at one time, it took
about 9the
while mintime
for the
is robot
up toto15process 90% ofmanual
min using the flywheel shell flash/burrs
operation. at one time,
Thus, robotic whileenhanced the
milling
the time is up to 15 min using manual operation. Thus, robotic milling enhanced the
machining efficiency by about 40%.
machining efficiency by about 40%.

Figure 14. Allowance height measurement after machining.

Figure 14. Allowance height measurement after machining.
Table 4. Comparison of the allowance height and time with two machining methods.

Robotic Milling Manual Operation

Path
Allowance Height (mm) Allowance Height (mm)
segment Time (s) Time (s)
Str. Cur. Cir. Str. Cur. Cir.
1 0.155 0.314 0.283 552 0.486 0.535 0.535 921
Metals 2022, 12, 1271 15 of 17

Table 4. Comparison of the allowance height and time with two machining methods.

Robotic Milling Manual Operation

Path Allowance Height (mm) Allowance Height (mm)
segment Time (s) Time (s)
Str. Cur. Cir. Str. Cur. Cir.
1 0.155 0.314 0.283 552 0.486 0.535 0.535 921
2 0.271 0.409 0.146 531 0.511 0.931 0.486 886
3 0.197 0.771 0.364 526 0.624 0.846 0.399 883
4 0.276 0.923 0.067 558 1.044 0.797 0.588 905
5 0.296 0.629 0.223 536 0.852 0.579 0.523 913
6 0.549 0.092 0.175 545 0.635 0.758 0.557 892
7 0.349 0.258 0.256 544 0.858 0.732 0.597 897
8 0.335 0.062 0.237 534 0.763 0.643 0.376 889
9 0.264 0.725 0.183 554 0.644 0.832 0.478 903
10 0.504 0.486 0.235 543 0.673 0.774 0.527 918

5. Conclusions
In this paper, a spatial path planning method for the robotic milling of an automotive
engine flywheel shell’s flash and burrs was proposed based on the optimal machining
posture and then verified with experiments. The following conclusions were achieved:
(1) The improved stereolithography slicing algorithm could intercept the complex com-
ponent contour in an arbitrary direction and eliminate the redundant points by sepa-
rating the inner and outer contours under certain constraints. The generated task path
accurately could match the model of the flywheel shell, showing excellent adaptability.
(2) The robotic posture evaluation function was established based on the stiffness perfor-
mance, supplemented by the robot’s collision, accessibility, and singularity indicators.
The optimal machining posture of the path point was determined according to the
frame rotation angle when the posture evaluation function value was the smallest.
(3) Both the simulation and experiment verified the feasibility of the proposed spatial
path planning method. Compared with manual operation, the robotic milling of the
flywheel shell’s flash and burrs exhibited better machining quality and consistency. In
particular, the machining time was reduced by 40%, compared with 15 min per piece
by hand.

Author Contributions: Conceptualization, H.W. and Y.W.; methodology, H.W.; software, H.W.;
investigation, H.W.; data curation, Y.W.; writing—original draft preparation, H.W.; writing—review
and editing, X.W.; supervision, D.Z.; funding acquisition, D.Z. All authors have read and agreed to
the published version of the manuscript.
Funding: This research was funded by the National Nature Science Foundation of China (No.51975443),
and Hubei Province Key R&D Program (No. 2020BAA025).
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: The raw/processed data required to reproduce these findings cannot
be shared at this time, as the data also form part of an ongoing study.
Conflicts of Interest: The authors declare no conflict of interest.

Nomenclature

Angi Rotation angle of the i-th intermediate point

Cfd 3 × 3 order translation compliance matrix
Cfδ 3 × 3 order coupling compliance matrix
Metals 2022, 12, 1271 16 of 17

Cmd 3 × 3 order coupling compliance matrix

Cmδ 3 × 3 order rotational compliance matrix
d Robot end coordinate system translation
E Overall posture evaluation function
Ereach Reachability posture evaluation function
Erigidity Stiffness posture evaluation function
Esingular Posture evaluation function based on singularity point
f Force vector on the end of the robot
KE Combined longitudinal and tangential directional stiffness
Kef Directional stiffness index
Kx Normal stiffness
Ky Tangential stiffness
m Torque at the end of the robot
STL Stereolithography
δ Robot end coordinate system rotation amount

References
1. Liu, Z.; Liu, Q.; Xu, W.; Wang, L.; Zhou, Z. Robot learning towards smart robotic manufacturing: A review. Robot. Comput.-Integr.
Manuf. 2022, 77, 102360. [CrossRef]
2. Tao, B.; Zhao, X.; Ding, H. Mobile-robotic machining for large complex components: A review study. Sci. China Technol. Sci. 2019,
62, 1388–1400. [CrossRef]
3. Jiang, Y.; Huang, Z.; Yang, B.; Yang, W. A review of robotic assembly strategies for the full operation procedure: Planning,
execution and evaluation. Robot. Comput.-Integr. Manuf. 2022, 78, 102366. [CrossRef]
4. Li, Y.; He, D.; Yuan, S.; Tang, K.; Zhu, J. Vector field-based curved layer slicing and path planning for multi-axis printing. Robot.
Comput.-Integr. Manuf. 2022, 77, 102362. [CrossRef]
5. Liang, F.; Kang, C.; Lu, Z.; Fang, F. Iso-scallop tool path planning for triangular mesh surfaces in multi-axis machining. Robot.
Comput.-Integr. Manuf. 2021, 72, 102206. [CrossRef]
6. Elber, G.; Cohen, E. Toolpath Generation for Freeform Surface Models. Comput. Aided Des. 1994, 26, 490–496. [CrossRef]
7. Suresh, K.; Yang, D. Constant Scallop-height Machining of Free-form Surfaces. J. Eng. Ind. 1994, 116, 253–259. [CrossRef]
8. Ding, S.; Mannan, M.; Poo, A.; Yang, D.; Han, Z. Adaptive iso-Planar Tool Path Generation for Machining of Free-Form Surfaces.
Comput. Aided Des. 2003, 35, 141–153. [CrossRef]
9. Yan, S.; Wang, S.; Tong, S. Uneven Offset Method of NC Tool Path Generation for Free-Form Pocket Machining. Comput. Ind.
2000, 43, 97–103. [CrossRef]
10. Ma, W.; But, W.; He, P. NURBS-based adaptive slicing for efficient rapid prototyping. Comput. Aided Des. 2004, 36,
1309–1325. [CrossRef]
11. Ma, X.; Lin, F.; Yao, B. Fast parallel algorithm for slicing STL based on pipeline. Chin. J. Mech. Eng. 2016, 29, 549–555. [CrossRef]
12. Zhang, Z.; Joshi, S. An improved slicing algorithm with efficient contour construction using STL files. Int. J. Adv. Manuf. Technol.
2015, 80, 1347–1362. [CrossRef]
13. Huang, S.; Zhang, L.; Han, M. An effective error-tolerance slicing algorithm for STL files. Int. J. Adv. Manuf. Technol. 2002, 20,
363–367. [CrossRef]
14. Xiao, W.; Huan, J. Redundancy and optimization of a 6R robot for five-axis milling applications: Singularity, joint limits and
collision. Prod. Eng. 2012, 6, 287–296. [CrossRef]
15. Zhu, W.; Qu, W.; Cao, L.; Yang, D.; Ke, Y. An off-line programming system for robotic drilling in aerospace manufacturing. Int. J.
Adv. Manuf. Technol. 2013, 68, 2535–2545. [CrossRef]
16. Xiong, G.; Ding, Y.; Zhu, L. Stiffness-based pose optimization of an industrial robot for five-axis milling. Robot. Comput.-Integr.
Manuf. 2019, 55, 19–28. [CrossRef]
17. Chen, Q.; Yang, J.; Yin, S.; Li, P. Optimization of Comprehensive Stiffness Performance Index for Industrial Robot in Milling
Process. In Proceedings of the 2019 IEEE 9th Annual International Conference on CYBER Technology in Automation, Control,
and Intelligent Systems (CYBER), Suzhou, China, 29 July–2 August, 2019.
18. Sun, J.; Zhang, W. Measured Stiffness Performance Evaluation Index of Industrial Robot in Side Milling. J. Phys. Conf. Ser. 2021,
1802, 022093. [CrossRef]
19. Chen, Q.; Zhang, C.; Hu, T.; Zhou, Y.; Ni, H.; Xue, X. Posture optimization in robotic machining based on comprehensive
deformation index considering spindle weight and cutting force. Robot. Comput.-Integr. Manuf. 2022, 74, 102290. [CrossRef]
20. Bcchet, E.; Cuillere, J.; Trochu, F. Generation of a Finite Element Mesh from Stereolithography (STL) Files. Comput. Aided Des.
2002, 34, 1–17. [CrossRef]
21. Piegl, L.; Tiller, W. The NURBS Book, 1st ed.; Springer: Berlin/Heidelberg, Germany, 1997; pp. 81–116.
22. Chen, C.; Peng, F.; Yan, R. Stiffness performance index based posture and feed orientation optimization in robotic milling process.
Robot. Comput.-Integr. Manuf. 2019, 55, 29–40. [CrossRef]
Metals 2022, 12, 1271 17 of 17

23. Huo, L.; Baron, L. The Self-Adaptation of Weights for Joint-Limits and Singularity Avoidances of Functionally Redundant
Robotic-Task. Robot. Comput.-Integr. Manuf. 2011, 27, 367–376. [CrossRef]
24. Fu, J.; Liu, H.; He, M.; Zhu, D. A Hand-Eye Calibration Algorithm of Binocular Stereo Vision Based on Multi-Pixel 3D Geometric
Centroid Relocalization. J. Adv. Manuf. Sci. Technol. 2022, 2, 2022005. [CrossRef]