materials

Article
An Efficient and Adaptable Path Planning Algorithm
for Automated Fiber Placement Based on Meshing
and Multi Guidelines
Hong Xiao 1, *, Wei Han 1 , Wenbin Tang 1,2 and Yugang Duan 1
1 State Key Lab for Manufacturing Systems Engineering, Xi’an Jiaotong University, Xi’an 710049, China;
jx18hanwei@[Link] (W.H.); tangwb@[Link] (W.T.); ygduan@[Link] (Y.D.)
2 School of Mechanical and Electrical Engineering, Xi’an Polytechnic University, Xi’an 710048, China
* Correspondence: xiaohongjxr@[Link]




Received: 11 August 2020; Accepted: 18 September 2020; Published: 22 September 2020


Abstract: Path planning algorithms for automated fiber placement are used to determine the directions
of the fiber paths and the start and end positions on the mold surfaces. The quality of the fiber paths
determines largely the efficiency and quality of the automated fiber placement process. The presented
work investigated an efficient path planning algorithm based on surface meshing. In addition,
an update method of the datum direction vector via a guide-line update strategy was proposed to
make the path planning algorithm applicable for complex surfaces. Finally, accuracy analysis was
performed on the proposed algorithm and it can be adopted as the reference for the triangulation
parameter selection for the path planning algorithm.

Keywords: fiber-reinforced polymers; automated fiber placement; path planning

1. Introduction
Fiber-reinforced polymers (FRPs), especially carbon fiber reinforced polymers (CFRPs), are widely
used in the aerospace and automobile industry as well as other fields because of their high
specific strength, high specific modulus, excellent corrosion- and fatigue-resistance, and outstanding
designability [1]. Automated fiber placement combines special manufacturing equipment with
computerized numerical control (CNC) systems, which enables good forming quality, strong
adaptability, and high efficiency. As an advanced composite manufacturing technology of large
and complex components, automated fiber placement has become one of the most advanced and
cutting-edge technologies for composite forming [2]. Path planning algorithms for automated fiber
placement are used to specify the direction of a fiber path and determine the start and end positions of
the fiber on the mold surface. The quality of the fiber path plays a crucial role for the efficiency and
quality of the fiber placement process.
At present, the typical fiber path planning algorithms include the geodesic method and the
meshing method based on parametric surfaces. The path planning algorithm of geodesic method based
on parametric surfaces includes analytical and numerical methods [3]. The fiber path, which is obtained
via the numerical solution of geodesics, is more aligned with practical engineering requirements.
Lewis et al. [4] proposed a path planning method based on natural path, which is also widely used
for path planning of tape laying. Shirinzadeh et al. [5–7] put forward the method of intersecting
lines and surfaces to build the initial path, which improves the adaptability of the path to the surface
curvature. However, there are certain disadvantages with the numerical solution of geodesics, e.g.,
the difficulty of finding the solution and the low efficiency associated with the calculation [8,9].
The meshing method, which uses the polygonal meshes to describe the complex parametric surfaces,

Materials 2020, 13, 4209; doi:10.3390/ma13184209 [Link]/journal/materials

Materials 2020, 13, 4209 2 of 18

Materials 2020, 13, x FOR PEER REVIEW 2 of 19

can reduce computational complexity. The deviation of the fiber path correlates with the accuracy of
thewith the accuracy
meshing [10,11].ofShinno
the meshing [10,11].
et al. [12] Shinno an
proposed et al. [12] proposed
iterative geodesic an algorithm
iterative geodesic algorithm
for a quadrilateral
for a quadrilateral mesh surface to obtain the fiber path. Li et al. [13] proposed
mesh surface to obtain the fiber path. Li et al. [13] proposed a path planning algorithm which a path planning
used
thealgorithm which used
mesh information the meshininformation
contained the STL [Link]
However, in the
thereSTLis file.
a lackHowever, there is analysis
of quantitative a lack of of
thequantitative
efficiency of analysis of the
the fiber efficiency
path planningof the fiber pathbased
algorithm planning algorithm
on meshed based onAlso,
surfaces. meshednosurfaces.
deviation
Also, no deviation analysis was performed for the generated fiber paths. Meanwhile,
analysis was performed for the generated fiber paths. Meanwhile, for complex components (surfaces), for complex
thecomponents (surfaces),
design of angle the direction
reference design ofdatum
angle reference direction
is too simple, which datum
leads is totoo
thesimple, which leadsoftothe
poor applicability
the poor applicability of the algorithm, fiber wrinkles, and eventually, affects the quality and
algorithm, fiber wrinkles, and eventually, affects the quality and mechanical properties of the fabricated
mechanical properties of the fabricated FRP components. In this paper, a path planning algorithm for
FRP components. In this paper, a path planning algorithm for automated fiber placement based on
automated fiber placement based on meshing and multi guide-lines was proposed. Both the
meshing and multi guide-lines was proposed. Both the efficiency and accuracy of the algorithm were
efficiency and accuracy of the algorithm were analyzed. The outline of the proposed algorithm is
analyzed. The outline of the proposed algorithm is shown in Figure 1.
shown in Figure 1.

Figure
Figure [Link]
Theoutline
outlineof
ofthe
theproposed
proposed algorithm
algorithm for
forautomated
automatedfiber
fiberplacement.
placement.

2. Efficiency Analysis
2. Efficiency and
Analysis andTopology
TopologyReconstruction
Reconstruction

[Link].
Efficiency Analysis
Efficiency Analysisofofthe
theProposed
ProposedAlgorithm
Algorithm
AA traditional
traditionalpathpathplanning
planningalgorithm
algorithm based
based on the geodesic
on the geodesic method
methodofofparametric
parametricsurface
surface
generates
generates fiber paths
fiber pathsbybysolving
solvingthe
thedirection
direction ofof each
each point on
on the
the surface.
[Link]
mainlyuses uses
thethe
twotwo
geometric
geometric numerical
numericalsolution
solutionoperations,
operations,i.e.,
i.e., the
the intersection betweenaasurface
intersection between surfaceandanda aplane
plane
andandthethe
parallel offset
parallel of of
offset a curve.
a [Link],
Essentially,the
theparallel
parallel offset
offset of a curve
curveisisrealized
realizedbybythe
theintersection
intersection of of
thethe
surface
surface at at
thethe samplingpoints
sampling pointsand
andthe
thecorresponding
corresponding offset direction
direction planes.
planes.
Commercial
Commercial CAD / CAM
CAD/CAM softwaresgenerally
softwares generallyuse useNURBS
NURBS (Non-Uniform
(Non-UniformRational
RationalB-Splines)
B-Splines)forfor
modelling, and the NURBS surface is represented by the following parameter
modelling, and the NURBS surface is represented by the following parameter equations: equations:

 = ( , )

 x = =sx (u,( v,) ) (1)

y ==s y (u,
( v, ) ) (1)




 z = s (u, v)
z
∑ ∑ , , , ( ) , ( )
s( , ) P
=Cu PCv W P N (u)N (v) (2)
i=0∑ j=∑0 i,j ,i,j ,i,p( ) , j,q ( )
s(u, v) = PC PC (2)
u v
where Ω is the domain of s, , is the control i=0 point,
j=0
Wi,j Ni,pis (uthe
)N weight,
j,q (v) ( ) is the basis function
, ,
of the pth-order B-spline in u direction, and , ( ) is the basis function of the qth-order B-spline in
where Ωs is the domain of s, Pi,j is the control point, Wi,j is the weight, Ni,p (u) is the basis function
v direction.
of the pth-order B-spline in u direction, and N j,q (v) is the basis function of the qth-order B-spline in
Let the plane be represented by an implicit surface equation:
v direction.
Let the plane be represented by an implicit , , ) =equation:
ℎ( surface 0 (3)
where Ω is the domain of h.
h(x, y, z) = 0 (3)
Combining Equation (1) with Equation (3) yields the intersection of plane h and surface s in
domain Ω = Ω ∩ Ω . The parameters u and v of the intersection line satisfy the equation for the
where Ωh is the domain of h.
intersection line [14]：
Combining Equation (1) with Equation (3) yields the intersection of plane h and surface s in
domain Ω = Ωh ∩ Ωs . The parameters ℎ ( , u ), and
( , v ),of the
( , intersection
) = ( , )= 0 satisfy the equation for(4)the
line
intersection line [14]:
Intersection line l is a plane
 curve in the plane of parameter  fields for u and v. If the surface s is
h sx (u, v), s y (u, v), sz (u, v) = l(u, v) = 0 (4)
the p × q-th order, the intersection line l is the p × q-th order. Because of the low efficiency of high-
Materials 2020, 13, 4209 3 of 18

Intersection line l is a plane curve in the plane of parameter fields for u and v. If the surface
s is the p × q-th order, the intersection line l is the p × q-th order. Because of the low efficiency of
high-precision floating-point calculations, the frequent and high-precision solution of the intersection
equation l consumes a lot of resources. This decreases both efficiency and stability of the numerical
algorithms for the intersection between surfaces and planes. In the path planning algorithm based on
the parametric surfaces, the intersection between surfaces and planes needs to be solved frequently.
This decreases the efficiency of the path calculation, and for complex surfaces, the equation order is too
high to be solved, which leads to the failure of path planning.
For triangular mesh surfaces, the surface is approximately a plane within the mesh cell [15].
The intersection between the surface and the plane can be converted into an intersection between
planes, and the straight line in the plane is the result of the intersection.
The intersection line equation of a triangular patch Ai (x, y, z) = 0 and the plane h is:

l0 (x, y, z) = 0 (5)


 x = x0 + mt

y = y0 + nt (6)




 z = z + pt
0

The intersection line l0 is the first order equation about t, and the solution complexity decreases
significantly. Since the triangulation algorithm will generate discrete grid planes of different sizes
with different model complexity, the algorithm will not increase the order and complexity of equation
solution during path generation due to the increase of model complexity. Therefore, the fiber path
planning algorithm based on triangular mesh surface can obtain stable numerical solution quickly
and efficiently.

2.2. Triangulation Algorithm for Parametric Surface

In the triangulation algorithm used in this paper, the edge and inner surfaces of the cell surface
are approximated by straight line segments. Hence, the surface is discretized into area strips, and it is
further divided into plane triangles. The triangulation results for the surface are generally given as
strips of triangles. A strip of triangles is a list of points, such that any three consecutive points define a
triangle. A parametric surface can be approximated by a series of triangle strips.
Due to the approximation process of replacing a curved surface with plane surfaces, the discrete
error occurs in the triangulation process for parametric surface. This error is mainly reflected in the
distance between the meshed plane surfaces and the original parametric surface. Two triangulation
parameters are used to constrain the approximation error: (a) Dl2s : The maximum distance from
the straight-line segment of the discrete triangular area to the original parametric surface. (b) Lseg :
The maximum length of the straight-line segment in the discrete triangular area.

2.3. Sub-surface Boundary Splicing and Surface Topology Reconstruction

NURBS curves and surfaces, which are widely used in CAD/CAM software, have exact
mathematical expressions, strong expression ability, good quality, and are easy to control. However,
for complex surfaces, NURBS surfaces need to be split, spliced and trimmed. A complex surface is
usually composed of multiple cellular surfaces. To reconstruct the topology of the whole complex
surface, in addition to the mesh reconstruction of the cellular parameter surfaces, boundary splicing
between cellular patches should be performed to obtain the global geometric topology.

2.3.1. Topology Reconstruction Algorithm for the Triangular Mesh of a Cellular Patch
As mentioned above, a cellular surface can be split into feature sampling points. To restore the
surface information, it is necessary to obtain the segmentation results and establish the data structure
for the triangular meshes by the vertex aggregation algorithm for triangulation (Algorithm 1).
Materials 2020, 13, 4209 4 of 18

Algorithm 1: Vertex aggregation algorithm for triangulation

Input: All vertex sets after the triangulation of a cellular parametric surface.
Output: Face ID List, Edge ID List and Point ID List for this cellular parametric surface.
1: Loop all strip, fans and triangles:
2: According to the right-hand rule, the vertices in the current discrete cell are stored to form the Point ID
List of the current cell parameter surface.
3: Loop all points in the Point ID List
4: Every three points in the point table form a triangle patch to summarize the Face ID list and store the
indexes of the three points of the current triangle patch.
5: Point ID List update, add triangle patch ID index.
6: Build the edge ID list, update the two-point indexes of the edge, update the Point ID list to add the edge
index, update the Face ID list to add the included edge index.
7: Calculate the normal vector of the current triangular patch according to the right-hand rule, and update it in
the Face ID list.

Through the above process, the local Point ID list, Edge ID list and Face ID list are established,
which reconstruct the topology information for the current NURBS cellular surface.

2.3.2. Algorithm for Subsurface-Boundary Splicing

Because the mesh discretization results for each cellular parameter surface are independent of
each other, and two adjacent cellular surfaces share the same edge, there are duplicate vertices for the
adjacent cellular surfaces. It is necessary to remove the duplicate vertices. The subsurface-boundary
splicing algorithm includes removing duplicate vertices and updating the global index of the vertices
in all the cellular surfaces.
The brute force algorithm, which traverses the whole Point ID List for the duplicate vertices
with the same coordinates with a given vertex and hence delete the duplicate ones, is the most
straightforward method. Suppose that the surface is made up of k quadrilateral cellular NURBS
surfaces with similar area and all the NURBS surfaces are smooth. Following the triangulation, it was
assumed that the triangulation results are all given as strips of triangles to obtain a uniform triangular
network. If the number of vertices on the boundary is a and b, the number of vertices on the cellular
surface is a × b, and the number of all vertices is N = k × a × b. The time complexity of the brute force
algorithm is O(N2 ).
However, all the duplicate vertices are located on the boundary line because they are formed by
two adjacent cellular surfaces sharing a common edge. For the vertices on the boundary line, they are
in a semi-closed state and not surrounded by all triangles. On the other hand, the vertices inside the
surfaces are in a fully-closed state, surrounded by several triangles (Figure 2). When the vertex is in
the fully-closed state, the number of adjacent patches is equal to the number of edges, while in the
semi-closed state, the number of patches is not equal to the number of edges. Therefore, all triangle
vertices can be divided into two types: boundary semi-closed vertices and internal fully-closed vertices.
Thanks to this feature, the boundary vertex set can be filtered out. Then, duplicate vertices can
be removed from the boundary vertex set, and Point ID list, Face ID list and Edge ID list can be
updated simultaneously.
 7: if [Link] (point2) < eps
8: delete point2 in Point ID List;
9: refresh Face ID List & Edge ID List;

The
Materials time
2020, 13, complexity
4209 of the algorithm is O (k × (a + b)2), which is far less than O (N2) of the 5brute
of 18
force algorithm, and the efficiency of the algorithm is greatly improved.

2. Internal fully-closed points and boundary semi-closed points

Figure 2. points in
in cell
cell parametric
parametric surfaces.
surfaces.

De-duplication
At the end of the algorithm
above based onboth
process, the vertex closureboundary
sub-surface checking was described
splicing in Algorithm
and surface 2.
topological
relation reconstruction are completed. This enables fast search of adjacent triangles and efficient path
Algorithm 2: De-duplication algorithm based on vertex closure checking
calculation. The topological information of the discrete mesh surfaces is contained in the Point ID list,
Input:
Face IDPoint ID list
list and Edgebefore
ID de-duplication.
list. The forms and requirements of the three data structure list are shown in
Output:
Table Point ID list after de-duplication.
1 [13].
1: Find semi-closed state vertices
Table 1. The forms and requirements of the three data structure lists.
2: Loop all points in all Point ID List:
3:Classif point->edgenum != point->facenum
Contents Features
4: save this point toint
duplicate
vIndex;Point ID List; Given the index number of a vertex, one can
5: Loop all point1 in duplicate Point ID
double List:
p [3]; quickly find the global index of the triangle
VertexLoop all point2 in duplicate Point ID List:
6: int EdgeIndexList [cur edge index]; patch and the global index for the edge,
7: if [Link] (point2) < eps
int FaceIndexList [cur face index]; where the current vertex belongs.
8: delete point2 in Point ID List;
Given the number of a side, one can quickly
9: intList
refresh Face ID eIndex;
& Edge ID List;
find the global index of the face, where the
Edge int VertexIndexList [2];
current edge belongs and the global indices of
The time complexity [2]; is O (k × (a + b)2 ), which is far less than O (N2 ) of the brute
of the algorithm
int FaceIndexList
the two vertices at the current edge.
force algorithm, and the efficiency of the algorithm is greatly improved.
int fIndex; Given the number of a triangle patch, one can
At the end of the above process, both sub-surface boundary splicing and surface topological
doubla n [3]; quickly find the global index and the global
relation
Face reconstruction are completed. This enables fast search of adjacent triangles and efficient path
int VertexIndexList [3]; index of the three edges of the three vertices
calculation. The topological information of the discrete mesh surfaces is contained in the Point ID list,
int EdgeIndexList [3]; on the triangle patch.
Face ID list and Edge ID list. The forms and requirements of the three data structure list are shown in
Table 1 [13].

Table 1. The forms and requirements of the three data structure lists.

Class Contents Features

int vIndex;
Given the index number of a vertex, one can quickly find the
double p [3];
Vertex global index of the triangle patch and the global index for
int EdgeIndexList [cur edge index];
the edge, where the current vertex belongs.
int FaceIndexList [cur face index];
int eIndex; Given the number of a side, one can quickly find the global
Edge int VertexIndexList [2]; index of the face, where the current edge belongs and the
int FaceIndexList [2]; global indices of the two vertices at the current edge.
int fIndex; Given the number of a triangle patch, one can quickly find
Face doubla n [3];int VertexIndexList [3]; the global index and the global index of the three edges of
int EdgeIndexList [3]; the three vertices on the triangle patch.

3. Main Path Planning Algorithm Based on Guidelines on Triangular Mesh Surfaces

The path planning process needs to generate the main motion paths of the head of the automated
fiber placement equipment and the corresponding fiber paths which are offset by the main motion
 The path planning process needs to generate the main motion paths of the head of the automated
fiber placement equipment and the corresponding fiber paths which are offset by the main motion
paths. The fiber path generation is relatively simple and this paper focuses on the main path planning
algorithm. The guidelines are extracted from the CAD model of the mold of the to-be-fabricated FRP
component to reflect the skeleton and unique appearance of the component. According to the guide-
Materials 2020, 13, 4209 6 of 18
lines, the direction vector of the main path can be calculated and adjusted adaptively, so that the main
path and the corresponding fiber paths can comply with the shape of the component. It is beneficial
to improve
paths. the mechanical
The fiber performances
path generation is relativelyofsimple
the FRP andcomponent
this paper while satisfying
focuses on the mainthe path
constraints
planningof
the minimum
algorithm. Thesteering radius
guidelines areofextracted
the fiber materials,
from the CADeven distribution
model of theofmold
the cutting
of thepoints, etc.
to-be-fabricated
FRP component to reflect the skeleton and unique appearance of the component. According to the
3.1. Main Motion
guide-lines, Path Generation
the direction vector Algorithm
of the main path can be calculated and adjusted adaptively, so that
the main path and the corresponding
According to the fiber placement process, fiber paths
threecan complyparameters
geometric with the shape
(initialofstarting
the component.
point P0,
It is beneficial to improve the mechanical performances of the FRP component while
guide line L, parametric surface Π) as well as the ply angle θ should be provided as input to the main satisfying the
constraints of the minimum steering radius of the fiber materials, even distribution
motion path generation algorithm, which outputs the main motion path li. Basically, the algorithm of the cutting
points,toetc.
needs calculate the direction vectors and subsequently generate the continuous points on the main
motion path, as described in Algorithms 3 and 4, respectively.
3.1. Main Motion Path Generation Algorithm
Algorithm 3: The direction vector calculation algorithm for the main motion path.
According
Step 1: Findtothe theprojection
fiber placement
point P0process,
’ for the three
initialgeometric parameters
starting point P0 on the (initial starting
triangle patchpoint P0 ,
A of the
guide line L,
triangular parametric
mesh surfacessurface Π) as well of
∑ (triangulation asthe ply angle θsurface
the parametric should Π)be provided
and obtain asthe
input to thevector
normal main
motion
n path patch
of triangle generation
A. Then,algorithm,
redefinewhich outputs
P0’ as the the point
starting main for
motion path li .main
the current Basically,
motionthepath.
algorithm
needsStepto calculate
2: Calculatethe direction
the tangent vectors
vector and subsequently
t’ for the projectiongenerate
point the
P0’’continuous
of P0’ on thepoints
guideon theL,main
line and
motion path,
generate as described
the parallel vectorint Algorithms
of t’ at point3Pand0’. 4, respectively.

Algorithm
Step 3: Thethe
3: Calculate direction vector
orthogonal calculation
vector k of thealgorithm for the
normal vector main
n and themotion
parallelpath.
tangent vector
Stepk1:
t, where × t. the projection point P0 ’ for the initial starting point P0 on the triangle patch A of the
= nFind
Π) and
P
triangular
Step 4:mesh surfaces
Calculate (triangulation
the orthogonal of the
vector m ofparametric surface
normal vector n and obtain
vector the normal
k, where m = k ×vector
n. Then
of triangle
vector m ispatch A. Then,direction
the datum P0 ’ as the
redefine vector for starting pointpath
the current for the current
point, main
based onmotion path.
the guide-line and
Step by
corrected 2: the
Calculate
curvaturethe feature
tangentofvector t’ for the
the surface, whereprojection pointPis
the pathpoint of P0 ’ on
0 ”located (as the guide
shown line L,
in Figure
and generate the parallel vector t of t’ at point P0 ’.
3).
Step 3:
5: Calculate
Using normal the orthogonal
vector n as vector k of the normal
the rotation then datum
vector
axis, rotate and thedirection
parallel tangent m for θt,
vector vector
degreekto=obtain
where n × t. the laying direction vector d for the current path point P0’.
Step 4: Calculate the orthogonal vector cosm of−sin normal0vector n and vector k, where m = k × n.
The vector m is the datum direction vector = sin for thecos × point, based on the guide-line and
current0 path (7)
corrected by the curvature feature of the surface, 0 where 0 the path
1 point is located (as shown in Figure 3).

Figure 3. Illustration of the direction vector calculation algorithm.

Figure 3. Illustration of the direction vector calculation algorithm.

Step 5: Using normal vector n as the rotation axis, rotate the datum direction vector m for θ degree
to obtain the laying direction vector d for the current path point P0 ’.

 cos θ − sin θ 0 
 

d =  sin θ cos θ 0  × m

 
(7)
 
0 0 1

In the current triangle patch A, use point P00 and vector d to construct a straight line and
hence obtain the intersection point P1 of the straight line and the three sides of triangle patch A.
Materials 2020, 13, 4209 7 of 18

The intersection point P1 is the next path point, and all the path points can be obtained after continuous
iteration for the main motion path generation.
Algorithm 4: The continuous path point generation algorithm.
Step 1: Take the starting point Pi and direction vector di in the current triangle patch Ai , construct
a ray Pa (λ), which is defined as the solution line for the path point. Its parameter equation is:

Pa (λ) = Pi + λ·d (8)

Step 2: The vertices Va and Vb of the triangle patch form a straight line Pb (λ), and the parameter
equation is:
Pb (λ) = (1 − λ)·Va + λ·Vb (9)

Step 3: Using Pa (λ) = Pb (λ), we can find the unique solution λ:

Va − Pi
λ= → (10)
Va −Vb
d+ →
kVa −Vb k

If 0 ≤ λ ≤ 1 is true, this means that the next path point is on the current sideline. If it is not on the
current sideline, then take two points to form a straight line for the calculation. Next, take Vc Vb and
Va Vc to form a straight line for the calculation. Repeat steps 1, 2, and 3 to get the next path point.
The solution line for the path point may overlap with the three side lines of the triangle in step 3,
and the next path point Pi+1 cannot be calculated using the above steps. At this time, the endpoint of
the edge Vm Vk (m, k ∈ {a, b, c}) of the triangle, where the point Pi is located, is used as the next path point
Pi+1 . The endpoint selection rules are described in Equation (11), where d is the laying direction vector.
(
Vk , d ∗ Vm Vk = 1
Pi + 1 = (11)
Vm , otherwise

Step 4: According to the topological relationship for the triangular mesh surface, obtain the
adjacent triangles on the other side of the edge, where Pi+1 is located, and update the patch index
as Ai+1 .
Step 5: Update the normal vector ni+1 . When the next path point Pi+1 lies on the triangle edge or
vertex, use the area weighing method to reduce the deviation.
Pm
Ak ·nk
n = Pkm=1 (12)
k=1 Ak ·nk

In Equation (12), m is the total number of adjacent triangles to which Pi+1 belongs, and Ak is the
area of the triangles. The area can be determined using Ak = kAB×ACk 2 , where points A, B, C are the
three vertices of a triangle Ak .
Step 6: Update the parallel tangent vector ti+1 , and solve the direction vector di+1 of the next path
point Pi+1 using Algorithm 3. The projection vector di+10 of di+1 on patch Ak+1 is the laying-direction
vector of point Pi+1 .
di+1 ·n
d0i+1 = di+1 − ·n (13)
knk2
Step 7: Repeat Step 1 to Step 6 to update and calculate the next path point until reach the triangular
mesh surface boundary and no adjacent patch can be found and the partial path on the one single
direction is constructed.
Step 8: Go back to the path initial point P0 , take the inverse of the laying-direction vector d0 and
repeat Step 1 to Step 7 until the whole main motion path is completely constructed.
Materials 2020, 13, 4209 8 of 18

The
Materials continuous
2020, pathREVIEW
13, x FOR PEER point generation algorithm is shown in Figure 4. After completing the above
8 of 19
steps, the main motion path can be derived using the cubic spline interpolation algorithm.
Materials 2020, 13, x FOR PEER REVIEW 8 of 19

Figure 4. Continuous path point generation algorithm. (a) Calculation of next point in single triangle
Figure 4. Continuous path point generation algorithm. (a) Calculation of next point in single triangle
Figure(b)
patch. 4. Continuous
Calculation path
of allpoint generation
points algorithm.
of one path (a) Calculation
on the triangular of next point in single triangle
mesh surfaces.
patch. (b) Calculation of all points of one path on the triangular mesh surfaces.
patch. (b) Calculation of all points of one path on the triangular mesh surfaces.
3.2.
3.2. Guide-line
Guide-line Update
Update Algorithm
Algorithm forfor Complex
Complex Surfaces
Surfaces
3.2. Guide-line Update Algorithm for Complex Surfaces
If
If the
the curvature
curvature ofof the
the surface
surface of of certain
certain components
components changes
changes substantially,
substantially, the
the path,
path, which
which was
was
If the
calculated
calculated curvature
and plannedof the surface
using a singleof certain
guide components
line, cannot changes
meet the substantially,
requirements the
for path,
the
planned using a single guide line, cannot meet the requirements for the laying ability in which
laying was
ability
calculated
in each
each surfaceand
surface planned
area.
area. This
This using
leads
leads a single
to wrinklingguide
to wrinkling of line,
of the
the cannot
fiber
fiber meet
during
during thethe
the requirements
laying
laying for
process,
process, the laying
which
which ability
degrades
degrades the
in
theeach surface
mechanical
mechanical area. This
properties
properties leads
of the thetofabricated
offinal wrinkling
final of the
fabricated fiber To
during
components.
components. Totheaddress
address laying
this process,
thiswhich
this problem,
problem, degrades
paperthis paper
proposed
the mechanical
proposed
a guide-line properties
a guide-line
algorithm of update
algorithm
for the theforfinaloffabricated
the update
the of components.
tangent the tangent
vector To
t andvectoraddress
datum t and this vector
datum
direction problem,
directionthis
m for thepaper
vector m
path
proposed
for a guide-line
the according
points path points algorithm
according
to the surface to for the
the and
shape update
surface of the tangent
shape and the
the distribution vector t and
distribution
of multi datum
of multi
guide-lines, direction vector
guide-lines, as
as demonstrated m
in
for the 5path
demonstrated
Figure points according
in Figure
and described 5 and
in to the5. surface
described
Algorithm shape and
in Algorithm 5. the distribution of multi guide-lines, as
demonstrated in Figure 5 and described in Algorithm 5.

Guide-line
Figure [Link]-line
Figure update
update algorithm
algorithm for complex
for complex surfaces:
surfaces: (a) Distribution
(a) Distribution diagramdiagram forguide-
for multi multi
guide-lines.
Figure(b)
lines. (b)
5. PolygonalPolygonal
Guide-linesection section
updateformed
algorithmformed by the
for complex
by the projection
surfaces:
projection point of(a) point of the
theDistribution point
point on thediagram on the guide-line.
for multi
guide-line. guide-
(c) Diagram
(c)tangent
of Diagram
lines. of tangent
(b) Polygonal
vector vector
section calculation.
formed
calculation. by the projection point of the point on the guide-line. (c) Diagram
of tangent vector calculation.
Algorithm 5: Guideline update algorithm for complex surfaces.
Materials 2020, 13, x FOR PEER REVIEW 9 of 19
Materials 2020, 13, 4209 9 of 18
Algorithm 5: Guideline update algorithm for complex surfaces.
Step 1: Project the current path point Pi on each guide-line Li to obtain projection points A, B, C
Step 1: Project the current path point Pi on each guide-line Li to obtain projection points A, B, C
and D. Connect the projection points to form a polygonal section plane Ω, the four sides of the
and D. Connect the projection points to form a polygonal section plane Ω, the four sides of the polygon
polygon on the section plane Ω are a, b, c and d.
on the section plane Ω are a, b, c and d.
Step 2: Project the current path point Pi to the edge of the polygon, and find the projection points
Step 2: Project the current path point Pi to the edge of the polygon, and find the projection points
, , and .
P1i , P2i , P3i and P4i .
Step 3: If the projection point is on the extension line of the edge, the edge is excluded. As shown
Step 3: If the projection point is on the extension line of the edge, the edge is excluded. As shown
in Figure 5b, the projection point 4 is not on the edge d, so the edge d is not considered in the
in Figure 5b, the projection point Pi is not on the edge d, so the edge d is not considered in the
algorithm below.
algorithm below.
Step 4: Calculate the distance from Pi to the projection points on the not excluded edge lines,
Step 4: Calculate the distance from Pi to the projection points on the not excluded edge lines,
where = ‖ ‖.
where dist = kPi Pxi k.
Step 5: Determine the minimum distance dist and its edge x. The two guide-lines La and Lb, at the
Step 5: Determine the minimum distance dist and its edge x. The two guide-lines La and Lb , at the
end of x, are the guide-lines for the current path point Pi.
end of x, are the guide-lines for the current path point Pi .
Step 6: As shown in Figure 5c and Equation (14), calculate the tangent vectors on La and Lb,
Step 6: As shown in Figure 5c and Equation (14), calculate the tangent vectors on La and Lb ,
respectively. The smooth transition for the path direction vector from La to Lb is realized using the
respectively. The smooth transition for the path direction vector from La to Lb is realized using the
distance-weighing method:
distance-weighing method:
∙ 1 ++t2 ·dist
t ·dist ∙ 2
t= 1 (14)
(14)
dist1 ++dist2
dist11 and
where dist
where dist22are
and dist arethe
thedistances
distancesof ofPPi ito
toLLa aand
andLLb.b .
Update the tangent vector t
Update the tangent vector t and calculate thedatum
and calculate the datumdirection vectormmofof
directionvector PiPin Algorithms
i in Algorithms3 3and and4.
The path, which is generated by Algorithm 5, is more adaptive to the changing
4. The path, which is generated by Algorithm 5, is more adaptive to the changing curvature of the curvature of the surface
and theand
surface scalability is improved.
the scalability FigureFigure
is improved. 6 shows the paths
6 shows generated
the paths basedbased
generated on single guideline
on single and
guideline
multi-guidelines for a panel and a curved surface models. It can be seen that the
and multi-guidelines for a panel and a curved surface models. It can be seen that the fiber directionfiber direction transits
smoothly
transits betweenbetween
smoothly the guidethe lines andlines
guide this makes
and this themakes
placedthefibers adapt
placed to the
fibers shape
adapt of the
to the moulds
shape of
of the
the finalof
moulds parts.
the final parts.

Figure 6.
Figure The paths
6. The paths generated
generated based
based on
on single
single guide
guide line
line and
andmulti-guide
multi-guidelines
linesfor
foraapanel
panel(a,b)
(a,b)and
anda
curved surface (c,d) models.
a curved surface (c,d) models.
Materials 2020, 13, 4209 10 of 18
Materials 2020, 13, x FOR PEER REVIEW 10 of 19

4. Accuracy Analysis of the Generated Path Based on Surface Meshing

In the triangulation
triangulation process
processin inSection
Section2,2,DDl2sl2sand
and Lseg
Lseg areare
thethe main
main parameters
parameters thatthat affect
affect bothboth
the
the mesh density and the approximation accuracy. When the two parameters
mesh density and the approximation accuracy. When the two parameters are set to larger values, theare set to larger values,
the
meshmesh density
density is small,
is small, thethe surface
surface approximation
approximation accuracyisislow,
accuracy low,the
the subsequent
subsequent path
path planning
process data volume is small, and the algorithm efficiency is high. high. If the two parameter values are
continuously reduced, on the other hand, hand, thethe efficiency
efficiency is is lower.
lower.
It is generally believed that the path generated using the geodesic method on the parametric
surface is used as the standard path when the planning accuracy and error of the algorithm need to
be verified. The
Thepoints
pointson onthe
thediscrete
discretesurface
surface ofofthethe mesh
mesh areare used
used to replace
to replace thethe
pathpath points
points on
on the
the original
original parametric
parametric surface.
surface. Furthermore,
Furthermore, the normal
the normal vectorvector
of theofpath
the points
path points
on theon the triangular
triangular mesh
mesh surface can be used to replace the normal vector on the original parametric
surface can be used to replace the normal vector on the original parametric surface. As a result, a surface. As a result,
a cumulativeerror
cumulative errorcan
canoccur
occurduring
duringthetheiteration
iterationprocess
processof ofthe
the proposed
proposed path path planning
planning algorithm,
algorithm,
which can cause the angle deviation from the design datum and the distance deviation from the path
generated on the original parametric surface.
A complex
complexsurface
surfacewith
withpositive andand
positive negative
negative curvature was adopted
curvature to conduct
was adopted accuracyaccuracy
to conduct analysis
of the proposed
analysis path planning
of the proposed path algorithm, as shown inasFigure
planning algorithm, shown 7. in
The surface
Figure consists
7. The of 29consists
surface independent
of 29
cellular
independent cellular patches. Each cellular parameter surface represents a NURBS surface, andeach
patches. Each cellular parameter surface represents a NURBS surface, and the order of the
cellular
order ofparameter surface
each cellular in both
parameter U andin
surface V both
directions
U andisV6directions
degrees. is 6 degrees.

Figure [Link]
Figure Illustrationof of
thethe
surface for accuracy
surface analysis.
for accuracy (a) Original
analysis. parametric
(a) Original surface. (b)
parametric Triangular
surface. (b)
mesh surface
Triangular meshwith
surface = 0.2Dl2s
Dl2s with and Lseg
= 0.2 and= L100. (c) Surface curvature variation diagram for the
seg = 100. (c) Surface curvature variation diagram for the
path-deviation-error test area.
path-deviation-error test area.
4.1. Distance Deviation Analysis
4.1. Distance Deviation Analysis
The distance deviation of the path generated on the triangular mesh surface can be decomposed
The normal
into the distancedistance
deviation of the path
deviation andgenerated on the
the geodesic triangular
distance mesh surface
deviation. can
In this be decomposed
paper, a uniform
orthogonal test was carried out for the appropriate ranges of Dl2s and Lseg . The path was auniformly
into the normal distance deviation and the geodesic distance deviation. In this paper, uniform
orthogonal test was carried out for the appropriate ranges of Dl2s and Lseg. The path was uniformly
sampled (with a distance of 5 mm) to evaluate the distance deviation, as shown in Figure 8.
sampled (with
During thea distance of 5 mm)
triangulation to evaluate
process, both Dthe distance deviation, as shown in Figure 8.
l2s and Lseg did not reach 0, which means that the
Duringvalue
parameter the triangulation process, both
(without approximation Dl2scould
error) and Lnot
seg did not reach 0, which means that the
be determined. Therefore, in the actual
parameter value (without approximation error) could
application, the Dl2s range was 0.2–1.0 mm, and the Lseg range not be determined.
was 20–200Therefore,
mm. Theinexperiment
the actual
application, the Dl2s range was  0.2–1.0 mm, and the Lseg range was 20–200 mm. The experiment was
was carried out by the L25 56 orthogonal design. The path-generation time, the distance deviation,
carried out by the (5 ) orthogonal design. The path-generation time, the distance deviation, the
the normal distance deviation and the geodesic distance deviation were recorded.
normal distance deviation and the geodesic distance deviation were recorded.
The Euclidean distance deviation d from the sample point on the generated path to the standard
reference path was calculated and decomposed into the normal distance deviation dN and the geodesic
distance deviation dT .
Materials 2020, 13, 4209 11 of 18
Materials 2020, 13, x FOR PEER REVIEW 11 of 19

Figure
Figure 8. Standardpath,
8. Standard path,verification
verificationpath
pathgenerated
generatedbyby the
the new
new algorithm,
algorithm, andand sample
sample points
points for
for D l2s
D l2s = 0.4 and L seg
= 0.4 and Lseg = 100. = 100.

The project vector from the sample point to standard reference path is V, and the normal vector
The Euclidean distance deviation d from the sample point on the generated path to the standard
for sample point Pi on the original parametric surface is n. dN and dT can be calculated as follows:
reference path was calculated and decomposed into the normal distance deviation dN and the
geodesic distance deviation dT.   

 d = DistBetween Pi , P0i
The project vector from the sample
 dpoint

 to standard reference path is V, and the normal vector
N = d ×q cos < V, n >

for sample point Pi on the original parametric


 surface is n. dN and dT can be calculated as follows:(15)
dT = d2 − d2N


= ( , )


  = × < , >
According to the L25 56 orthogonal design, 25 sets of experiments were carried out. Four(15) of
them, which were Dl2s = 0.2 and Lseg = 65, Dl2s = 0.4 − and Lseg = 110, Dl2s = 0.6 and Lseg = 155 and
Dl2s = 0.8 and Lseg = 200, were selected to analyze the distribution of distance deviation along the
According to the (5 ) orthogonal design, 25 sets of experiments were carried out. Four of
generated paths.
them, which were Dl2s = 0.2 and Lseg = 65, Dl2s =0.4 and Lseg = 110, Dl2s = 0.6 and Lseg = 155 and Dl2s = 0.8
The initial path point was located near the 350th sample point. According to Figure 9, the closer the
and Lseg = 200, were selected to analyze the distribution of distance deviation along the generated
sample point was to the initial point, the smaller was the distance deviation. During path generation,
paths.
the normal vector for the path points on the triangular mesh surface was used to (approximately)
The initial path point was located near the 350th sample point. According to Figure 9, the closer
replace the normal vector of the path points on the original parameter surface. Therefore, the direction
the sample point was to the initial point, the smaller was the distance deviation. During path
datum of the path point was deviated, and it caused geodesic distance deviation between the generated
generation, the normal vector for the path points on the triangular mesh surface was used to
path and the standard reference path. It also shows that the geodesic distance deviation dT was the
(approximately) replace the normal vector of the path points on the original parameter surface.
main deviation and it increased as the parameters Dl2s and Lseg increased.
Therefore, the direction datum of the path point was deviated, and it caused geodesic distance
To analyze the relationship between the parameters (i.e., Dl2s and Lseg ) and the distance deviation of
deviation between the generated path and the standard reference path. It also shows that the geodesic
the generated path and the algorithm efficiency, the mean distance deviation, the normal mean distance
distance deviation dT was the main deviation and it increased as the parameters Dl2s and Lseg increased.
deviation, the geodesic mean distance deviation, the maximum distance deviation, the maximum
To analyze the relationship between the parameters (i.e., Dl2s and Lseg) and the distance deviation
normal distance deviation, the maximum geodesic distance deviation, and the path generation time
of the generated path and the algorithm efficiency, the mean distance deviation, the normal mean
were calculated. This was done when the generation time started from the beginning of the surface
distance deviation, the geodesic mean distance deviation, the maximum distance deviation, the
triangulation to the end of the path generation, as shown in Table 2.
maximum normal distance deviation, the maximum geodesic distance deviation, and the path
generation time were calculated. This was done when the generation time started from the beginning
of the surface triangulation to the end of the path generation, as shown in Table 2.
Materials 2020, 13, 4209 12 of 18

Materials 2020, 13, x FOR PEER REVIEW 12 of 19

Figure 9. The distance deviation distribution (N Distance is dN and T Distance is dT.) for different
Figure 9. The distance deviation distribution (N Distance is dN and T Distance is dT .) for different
parameters (a) Dl2s = 0.2,Lseg = 65, (b) Dl2s = 0.4,Lseg = 110, (c) Dl2s = 0.6,Lseg = 155 and (d) Dl2s = 0.8, Lseg =
parameters (a) Dl2s = 0.2, Lseg = 65, (b) Dl2s = 0.4, Lseg = 110, (c) Dl2s = 0.6, Lseg = 155 and (d) Dl2s = 0.8,
200.
Lseg = 200.
[Link]
Table Distance deviation
deviationand
andgeneration
generationtime.
time.
Normal Geodesic Maximum
Maximum
Maximum Maximum
Mean Normal Mean
Mean Geodesic Maximum
Mean MeanMean Maximum Normal
NormalGeodesicGeodesic
GenerationGeneration
Dl2s Dl2s LsegLseg Distance Distance
Distance Distance Distance
Distance
DistanceDistanceDistanceTime Time
/mm Distance Distance Deviation Distance
/mm /mm /mm Deviation
Deviation Deviation Deviation
Deviation Deviation Deviation /s
/mm /mm
Deviation /mm
Deviation /mm Deviation
/mm /mm /mm Deviation /mm /s
0.2 20 0.41072
/mm
0.00824
/mm
0.40834 0.81222
/mm 0.05218 /mm 0.81219 5.506
0.2 0.2 65 20 0.41072
0.45268 0.00824
0.02694 0.40834
0.4329 0.81222
1.29759 0.05218
0.16468 0.81219 1.29745 5.506 0.549
0.2 0.2 11065 0.57411
0.45268 0.02534
0.02694 0.55406
0.4329 1.60354
1.29759 0.16447 1.29745 1.60342 0.549
0.16468 0.508
0.2 155 0.57298 0.02557 0.55294 1.60354 0.16447 1.60342 0.508
0.2 0.2 200110 0.57411
0.57309 0.02534
0.02547 0.55406
0.55304 1.60354
1.60354 0.16447
0.16447 1.60342 1.60342 0.508 0.505
0.4 0.2 20155 0.57298
0.41072 0.02557
0.00824 0.55294
0.40834 1.60354
0.81222 0.16447
0.05218 1.60342 0.81219 0.508 5.514
0.4 0.2 65200 0.54446
0.57309 0.09072
0.02547 0.47922
0.55304 2.12237
1.60354 0.36747 1.60342 2.12172 0.505
0.16447 0.296
0.4 110 0.78375 0.09412 0.71355 3.00629 0.46252 3.00561 0.214
0.4 0.4 15520 0.41072
0.78403 0.00824
0.09052 0.40834
0.71476 0.81222
2.94274 0.05218
0.46246 0.81219 2.94204 5.514 0.209
0.4 0.4 20065 0.54446
0.78403 0.09072
0.09052 0.47922
0.71476 2.12237
2.94274 0.36747
0.46246 2.12172 2.94204 0.296 0.213
0.6 0.4 20110 0.41072
0.78375 0.00824
0.09412 0.40834
0.71355 0.81222
3.00629 0.05218 3.00561 0.81219 0.214
0.46252 5.495
0.6 65 0.62143 0.1036 0.55255 2.14027 0.55588 2.13962 0.231
0.6 0.4 110155 0.78403
0.65292 0.09052
0.20369 0.71476
0.4915 2.94274
1.62566 0.46246
0.77333 2.94204 1.62391 0.209 0.142
0.6 0.4 155200 0.78403
0.86012 0.09052
0.18893 0.71476
0.69804 2.94274
3.09987 0.46246
0.77349 2.94204 3.0993 0.213 0.13
0.6 0.6 20020 0.8602
0.41072 0.18895
0.00824 0.6985
0.40834 3.08182
0.81222 0.77349 0.81219 3.08124 5.495
0.05218 0.128
0.8 20 0.41072 0.00824 0.40834 0.81222 0.05218 0.81219 5.499
0.8 0.6 65 65 0.62143
0.62315 0.1036
0.10188 0.55255
0.55424 2.14027
2.19272 0.55588
0.55588 2.13962 2.19207 0.231 0.213
0.8 0.6 110110 0.65292
0.78546 0.20369
0.20771 0.4915
0.63916 1.62566
1.94612 0.77333
0.84012 1.62391 1.94459 0.142 0.115
0.8 155 1.18547 0.19293 1.03511 3.91933 0.84004
0.6 155 0.86012 0.18893 0.69804 3.09987 0.77349 3.0993 3.91922 0.13 0.098
0.8 200 1.15638 0.19065 1.00657 3.8031 0.84004 3.80299 0.103
1 0.6 20200 0.8602
0.41072 0.18895
0.00824 0.6985
0.40834 3.08182
0.81222 0.77349
0.05218 3.08124 0.81219 0.128 5.491
1 0.8 65 20 0.41072
0.62147 0.00824
0.10314 0.40834
0.55262 0.81222
2.14027 0.05218
0.55588 0.81219 2.13962 5.499 0.203
1 110 0.83477 0.2216 0.68092 2.01082 1.10312
0.8 65 0.62315 0.10188 0.55424 2.19272 0.55588 2.19207 2.00955 0.213 0.105
1 155 1.33796 0.21533 1.16792 3.97448 1.10346 3.97407 0.082
1 0.8 200110 0.78546
1.29659 0.20771
0.20791 0.63916
1.12927 1.94612
3.81894 0.84012
1.10346 1.94459 3.8185 0.115 0.083

The distribution of mean distance deviation are shown in Figure 10, while the maximum
distance deviation are shown in Figure 11, and the generation time are shown in Figure 12.
 0.8
1 200
65 1.15638
0.62147 0.19065
0.10314 1.00657
0.55262 3.8031
2.14027 0.84004
0.55588 3.80299
2.13962 0.103
0.203
11 20
110 0.41072
0.83477 0.00824
0.2216 0.40834
0.68092 0.81222
2.01082 0.05218
1.10312 0.81219
2.00955 5.491
0.105
11 65
155 0.62147
1.33796 0.10314
0.21533 0.55262
1.16792 2.14027
3.97448 0.55588
1.10346 2.13962
3.97407 0.203
0.082
11 110
200 0.83477
1.29659 0.2216
0.20791 0.68092
1.12927 2.01082
3.81894 1.10312
1.10346 2.00955
3.8185 0.105
0.083
1 155 1.33796 0.21533 1.16792 3.97448 1.10346 3.97407 0.082
1 The 200 1.29659of mean
distribution 0.20791
distance1.12927
deviation 3.81894
are shown 1.10346
in Figure 10,3.8185
while the 0.083
maximum
Materials 2020, 13, 4209 13 of 18
distance deviation are shown in Figure 11, and the generation time are shown in Figure 12.
The distribution of mean distance deviation are shown in Figure 10, while the maximum
distance deviation
Supplementary VideoareS1shown indemonstrates
further Figure 11, andthe
thehigh
generation time
efficiency of are
the shown in Figure
proposed 12. which can
algorithm,
complete the path planning for one layer of a complex surface in just only tens of seconds.

Figure 10. Mean distance deviation and its decomposition.

Figure [Link]
Figure10. Mean distance deviationand
distance deviation andits
itsdecomposition.
decomposition.

Figure [Link]
Figure11. Maximumdistance
distance deviation and its
deviation and its decomposition.
decomposition.

According to Figures 10 and

Figure 11, for different
11. Maximum distance parameters,
deviation and the distance deviation mainly consisted
its decomposition.
of geodesic distance deviation, and the effect of normal distance deviation on the overall distance
deviation is relatively small. Meanwhile, it can be observed that Dl2s and Lseg restricted each other in
the triangulation process. When the two parameters cannot be satisfied simultaneously, the algorithm
will adopt the parameter that makes the mesh more precise. For example, when Lseg is 20 mm,
the triangulation algorithm generated the same triangular mesh surfaces for a Dl2s of 0.4 mm, 0.6 mm,
0.8 mm and 1.0 mm as for the Dl2s of 0.2 mm.
Table 2 shows that, when Dl2s exceeds 0.8mm and Lseg exceeds 155 mm, the mean distance
deviation surpasses 1mm, and the maximum distance deviation exceeds 2 mm. When Dl2s is between
 Materials 2020, 13, 4209 14 of 18

0.2 and2020,
Materials 0.6mm, and PEER
13, x FOR Lseg is 65 to 155 mm, the average distance deviation remains within 1mm, while
REVIEW 14 of the
19
maximum distance deviation is within 2 mm. The generation time of the algorithm is within 0.5 s.

Figure 12. The generation times.

Figure 12. The generation times.
4.2. Angle Deviation Analysis
According
To ensuretotheFigures 10 and 11,
orthogonal for differentof
arrangement parameters, the distance
fibers in different deviation
layers, mainlydistribution
the angle consisted
ofbetween
geodesic distance
different deviation,
layers anddefined.
is strictly the effect of ensures
This normal distance deviation
quasi-isotropic on the overall
mechanical behavior distance
for the
deviation is relatively small. Meanwhile, it can be observed that D and
components [3]. To analyze the angle deviation between the generated paths and the design datum,
l2s Lseg restricted each other in
the
antriangulation
angle deviationprocess. When
analysis wasthe two parameters
performed cannot be satisfied simultaneously, the algorithm
in this section.
will adopt the parameter that makes the
The generated path was sampled uniformly with mesh more precise. For of
a step example,
5 mm. when Lseg issample
For each 20 mm,point,
the
triangulation algorithm generated the same triangular mesh surfaces for a D
the forward direction of the current path and the design direction datum of the path were calculated
l2s of 0.4 mm, 0.6 mm, 0.8
mm andboth
using 1.0 mm as for the D
the guide-line l2s of
and 0.2 mm.
surface normal vector at the point. Subsequently, the angle deviation
Table 2 shows that, when D
was obtained. The forward direction vector di0.8mm
l2s exceeds at eachand sample Lseg point
exceeds 155
is the mm, the
tangent mean
vector distance
of the sample
deviation
point on surpasses
the path. The1mm, and the
datum maximum
direction vectordistance
0
di can be deviation exceeds
calculated using:2 mm. When Dl2s is between
0.2 and 0.6mm, and Lseg is 65 to 155 mm, the average distance deviation remains within 1mm, while
( 0
the maximum distance deviation is within 2 mm. ki =Then0i ×generation
t0i time of the algorithm is within 0.5
0 0 0 (16)
s. di = ni × ki

where
4.2. n0i Deviation
Angle is the normal vector of the sample point on the surface, and t0i is the tangent of the projection
Analysis
point for the sample point on the guide-line. The angle deviation δi is calculated using:
To ensure the orthogonal arrangement of fibers in different layers, the angle distribution
between different layers is strictly defined.δ This ensures quasi-isotropic
0 mechanical behavior for the
i = arcos< di , di > (17)
components [3]. To analyze the angle deviation between the generated paths and the design datum,
an angle
Thedeviation
four setsanalysis was performed
of experiments selected in in
thisSection
section.4.1 were also adopted here to analyze the
The generated path was sampled
distribution of the angle deviation. uniformly with a step of 5 mm. For each sample point, the
forwardAs direction
mentioned of above,
the current path and
the initial paththe design
point was direction datumsampling
near the 350th of the path were
point. calculatedto
According
using both the guide-line and surface normal vector at the point. Subsequently,
Figure 13, for different triangulation parameters, the angle deviation does not depend on the the angle deviation
distance
was obtained.
between The forward
the sampling pointdirection
and the vector
initial pathatpoint
each but
sample
dependpoint
on isthe
the tangent of
curvature vector of the
the original
sample
[Link]
Theon the path.
larger The datum
the curvature ofdirection vector
the surface, wherecan the be calculated
sample points using:
are located on the path,
= ×
(16)
= ×
 The four sets of experiments selected in Section 4.1 were also adopted here to analyze the
distribution of the angle deviation.
As mentioned above, the initial path point was near the 350th sampling point. According to
Figure 13, for different triangulation parameters, the angle deviation does not depend on the distance
between
Materials 2020,the sampling point and the initial path point but depend on the curvature of the 15
13, 4209 original
of 18
surface. The larger the curvature of the surface, where the sample points are located on the path, the
greater is the angle deviation. This confirms that the angle deviation is due to the approximation of
thethe
greater
normalis the angle
vector ofdeviation. This mesh
the triangular confirms that the
surface angle
during deviation
the executionis due to the
of the approximation
algorithm, and thereof is
theno
normal vector
cumulative [Link] the triangular mesh surface during the execution of the algorithm, and there is no
cumulative error.

Figure 13. The angle deviation with different parameters: (a) Dl2s = 0.2, Lseg = 65; (b) Dl2s = 0.4,
Figure 13. The angle deviation with different parameters: (a) Dl2s = 0.2, Lseg = 65; (b) Dl2s = 0.4, Lseg =
Lseg = 110; (c) Dl2s = 0.6, Lseg = 155 and (d) Dl2s = 0.8, Lseg = 200.
110; (c) Dl2s = 0.6, Lseg = 155 and (d) Dl2s = 0.8, Lseg = 200.

To study the effect of Dl2s and Lseg on the angle deviation of the path, the mean angle deviation,
To study the effect of Dl2s and Lseg on the angle deviation of the path, the mean angle deviation,
the mean square error, and the maximum angle deviation were calculated using the experimentally
the mean square error, and the maximum angle deviation were calculated using the experimentally
obtained data—see Table 3.
obtained data—see Table 3.
Table 3. Data used for the calculation of the angle deviation.

Mean Angle Mean Square Maximum Angle

Dl2s Lseg
Deviation Error Deviation
/mm /mm
/deg /deg2 /deg
0.2 20 0.10162 0.09583 0.43632
0.2 65 0.19909 0.1562 1.019
0.2 110 0.21098 0.15039 1.01133
0.2 155 0.21159 0.15056 1.01142
0.2 200 0.21204 0.15059 1.01141
0.4 20 0.10162 0.09583 0.43632
0.4 65 0.22812 0.16134 1.1537
0.4 110 0.27022 0.21115 1.47317
0.4 155 0.27248 0.20807 1.47248
0.4 200 0.27248 0.20807 1.47248
0.6 20 0.10162 0.09583 0.43632
0.6 65 0.2696 0.24432 1.47683
0.6 110 0.33184 0.27956 1.41166
 Materials 2020, 13, 4209 16 of 18

Table 3. Cont.

Mean Angle Mean Square Maximum Angle

Dl2s Lseg
Deviation Error Deviation
/mm /mm
/deg /deg2 /deg
0.6 155 0.36801 0.27632 1.4119
0.6 200 0.36668 0.27738 1.4119
0.8 20 0.10162 0.09583 0.43632
0.8 65 0.27083 0.2446 1.47683
0.8 110 0.3184 0.25209 1.50431
0.8 155 0.39048 0.25162 1.50424
0.8 200 0.38589 0.2524 1.50424
1 20 0.10162 0.09583 0.43632
1 65 0.26882 0.24516 1.47683
1 110 0.33021 0.29158 2.0499
1 155 0.45625 0.32954 2.05242
1 200 0.44346 0.32784 2.05242
Materials 2020, 13, x FOR PEER REVIEW 17 of 19
The distribution of average angle deviation data is shown in Figure 14, and the maximum angle
deviation data distribution is shown in Figure 15.

Figure 14. Distribution of the average angle deviation.

Figure 14. Distribution of the average angle deviation.
According to Figures 14 and 15, as the Dl2s and Lseg increased, both the mesh density and
approximation accuracy of the triangular mesh surface decrease. In addition, the normal vector of the
path points on the triangular mesh surface deviate significantly from the normal vector on the original
parametric surface. Hence, the angle deviation increases. According to Table 3, when Dl2s exceeds
0.6 mm and Lseg exceeds 65 mm, the mean angle deviation surpasses 0.25 deg, and the maximum
deviation is more than 1.4 deg.
Based on the above distance deviation analysis, when Dl2s ranges between 0.2 and 0.6 mm and
Lseg is 65 to 110 mm, the mean distance deviation remains within 1mm. Furthermore, the maximum
distance deviation stays within 2 mm, the mean angle deviation is less than 0.25 deg, and the maximum
Materials 2020, 13, 4209 17 of 18

angle deviation is below 1.4 deg. At the same time, the path generation time stays within 0.5s. In this
case, Dl2s and Lseg can be Figure
selected
[Link] toofthe
Distribution thecomplexity of deviation.
average angle the surface and the acceptable path
error. Within the above range of the parameters, a high-precision triangular mesh surface and fiber
path with small error can be obtained, while the generation efficiency of the algorithm is high.

Figure 15. Distribution of the maximum angle deviation.

Figure 15. Distribution of the maximum angle deviation.
5. Conclusions
To improve the efficiency of automated fiber path planning process, a new path planning algorithm
based on meshing and multi guide-lines were investigated. The original parameter surface of the CAD
model of the FRP component was discretized into triangular mesh surface via surface discretization
and triangulation. Sub-surface boundary splicing and surface topology reconstruction algorithm
was proposed, and both the computational complexity reduction and the efficiency improvement
of the algorithm were analyzed. The proposed automated fiber path planning algorithm consists of
a main motion path direction vector algorithm and a continuous path point generation algorithm.
An updating method for the datum direction vector via the guide-lines update algorithm was also
introduced for complex surfaces. It improves the laying ability of the fibers and surface adaptability
for the planned path. Accuracy analysis was conducted to investigate the relationship between the
triangulation parameters and distance deviation, angle deviation and algorithm efficiency. The analysis
indicated that by choosing appropriate triangulation parameters, the fiber path can be generated with
high accuracy and efficiency.
More research efforts in the future work should be devoted to conduct experiments to test the
mechanical properties of the fabricated FRP components by using the multi-guide-line planned paths.

Supplementary Materials: The following is available online at [Link]

Supplementary Video S1 demonstrates the high efficiency of the proposed algorithm, which can complete the
path planning for one layer of a complex surface in just only tens of seconds.
Author Contributions: Conceptualization, H.X., W.H., W.T. and Y.D.; methodology, W.H. and W.T.; software,
H.X., W.H. and W.T.; validation, H.X. and W.H.; investigation, W.H.; writing—original draft preparation, H.X. and
Materials 2020, 13, 4209 18 of 18

W.H.; writing—review and editing, H.X. and W.H.; supervision, H.X.; project administration, Y.D. All authors
have read and agreed to the published version of the manuscript.
Funding: This research was funded by the National Natural Science Foundation of China (Grant No. 51875440),
China Postdoctoral Science Foundation (Grant No. 2019M663686) and the Open Fund of the State Key Laboratory
for Manufacturing Systems Engineering (Grant No. sklms2020003).
Conflicts of Interest: The authors declare no conflict of interest.

References
1. Kozaczuk, K. Automated fiber placement systems overview. Trans. Inst. Aviat. 2016, 245, 52–59. [CrossRef]
2. August, Z.; Ostrander, G.; Michasiow, J.; Hauber, D. Recent developments in automated fiber placement of
thermoplastic composites. SAMPE J. 2014, 50, 30–37.
3. Rousseau, G.; Wehbe, R.; Halbritter, J.; Harik, R. Automated Fiber Placement Path Planning: A state-of-the-art
review. Comput.-Aided Des. Appl. 2018, 16, 172–203. [CrossRef]
4. Lewis, H.; Romero, J. Composite Tape Placement Apparatus with Natural Path Generation Means. U.S.
Patent 4,696,707, 29 September 1987.
5. Shirinzadeh, B.; Foong, C.W.; Tan, B.H. Robotic fibre placement process planning and control. Assem. Autom.
2000, 20, 313–320. [CrossRef]
6. Shirinzadeh, B.; Alici, G.; Foong, C.W.; Cassidy, G. Fabrication process of open surfaces by robotic fibre
placement. Robot. Comput.-Integr. Manuf. 2004, 20, 17–28. [CrossRef]
7. Shirinzadeh, B.; Cassidy, G.; Oetomo, D.; Alici, G.; Ang, M.H. Trajectory generation for open-contoured
structures in robotic fibre placement. Robot. Comput.-Integr. Manuf. 2007, 23, 380–394. [CrossRef]
8. Peng, Z.; Ronglei, S.; Xueying, Z.; Lingjin, H. Placement suitability criteria of composite tape for mould
surface in automated tape placement. Chin. J. Aeronaut. 2015, 28, 1574–1581.
9. Zhang, P.; Sun, R.; Huang, T. A geometric method for computation of geodesic on parametric surfaces.
Comput. Aided Geom. Des. 2015, 38, 24–37. [CrossRef]
10. Savio, G.; Meneghello, R.; Concheri, G. Geometric modeling of lattice structures for additive manufacturing.
Rapid Prototyp. J. 2018, 24, 351–360. [CrossRef]
11. Zhang, Q.; Sabin, M.A.; Cirak, F. Subdivision surfaces with isogeometric analysis adapted refinement weights.
Comput.-Aided Des. 2018, 102, 104–114. [CrossRef]
12. Shinno, N.; Shigemat, T. Method for Controlling Tape Affixing Direction of Automatic Tape Affixing
Apparatus. U.S. Patent 5,041,179, 20 August 1991.
13. Li, L.; Wang, X.; Xu, D.; Tan, M. A Placement Path Planning Algorithm Based on Meshed Triangles for Carbon
Fiber Reinforce Composite Component with Revolved Shape. Int. J. Control Syst. Appl. 2014, 1, 23–32.
14. Shen, J.; Buse, L.; Alliez, P.; Dodgson, N.A. A line/trimmed NURBS surface intersection algorithm using
matrix representations. Comput. Aided Geom. Des. 2016, 48, 1–16. [CrossRef]
15. Lo, S.H.; Wang, W.X. An algorithm for the intersection of quadrilateral surfaces by tracing of neighbours.
Comput. Methods Appl. Mech. Eng. 2003, 192, 2319–2338. [CrossRef]

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access
article distributed under the terms and conditions of the Creative Commons Attribution
(CC BY) license ([Link]