Feature extraction and process planning of integrated hybrid additive-subtractive system for remanufacturing

1.   Introduction

Nowadays, the increasing developments and over-exploitation of resources outcome vast amount "end-of-life" products. Remanufacturing and repairing are identified as environmentally friendly approaches to deal with them since the energy consumption to produce a new part can be minimized by maintaining the intrinsic energy of the legacy part ^[1,2,3]. Repairing process is the correction of specified faults or restoring its original form of product. Whereas, remanufacturing, also called "reincarnation", is defined as being able to produce new products directly from the end-of-life product ^[4]. In different publications, the words "remanufacturing" and "repairing" have been used interchangeably with blurred definitions ^[4,5,6]. In this study, the authors differentiate them by defining "repairing" as restoring the part to its original form and "remanufacturing" as upgrades the part to have different functional features.

Recently, hybrid manufacturing (HM) system synergistically integrates additive manufacturing (AM) and subtractive manufacturing (SM) processes within a single workstation (Figure 1), has gained a lot of attention from academia and industry ^[7]. HM can capitalize on the strengths of independent techniques, whilst minimizing their disadvantages ^[8]. At the same time, it has the potentials to enable remanufacturing technologies to achieve further improvement because it can to remove and add features flexibly.

Most of the single platform workstations for HM processes are integrating subtractive CNC machining and additive directed energy deposition (DED) ^[9]. There are some commercial process planning software tools that have already been developed for industrial hybrid systems, such as Siemens NX used in DMG MoriSeiki, LaserTec 65 3D, and hyperMill used in Replicator and Cybaman ^[9]. However, the major limitation of these commercial tools is the poor support for automation. For example, AM and SM feature recognition is manually conducted by the user's knowledge. The consequence is that the process planning of complex remanufacturing problems is very tedious, and the quality of the derived process plan can hardly be evaluated. For this reason, an increasing level of interest in research on process planning for HM has been witnessed over recent years ^[9,10]. Joshi & Anand ^[11] developed a novel metric, called complexity score to quantify part's geometric complexities for decision-making among AM, SM, and HM. Chen et al. ^[12] discovered that complex geometry might have the tool accessibility issue for the SM process, in which the cutter cannot access the part's interior. They then developed an optimization algorithm for HM process planning, which derived the optimal HM process plan of an arbitrary geometry that was free of machinability issues. They also extended the concept to address the problem for quasi-rotational parts with columnar features in their other works ^[7]. A sequence planning method for five-axis additive-subtractive hybrid manufacturing is also presented by ^[13], and it greatly reduced the number of tool changes between different operations. Zhu et al. ^[14,15] and Newman et al. ^[16] presented a hybrid process planning framework, called iAtractive, which integrates CNC machining, fused filament fabrication, and inspection process for plastic components with internal structures. Behandish et al. ^[17] introduced an HM process planning method that formulated HM processes through logic representations, which has the capability to enumerate all the feasible sequences of operations to find the optimal plan. ElMaraghy & Moussa ^[18] developed a method of process planning for HM that included manual feature extraction, product platform design, and process determination. Basinger et al. ^[19] outlines a feature-based advanced HM process planning system which use feature, tolerance and material data as inputs to construct process plans. From these research results, the authors identify that SM is mostly playing a role as post-machining for AM in the HM process. Therefore, it is complicated to be applied directly in the remanufacturing, since SM is not only post-machining but also including geometric forming as AM process in remanufacturing.

Several studies have explored the potentials of AM for repairing or remanufacturing applications. Different aspects of AM techniques have been widely investigated, such as laser parameters ^[20], damaged volume reconstruction ^{[21,22,23,24]}, energy and environmental impact analysis ^[6,25], and microstructure ^[26,27,28] and mechanical property analysis ^[29,30]. Recently, HM for repairing application is also actively investigated. Zheng et al. ^[22] introduced a 3D reconstruction and triangle intersection algorithm to identify additive repair and subtractive repair areas of a broken turbine surface. Hascoët et al. ^[31] proposed a method to automate the repairing process of metallic parts partially. In their practice, defects were initially machined into a surface cavity, and the cavity was refilled by laser metal deposition with the aid of an inspection system. An increasing level of development in research on HM-based repairing technology has been witnessed by those publications. However, in terms of remanufacturing, there are comparatively fewer studies. In comparison with repairing, remanufacturing requires more decision-making support ^[32] since the process planning result is not unique, and AM and SM feature extraction relies on algorithms for the automation process.

Newman et al. ^[16] and Zhu et al. ^[33] features a remanufacturing framework that consisted of fused filament fabrication, CNC machining, and inspection. The method enabled the remanufacturing of an existing part or even a recycled and legacy part into a new part with new functional features. Le et al. ^[34] proposed an HM process planning for remanufacturing based on feature extraction and knowledge interpretation. In a subsequent study ^[35], they have extended the process planning framework by discussing the environmental impact of the proposed remanufacturing strategy compared with traditional approaches (material recycling, casting, and machining). Liu et al. ^[3] developed a novel design-for-remanufacturing method under a level-set framework, which provides a solution for upgrading broken parts. A cost-driven process planning for PBF-CNC remanufacturing was proposed in the algorithmic framework laid down in our previous work ^[1]. In that work, an automated additive-subtractive feature extraction method is developed and process planning sequencing is formulated as a cost-minimization optimization problem. However, there are three significant limitations: 1). this work is focusing on primitives features and free-form features are not covered; 2). the inspection process for the used part or legacy part is not included; 3). commonly, there are defects or damaged area on the used part, the process for dealing with the defects should be investigated.

To sum up, the research gaps and objectives of this study is the development of a method that provides automated feature extraction and cost-driven process planning for an integrated HM machine. The main contributions are list as follows:

● Both primitive and free-form features are modelled in level set-based representations for the automated feature extraction which facilitate the process planning for HM remanufacturing;

● A collision-free DED-CNC process planning method is developed, resulting in the minimal cost in HM remanufacturing process;

● The defect and damaged area of the used part are considered to form a pre-machining feature in the process planning, which is an issue rarely addressed in the previous studies of process planning for remanufacturing.

The rest of the paper is structured as follows. In Section 2, the automated additive and subtractive feature extraction method is presented, which includes level set function representation for CAD and point clouds, extraction for pre-machining feature, intersection part extraction, intersection part modification for collision-free remanufacturing and individual feature extraction. Section 3 shows the cost-driven process planning methodology for DED-CNC remanufacturing, which includes cost estimation for DED-CNC remanufacturing, and sequence optimization. Section 4 demonstrates a case study to validate the efficacy of the proposed method. Finally, a conclusion is given in Section 5.

2.   Additive and subtractive feature extraction

In geometric modelling, constructive solid geometry (CSG) and boundary representation (B-rep) are widely adopted. With CSG modelling, a physical object can be decomposed into multiple primitives and a sequence of Boolean operations. With the B-rep method, the solid is bounded by a set of closed and directional faces, which are bounded by edges and vertices. CSG modelling has the merits of supporting efficient Boolean operations and topology optimization ^[36]. Previous works on feature recognition/extraction in machining process planning have been developed for three decades ^[37,38]. However, these methods cannot be applied in remanufacturing or HM process. The reason is that most of the current machining feature extraction methods uses boundary representation (B-rep) format for feature modeling because it uniquely defines the faces and their topological patterns ^[38]. However, it has issues to deal with numerical calculations between two solid parts by adopting B-rep models for extracting features between the used part and final part. Therefore, in this work, a novel feature extraction method based on the level set function that implicitly represents models is developed for hybrid manufacturing in remanufacturing context.

The overall framework of the level set-based feature extraction method is presented in Figure 2. First the modeling history is extracted from the CAD model of the final part to be modelled in level set function representation (Section 2.1). The point cloud data is acquired from the used part and converted to be the level-set function representation (Section 2.2). Then, the pre-machining feature is extracted from the used part (Section 2.3). The intersection volume is then calculated by optimally overlapping the used part and the final part (Section 2.4). The intersection part is then modified with considering the collision problem and DED manufacturing constraints (Section 2.5). With The level set represented intersection part. The subtractive manufacturing volume (SMV) and additive manufacturing volume (AMV) can be identified and the individual AFs and SFs are extracted from the volumes (Section 2.6).

2.1.   Level set function representation for the CAD model

Level set function $\mathrm{\Phi }\left({\bf{X}}\right)({R}^{n}\to R)$ describes the geometry in an implicit form, as shown in Eq (1).

where $\mathrm{\Omega }/\partial \mathrm{\Omega }$ is the material domain, $D/\mathrm{\Omega }$ is the void, $\partial \mathrm{\Omega }$ is the structural boundary.

In the level set function-based modelling approach, the 3D model is constructed by bounding the boundary surfaces, as:

As an example, the cube with ${(x}_{0}, {y}_{0}, {z}_{0})$ as the center coordinates and $(Hx, Hy, Hz)$ as the lengths on the x, y, z axis can be represented by bounding six planer surfaces of ${\mathrm{\Phi }}_{1} = \frac{Hx}{2}-\left(x-{x}_{0}\right) = 0, {\mathrm{\Phi }}_{2} = \frac{Hx}{2}+\left(x-{x}_{0}\right) = 0, {\mathrm{\Phi }}_{3} = \frac{Hy}{2}-\left(y-{y}_{0}\right) = 0, {\mathrm{\Phi }}_{4} = \frac{Hy}{2}+\left(y-{y}_{0}\right) = 0, {\mathrm{\Phi }}_{5} = \frac{Hz}{2}-\left(z-{z}_{0}\right) = 0, {\mathrm{\Phi }}_{6} = \frac{Hz}{2}+\left(z-{z}_{0}\right) = 0$, as shown in Figure 3.

In our previous work ^[1], the examples of primitive geometry (cube, sphere, cone and cylinder) are given. This work will focus on the level set function modelling for freeform geometries. From Eq (2), it can be manifested that finding the implicit forms for boundary surfaces is the most crucial work for the level set-based modelling approach. Algebraic techniques based on elimination theory enable the conversion of parametric expression to its implicit expression $\mathrm{\Phi }\left({\bf{X}}\right) = 0$. Elimination theory investigates the conditions under that the sets of parametric expressions have common roots. The vanishing of the resultant is a necessary and sufficient condition for the parametric expressions to have a common non-trivial root. The implicitization of parametric geometry is based on the construction of these resultants.

A freeform 2.5D geometry can be constructed by extruding its freeform profile. As an example, in this study, the Bezier curve is implemented to represent the freeform profile. The parametric form of the Bezier curve is shown as:

where $\mathit{\boldsymbol{F}}\left(u\right) = \left[{f}_{x}\left(u\right), {f}_{y}\left(u\right)\right]$, ${\bar {\mathit{\boldsymbol{P}}}}_{i} = \left[{\bar {p}}_{i}^{x}, {\bar {p}}_{i}^{y}\right]$, ${B}_{i, n}\left(u\right) = C\left(n, i\right){u}^{i}{\left(1-u\right)}^{n-i}$, $C\left(n, i\right)$ is the binomial coefficient: $C\left(n, i\right) = \frac{n!}{i!\left(n-i\right)!}$, n is the degree of the curve, and i is the number of control points, ${\bar {P}}_{ix}$ and ${\bar {P}}_{iy}$are x and y coordinates of the control point.

By following elimination theory, the implicit form of the Bezier curve ${\mathrm{\Phi }}_{\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{v}\mathrm{e}}\left({\bf{X}}\right)$ can be obtained by eliminating the parameter u between the parametric expressions in Eq (3) letting the resultant of them to be equal to zero.

As an example, a cubic Bezier curve is constructed by three control points: ${\bar {P}}_{0} = \left(\mathrm{0, 0}\right), {\bar {P}}_{1} = \left(\mathrm{40,220}\right), {\bar {P}}_{2} = \left(\mathrm{200, 40}\right)$ and ${\bar {P}}_{3} = \left(\mathrm{0, 0}\right)$. By implementing the elimination theory, the implicit form for the Bezier curve can be obtained and the contour figure is shown in Figure 4a. Then, the 2.5D freeform shape can be modelled by $\mathrm{\Phi }\left({\bf{X}}\right) = \mathrm{m}\mathrm{i}\mathrm{n}({\mathrm{\Phi }}_{\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{v}\mathrm{e}}, {\mathrm{\Phi }}_{\mathrm{t}\mathrm{o}\mathrm{p}}, {\mathrm{\Phi }}_{\mathrm{b}\mathrm{o}\mathrm{t}\mathrm{t}\mathrm{o}\mathrm{m}})$, where ${\mathrm{\Phi }}_{\mathrm{t}\mathrm{o}\mathrm{p}} = 157.5-z$ and ${\mathrm{\Phi }}_{\mathrm{b}\mathrm{o}\mathrm{t}\mathrm{t}\mathrm{o}\mathrm{m}} = z-122.5$. (see Figure 4b.)

The parametric Bezier surface is extended from Bezier curve to u and v directions, as follows:

where $\mathit{\boldsymbol{F}}\left(u, v\right) = \left[{f}_{x}\left(u, v\right), {f}_{y}\left(u, v\right), {f}_{z}\left(u, v\right)\right]$ and ${\bar {\mathit{\boldsymbol{P}}}}_{i, j} = [{\bar {p}}_{i, j}^{x}, {\bar {p}}_{i, j}^{y}, {\bar {p}}_{i, j}^{z}]$.

Analogously, the elimination theory can help to find the implicit form for the parametric expression in Eq (4). A Bezier surface is modelled by 3 × 3 control points; the level set form of this surface is shown in Figure 5a. The Bezier shape is bounded with planar surfaces ${\mathrm{\Phi }}_{1} = x-200, {\mathrm{\Phi }}_{2} = 600-x, {\mathrm{\Phi }}_{3} = y-200, {\mathrm{\Phi }}_{4} = 600-y$, see Figure 5b.

The complex geometry can be constructed through Boolean operations on the level-set functions ^[39]. However, this representation will cause non-differentiable problems in numerical calculation. R-functions can combine level functions of a complex structure into a new smooth level set function by operations of R-conjunction $\bigwedge$ and R-disjunction $\bigvee$, which are equivalent to Boolean operations $\cap$ and $\cup$ ^[40]. The operations of R-functions are defined as:

2.2.   Level set function representation for point clouds

For a given used part, data acquisition is the first step to digitalize the part to point clouds. Therefore, the other situation is that we have the point cloud as the geometry input. To converting the point cloud data to the level set function model, the surface fitting techniques can be applied to obtain the parameters of surfaces and these parameters are used to forming the level set function representations. In this study, random sample consensus (RANSAC) surface fitting ^[41] is employed for surface fitting. The pseudocode for the forming level set function representation of the point clouds is given in Table 1. Figure 6 provides two examples of the surface fitting for cylindrical surfaces and sphere from point clouds from our previous results ^[21]. The color scale bars indicate the distance of each point to the fitted surface. The parameters of fitted surfaces are utilized for forming the level set functions through the proposed algorithm.

2.3.   Extraction of pre-machining feature

In the repairing process, the defects on the damaged part need to be machined into a surface cavity. This cavity is to be refilled by the deposition of materials to recover the local geometry of the part. Equivalently, in terms of the remanufacturing process, it is necessary to carve out defects from the damaged part to eliminate the perturbation for the remanufacturing process planning caused by the defects. In this study, the authors define the machining feature to carve out defects as a pre-machining feature since it has a similar concept as the pre-machining in the traditional machining process to remove the imperfections of the stock. It is important to mention that this study is focusing on remanufacturing end-of-life part with small defects which do not cover the any entire features on the end-of-life part.

To construct the pre-machining feature, the first step is to segment surface defects from 3D scan data. Many methods for surface defect segmentation have been presented in the literature ^[21,42,43]. In this work, random walks for unorganized point cloud segmentation ^[42] is adopted, since it does not rely on strong assumptions made on the characteristics of the expected defect or the geometry of the surrounding area. The algorithm segments defect areas on a weighted, undirected k-nearest neighbour graph (k-NNG) defined by local changes in point cloud properties. Since it is not the contribution of the authors' work, the detailed algorithm is not demonstrated. The interested reader can refer to the original work ^[42]. An example of the defect segmentation from point clouds by the random walks algorithm is shown in Figure 7. In the figure, the blue points indicate the points without defects and red points refer the points on the defective area.

Different machining features can be applied to carve out the defect area. In this study, the authors only investigate two basic machining features for the sake of simplicity. It is trivial to extend the method for other machining features. The hole and rectangular pocket are shown in Figure 8, and the level set function to represent these features are presented in Eq (6).

It is worth investigating the type and parameters of the machining feature, which leads to minimal materials being carved out. Meanwhile, the machining feature must remove all the defects. The problem can be mathematically formulated as a constrained optimization problem. For a given machining feature, the optimization problem is solving the optimal parameters. The objective function is minimizing the volume of the given machining feature by integrating $\mathrm{d}x\mathrm{d}y\mathrm{d}z$ in the material domain, referred to Eq (7).

where $\mathrm{H}\left(\right)$ is the Heaviside function, $\mathrm{\Phi }\left(\right)$ is the level set function for the given machining feature, ${\bf{X}} = (x, y, z)$, ${\bf{a}} = \left({a}_{1}, \dots, {a}_{n}\right)$ which indicates the parameters for the given machining feature.

To satisfy the condition that all defects are removed, the defective points need to be enclosed by the machining feature. Assuming there are defect points $({\mathit{\boldsymbol{p}}}_{1}, {\mathit{\boldsymbol{p}}}_{2}, \dots, {\mathit{\boldsymbol{p}}}_{M})$, M series of constraints are formulated as in Eq (8) and the schematic plot is presented in Figure 9.

Lagrange formulation of this problem is written as:

where ${\lambda }_{i}$ is the i-th Lagrange multiplier and ${\eta }_{i}$ is the i-th slack variable.

The sensitivity $\frac{\partial \mathcal{L}}{\partial {a}_{1}}$ of the Lagrange formulation is derived through Eq (10); others can be calculated similarly.

where $\delta$ is the Dirac delta function.

This problem can be solved with a gradient-based optimization solver by updating the variables ${a}_{1}, \dots, {a}_{n}, {\lambda }_{1}, \dots, {\lambda }_{n}, {\eta }_{1}, \dots, {\eta }_{n}$ with their corresponding sensitivity.

Examples of the optimal parameters of a hole and a rectangular pocket features are illustrated in Figure 10. From the two pre-machining features, the minimal volume of the feature can be determined. In this example, the hole feature has a volume of $3.657\;{\mathrm{m}\mathrm{m}}^{3}$ and rectangular pocket has the volume of $4.3665\;{\mathrm{m}\mathrm{m}}^{3}$. Therefore, the hole feature is optimal for pre-machining.

2.4.   Intersection part extraction

As the level set function of the used part and final part are formulated, the relative position needs to be identified between two solid models to prepare for feature extraction. The objective is to maximize the overlapping material volume because the cost of AM is strongly affected by the volume of the new material to deposit.

${\mathrm{\Phi }}_{u}$ and ${\mathrm{\Phi }}_{f}$ are the level set functions of the used part and the final part, respectively, within the global coordinate system ${\bf{X}} = (x, y, z)$. $\tilde {{\bf{X}}}$ is a local coordinate system attached to the used part. Through Eq (11), an optimization problem can be formulated to figure out the translation and rotation of the local coordinate system ($\tilde {{\bf{X}}}$) needed to maximize the overlapping volume between ${\mathrm{\Phi }}_{u}$ and ${\mathrm{\Phi }}_{f}$.

where ${{\bf{R}}}_{\mathit{\boldsymbol{x}}}$, ${{\bf{R}}}_{\mathit{\boldsymbol{y}}}, {{\bf{R}}}_{\mathit{\boldsymbol{z}}}$ are the rotation matrix along x, y, z direction with the variables of ${\theta }_{x}, {\theta }_{z}, {\theta }_{y}$, ${\bf{T}}$ is the translation matrix with the variables of ${t}_{x}, {t}_{y}{, t}_{z}$.

The used part ${\mathrm{\Phi }}_{i}({\theta }_{x}, {\theta }_{z}, {\theta }_{y}, {t}_{x}, {t}_{y}{, t}_{z})$ has translation variables and rotation variables. The final part is fixed by ${\mathrm{\Phi }}_{f}$. The intersection part (${\mathrm{\Phi }}_{i}$) is the intersection of the used part and the final part, given as:

The maximization optimization problem can be formulated of maximizing the intersection part by optimizing $\mathrm{t}\mathrm{h}\mathrm{e}\;\mathrm{ }\mathrm{v}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{s}\;{\theta }_{x}, {\theta }_{z}, {\theta }_{y}, {t}_{x}, {t}_{y}{, t}_{z}$, through in Eq (13):

The sensitivity $\frac{\partial f}{\partial {\mathrm{t}}_{x}}$ of the objective function is derived through Eq (14); others can be calculated similarly.

This problem can be solved with a gradient-based optimization solver. In the gradient-based optimization algorithm, the $({\theta }_{x}, {\theta }_{z}, {\theta }_{y}, {t}_{x}, {t}_{y}{, t}_{z})$ = (0, 0, 0, 0, 0, 0) is set as the initial values. Generally, to prevent a local optimum issue, a multi-start strategy is suggested, i.e., to parallelly run the optimization program with a different initial guess of variables.

Figure 11 demonstrates an example of the intersection part extraction. The used part and the final part are represented by its level set function ${\mathrm{\Phi }}_{u}$ and ${\mathrm{\Phi }}_{f}$ respectively. Eq (13) helps to find the optimal transformation to transform ${\mathrm{\Phi }}_{u}$ to ${\tilde {\mathrm{\Phi }}}_{u}$ with the optimized rotation angle (${\theta }_{x} = -9{0}^{o}, {\theta }_{y} = {0}^{o}, {\theta }_{z} = {0}^{o}$) and translation (${t}_{x} = -20.1, {t}_{y} = -110.7, {t}_{z} = 15$). The intersection part ${\mathrm{\Phi }}_{i}$ can be obtained by ${\tilde {\mathrm{\Phi }}}_{u}\cap {\mathrm{\Phi }}_{f}$.

2.5.   Intersection part modification for collision-free remanufacturing

The geometry of intersection part geometry is not generally acceptable for AM processes. Figure 12 provides an illustration of collision problems in the DED process, and the material deposition nozzles may have collisions with the intersection part. Therefore, it is crucial to modify the intersection part by analyzing the tool accessibility constraints of the DED process.

Initially, the intersection part and DED nozzle are formulated as ${\mathrm{\Phi }}_{i}$ and ${\mathrm{\Phi }}_{n}$ by the level-set representation, respectively. In addition, the deposited material volume by AM is represented as ${{\mathrm{\Phi }}_{d} = \mathrm{\Phi }}_{f}/{\mathrm{\Phi }}_{i}$, as shown in Figure 13. The voxel representations for ${\mathrm{\Phi }}_{i}$, ${\mathrm{\Phi }}_{n}$ and ${\mathrm{\Phi }}_{d}$ are computing by applying Heaviside functions ${V}_{i} = \mathrm{H}\left({\mathrm{\Phi }}_{i}\right)$, ${V}_{n} = \mathrm{H}\left({\mathrm{\Phi }}_{i}\right)$, and ${V}_{d} = \mathrm{H}\left({\mathrm{\Phi }}_{d}\right)$, separately.

As far as the collision detection is concerned, the collision-free rigid motion of the DED nozzle in rigid motion needs to be calculated. In this study, dilation as a morphology operation is adopted to analyze the spatial planning of the DED nozzle. As one of the basic operations in mathematical morphology, dilation operation $\oplus$ is defined as:

where ${A}_{b}$ represents the solid $A$ transformed by a rigid transformation $b$, $B$ is a structuring element, which is termed filters. For level set function represented models, the dilation operation is expressed as:

For the DED process, the materials are deposited on the deposition volume layer by layer, which indicates that the tip of the deposition nozzle requires going through each point of the deposition volume (Figure 14a). It is crucial to mention that in a practical case, the vertical distance between the tip of the DED nozzle and the deposition area is not zero. The vertical distance is assumed to be zero for the sake of simplifying in explaining the proposed method. The technical implantation of the approach will be discussed with considering the vertical distance later in this section. It is meaningful to explore all motions of the DED nozzle to deposit the materials in the deposition volume. In the proposed method, the motions of the DED nozzle during deposition is calculated by dilating the deposition volume by the DED nozzle: ${\mathrm{\Phi }}_{m} = {\mathrm{\Phi }}_{n}\oplus {\mathrm{\Phi }}_{d}$, and the authors define ${\mathrm{\Phi }}_{m}$ as motion space of the DED nozzle, as presented in Figure 14b.

Theorem:

The modified intersection part ${\tilde {\mathrm{\Phi }}}_{i}$ can be derived from: ${\tilde {\mathrm{\Phi }}}_{i} = {\mathrm{\Phi }}_{i}\setminus {\mathrm{\Phi }}_{m} = {\mathrm{\Phi }}_{i}\setminus ({\mathrm{\Phi }}_{n}\oplus {\mathrm{\Phi }}_{d})$, which leads no collision with the DED nozzle in operation, as shown in Figure 14d.

Proof:

In the condition of no collision occurring during deposition, ${\tilde {\mathrm{\Phi }}}_{i}\cap {\tilde {\mathrm{\Phi }}}_{m} = \mathrm{\varnothing }$ is required to be satisfied, where ${\tilde {\mathrm{\Phi }}}_{m}$ is the new motion space derived from the modified intersection part ${\tilde {\mathrm{\Phi }}}_{m} = {{\mathrm{\Phi }}_{n}\oplus \tilde {\mathrm{\Phi }}}_{d}$. ${\tilde {\mathrm{\Phi }}}_{m}$ can be expanded into:

where ${\mathrm{\Phi }}_{f}$ indicates the level set function of the final part. According to the associativity, commutativity of Boolean operation and distribution of dilation, Eq (17) can be rearranged as:

Therefore, no collusion condition is derived as:

Since ${\mathrm{\Phi }}_{n}\oplus {\mathrm{\Phi }}_{d} = {\mathrm{\Phi }}_{m}$, $\left(-{\mathrm{\Phi }}_{m}\right)\cap \left({\mathrm{\Phi }}_{n}\oplus {\mathrm{\Phi }}_{d}\right) = \mathrm{\varnothing }$. So, ${\tilde {\mathrm{\Phi }}}_{i}\cap {\tilde {\mathrm{\Phi }}}_{m} = {\mathrm{\Phi }}_{i}\cap \mathrm{\varnothing }\cap \left({\mathrm{\Phi }}_{n}\oplus {\mathrm{\Phi }}_{m}\right).$ Due to the annihilator law for $\cap$, ${\tilde {\mathrm{\Phi }}}_{i}\cap {\tilde {\mathrm{\Phi }}}_{m} = \mathrm{\varnothing }$.

The authors develop an algorithm to implement the proposed method for discrete level set function representations, also considering the vertical distance (d) between the nozzle tip and deposition area. The pseudocode for the proposed algorithm is presented in Table 2. An example of the implantation of the proposed method is shown in Figure 15. Figure 15a presents the original intersection part and Figure 15b shows deposition volume derived from the original intersection part by ${\mathrm{\Phi }}_{f}\backslash {\mathrm{\Phi }}_{i}$. The motion space calculated from Eq (16) is shown in Figure 15c, and Figure 15d gives the modified intersection part ${\tilde {\mathrm{\Phi }}}_{i}$.

2.6.   Individual feature extraction

As the modified intersection part has been calculated, a subtractive manufacturing volume (SMV) and an additive manufacturing volume (AMV) can be calculated in their level set function representations as:

Individual features are required to be recognized and extracted from SMV and AMV. The research of CSG-based feature recognition has been developed from the 1990s, but this technique did not go far primarily due to the non-uniqueness of CSG trees ^[38]. Recently, there are some research efforts have been devoted to address non-uniqueness problems in CSG and show the strength in recognize sophisticated machining features by 3D convolution neuron network ^[44,45]. Therefore, developing a new machining feature recognition method is not a contribution to this study. The level set function representation of SFV is converted into 3D voxel grid information ${\mathrm{V}}_{\mathrm{S}\mathrm{F}\mathrm{G}}$ by Heaviside function$\mathrm{H}$, and then subtractive features (SFs) can be extracted by any CSG-based feature recognition techniques; see Eq (21).

AFV comprises both additive features (AFs) and SFs. In the AM process, leaving a sufficient over-thickness to have a finishing operation is vital for meeting the tolerance and surface roughness requirements. The over-thickness value is estimated by the required specifications of the final feature, the surface roughness generated by the AM processes and the machining conditions ^[34]. With considering over-thickness, the modified additive feature volume ($\overline{\mathrm{A}\mathrm{F}\mathrm{V}}$) is modified. The level set function representation of $\overline{\mathrm{A}\mathrm{F}\mathrm{V}}$ can be derived via Eq (22), where $\mathrm{t}$ represents the over-thickness value. The residual subtractive feature volume is obtained via Eq (23). Similarly, the individual SFs could be recognized by a CSG-based feature recognition method by Eq (24).

3.   Cost-driven process planning for remanufacturing

A large number of AFs and SFs are obtained from the proposed feature extraction method. Each feature represents an operation in the AM or SM process. Under considering the topological relationship of different features, the sequences of some operations are forced. Therefore, precedence constraints between operations are required to formulated to respect the hybrid additive-subtractive manufacturing rules. The precedence constraints are comprehensively summarized, and interested readers can refer to ^[1,34].

Although precedence constraints are formulated, some residual process sequences are still undermined. In order to develop the process planning problem into a process sequence optimization, the cost for each operation/feature is required to be estimated. With precedence constraints between features and cost model for each feature, an integer programming model is formulated to calculate the optimal process plan that minimizes the overall remanufacturing cost.

3.1.   Cost estimation for DED and CNC machining process

The motivation of the cost estimation for this study is approximating the cost models for the integer programming model to determine the optimal process plan, rather than exploring the precise cost for each operation. Although there are numbers of publications are relevant to cost estimation for the subtractive and additive manufacturing process ^[46,47,48], few researchers focus on the constructing cost model for hybrid additive-subtractive manufacturing system, explicitly considering the change cost between two different operations. Therefore, in this section, rough and fast estimations are given for SFs and AFs, and the change costs between various features are also introduced.

3.1.1.   Cost of SF

The total cost comprises operation cost and tool cost as:

where ${\mathrm{C}}_{\mathrm{S}\mathrm{F}}$ is the total cost of an SF, ${\mathrm{C}}_{\mathrm{h}\mathrm{h}}$ is the hourly operation cost for HM machine operation, ${\mathrm{t}}_{\mathrm{S}\mathrm{F}}$ is the machining time for the SF, ${\mathrm{C}}_{\mathrm{t}\mathrm{o}\mathrm{o}\mathrm{l}}$ is the cost of each cutting tool, and ${\mathrm{n}}_{\mathrm{t}\mathrm{o}\mathrm{o}\mathrm{l}}$ is the number of tool changes. Most commercial CAM software systems can estimate the machining time ${\mathrm{t}}_{\mathrm{S}\mathrm{F}}$ by dividing the tool path in the milling process by the programmed feed rate.

3.1.2.   Cost of AF

The cost of AF is determined by the machine cost and material consumption cost, see Eq (26).

where ${\mathrm{C}}_{\mathrm{A}\mathrm{F}}$ is the total cost of an AF, ${\mathrm{C}}_{\mathrm{h}\mathrm{h}}$ is the hourly operation cost for HM machine operation, ${\mathrm{t}}_{\mathrm{A}\mathrm{F}}$ is the building time of the AF, and ${\mathrm{C}}_{\mathrm{m}-\mathrm{A}\mathrm{F}}$ is the cost of material consumption.

The cost of material consumption in DED process is approximately calculated from the volume of the feature and its support structure, as:

where, ${\mathrm{V}}_{\mathrm{A}\mathrm{F}}$ and ${V}_{\mathrm{A}\mathrm{F}-\mathrm{s}}$ refer to the volume of the building part and support structure, $\rho$ is the material density, $\overline{\rho }$ is the material density of the support structure, and ${\mathrm{C}}_{\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}^{\mathrm{u}\mathrm{n}\mathrm{i}\mathrm{t}}$ indicates the price per unit of material. In DED process, complex gas flow leads a diffusion of powder distribution, which results the low powder efficiency because some powder cannot reach to the meltpool. Powder efficiency rate $\mu$ varies between different machines, and in this research the value is referred to ^[49,50], as 70%.

For the DED process, the building time ${\mathrm{t}}_{\mathrm{A}\mathrm{F}}$ estimation adopted an analytical build time model which is proposed in ^[51], and the general equation is:

This model uses G-code of the part as input and an algorithm extracts the kinematic characteristics of the nozzle to estimate very accurate build time results, since the acceleration and deceleration of the machine head is considered.

3.1.3.   Cost of change

For the hybrid additive-subtractive system, the AM/SM operations are switching by changing different tools. It is crucial to discuss the cost of change between two consecutive operations since it is costly due to frequent tool changing. Besides the tool change cost, the re-orientation of the workpiece also results in costs. In this study, CR and CT represent indexes of the re-orientation cost, and the tool change cost, respectively. The details of the calculation are given below.

1) Re-orientation cost:

While the orientation of the workpiece is switched in the HM machine, a re-orientation change occurs that requires workpiece fixing and laser calibration.

For the index of feature $\mathrm{i}, \mathrm{j}\in \mathrm{F} = \{\mathrm{1, 2}, ..., \mathrm{F}\}$, the index of cost for orientation $\mathrm{C}\mathrm{R}$ is formed by grouping the cost of orientation between any two features.

2) Tool change cost:

For HM machine the tool change cost is formulated as:

Similarly, the index of tool change cost $\mathrm{C}\mathrm{T}$ is constructed by grouping the cost of tool change cost between any two features.

3.2.   Sequence optimization

As the cost model for each feature, change cost, and precedence constraints are obtained. The process sequence optimization problem is formulated as an integer programming model and solved by branch and fathoming algorithm ^[52]. In a directed graph visualization of the model, a node represents a feature/operation and the cost model for them is calculated from Section 4.1.1 and 4.1.2, and the directed line indicates the precedence relation between two features/operations and the change costs are derived from Section 4.1.3. As an example, AFs (AF1 and AF2) and SFs (SF1, SF2, SF3) are extracted from the proposed feature recognition method. The cost for each node is calculated and the change cost between any two nodes are formulated. The precedence constraints are applied to these nodes initially as dash lines (Figure 16a). With the optimization of the integer programming model, the output is a sequence of operations that results in minimal cost (Figure 16b).

4.   Case study

In this section, the proposed method is verified by the correctness and efficiency of a virtual case study.

Figure 17a shows the point cloud scanned from the used part, which is required to be remanufactured. Figure 17b represents the CAD model of the final part, which has different functionalities compared to the used part. Especially, the final part has primitive features and 2.5D freeform feature.

The level set function for the final part is built from the CAD modelling on a design domain of size 150 × 150 × 150 with grid size $∆x$ (0.5 mm), as Figure 18 shown. In the figure, the modelling history is also given, and the mathematical formulation of each level set function is provided in the Appendix.

In terms of the used part, the RANSAC surface fitting technique is adapted and the max distance and angular distance variants are set as 0.002 m and 5 degrees, respectively. The surface fitting results are shown in Figure 19a. By implementing the k-NNG-based defect identification method, the defect points are differentiated from the point clouds Figure 19b. The level set function for the used part is constructed in Figure 20a with modelling history: ${{\mathrm{\Phi }}_{u} = \mathrm{m}\mathrm{i}\mathrm{n}(\mathrm{\Phi }}_{\mathrm{c}\mathrm{u}\mathrm{b}\mathrm{e}1}, {-\mathrm{\Phi }}_{\mathrm{c}\mathrm{u}\mathrm{b}\mathrm{e}2}, {-\mathrm{\Phi }}_{\mathrm{c}\mathrm{u}\mathrm{b}\mathrm{e}3}, {-\mathrm{\Phi }}_{\mathrm{c}\mathrm{u}\mathrm{b}\mathrm{e}4})$and formulations for each level set functions are given in the Appendix. The optimal pre-machining feature can be obtained from the defect points, and the used part is updated with the pre-machining feature, as shown in Figure 20b.

The proposed intersection part extraction algorithm is applied to find the relative position between the used part and the final part. The optimization results of the translation and rotation variables are listed in Table 3 and the output of the intersection part is shown in Figure 21.

For the next step, the intersection part is modified by considering the collision problem in the DED process. The DED nozzle is modelled in level set function as shown in Figure 22a and the vertical distance d is set as 10 mm. Figure 22b shows the deposition volume, which is obtained by subtracting the final part by the original intersection part. The motion space that represents all possible motions of the DED nozzle during deposition is derived from the proposed algorithm (Table 2); see Figure 22c. The result of the modified intersection part through subtracting the original intersection part by motion space is presented in Figure 22d.

By subtracting the final part and the used part by the modified intersection part, SFV and AFV can be collected, respectively. Then, each SF and AF are extracted, as presented in Figure 23. With respecting hybrid manufacturing rules, precedence constraints are applied to all features.

The parameters and machine resources that are used in the cost estimation are listed in Table 4. In this study, Ti-6Al-4V is used as the material for the DED process. The sequence optimization problem is solved by a branch-and-bound solver and the optimized process plans are presented in Table 5, and the optimal remanufacturing plan costs ＄ 1835.24.

5.   Conclusions

Nowadays, the combination of AM and SM in a single workstation is emerging to provide a more flexible and productive and capable manufacturing approach comparing with traditional manufacturing strategies. Because it utilizes the merits of AM and SM to add and remove features flexibly, HM has the potentials to raise remanufacturing technology to a higher level. In this study, taking to account the benefits of the HM process, hybrid DED-CNC manufacturing technology in a single workstation is investigated to remanufacture an end of life part (used part) to a new part (final part) with new functionalities, avoiding the material recycling process.

This paper has demonstrated a novel feature extraction algorithm and a cost-driven process planning method for hybrid DED-CNC manufacturing in a remanufacturing context. Specifically, starting from point clouds for the used part and a solid CAD model for the final part, geometry modeling is performed to transform the input to level set representations. Also, the defects on the used part are investigated and an optimal pre-machining feature is derived from the proposed method and applied to remove the perturbation caused by the defects. The feature extraction method, developed under the level set framework, is proposed as an automated process to extract the AM and SM features for remanufacturing process planning, which addresses the numerical calculation between two solid models. Moreover, the collision problems of DED nozzle during deposition are considered and this method provides collision-free motions. With the carefully developed hybrid DED-CNC cost model, the process planning work is converted to an -integer programming model as an optimization problem. Finally, the optimal process plan can be determined by solving the optimization problem.

For future work, there are more research works need to be done to expand the current methodology. A practical case study is required to validate the work presented in this paper. In the next stage, the relevant use case can be remanufacturing for high-value component, such as mold, die, and turbine blade. Moreover, in the present study, the AM and SM processes are considered as a 3-axis type. However, the feature extraction and process planning problems for multi-axis capability need to be addressed. In addition, the fixture design, AM support structures ^[53] in the real cases need to be investigated in the future work. In the metal additive manufacturing, the residual deformation is a significant problem and it also may affect the feature extraction and process planning, so it can also be explored in the further studies.

Acknowledgments

We express our appreciation to the other team members in the Laboratory of Intelligent Manufacturing, Design and Automation (LIMDA) group for sharing their wisdom during the research. The authors acknowledge the NSERC (Grant No. RGPIN-2017-04516 and CRDPJ 537378-18) and Minister of Economic Development, Trade and Tourism (through Major Innovation Project) for funding this project.

Conflict of Interests

All authors declare no conflicts of interests in this paper.