Trajectory prediction of three-dimensional forming tube based on Kalman filter

This article focuses on the prediction of forming trajectory and process optimization during the forming process for the variable curvature tubes. Firstly, through cubic B-spline interpolation, the geometric characteristics of the axis of the target tube are obtained. The whole tube is “separated and then integrated,” and the relationship between geometric parameters and processing parameters is established to obtain the initial process parameters. Based on the extended Kalman filter (EKF) algorithm, the motion model and observation model of tube forming are presented section by section, and the relevant calculation and analysis are carried out. The forming trajectory has been predicted; moreover, the processing parameters are optimized during the processing process, in which the effectiveness of the processing optimum scheme is illustrated.


Introduction
Metallic tubular parts are widely used in many high-end industries such as aviation, aerospace, and automobile. In modern turbofan aero-engines, it is reported that the weight of an external piping system is nearly 60 ~ 70% of the total weight of external structures [1][2][3]. According to statistics, 50% of all aero-engine failures are caused by the external piping system, electric wire, and sensors [4,5]. Therefore, improving the usability of the external piping system is an important strategy for extending the life span of an aeroengine. When designing the structures of the external piping system, the influences of general layouts, rigid strength, flow resistance characters, and vibration performance should be considered.
A good tubular product should possess some special characteristics, such as lightweight, high strength, low impedance, and high vibration resistance. In geometrical characters, the tube axis should be smooth and uniform in curvature distribution. However, under current conditions, the tubular parts in aero-engines are mainly formed by the rotary-draw-bending (RDB) process, which can form 2D shapes precisely but can not be applied for manufacturing complex spatial tubular products. As a result, the maximum usability of tubular parts is limited in the RDB process.
To overcome the deficiencies in the RDB process, free bending technology is proposed and considered an innovative process for manufacturing freeform bending tubes and profiles. The new bending technology inherits the high flexible characters in the traditional push bending process. Without changing tools, different bending curvatures and bending angles can be obtained easily [6][7][8]. Besides, the free bending technology extends the forming scope of the push bending process greatly. A typical model for the free bending process is shown in Fig. 1. The mobile die can shift around the fixed die freely. The clamp can provide a pushing force and a torque on the tube axis. When forming process begins, the mobile die deflects from the original position and forces the tube blank to be plastic. At the same time, the clamp forces the tube blank to move along the feeding direction. If torque is added to the clamp, a spatial tubular product can be formed easily.
In the current literature, Guo et al. [9,10] established several mathematical models for determining the bending die positions during free bending and investigated the strain stress distributions along bending lines numerically. Plettke et al. [11][12][13] investigated the relationship between mobile tool position and forming curvatures for push bending numerically and established the curvature distribution along bending lines when manufacturing arcshaped products. Vatter and Plettke [14] investigated the influence of mobile die position, tools stiffness, and clamp motion on bending curvature distributions and finally built a series of forming theories under single variable conditions. Engel et al. [15] investigated the influence of mobile tool loci on forming curvatures distribution and designed a series of strategies for reducing the length of transition zones on bending lines. Based on the existing experimental datum, Groth et al. [16] summarized the forming characters of curvature distribution along bending lines. When the bending process begins, the bending curvature would experience a viscous damping vibration before reaching constant. At the end of the bending process, there would be a sudden decrease in curvature values. Based on such a phenomenon, the curvature distribution along the bending line is divided into three parts: transition zone, constant curvature zone, and outfeed zone in Fig. 2.
It should be pointed out that the mentioned references above are mainly focused on arc-shaped bending instances. For 3D complex tube bending, the existing researches still need to be improved. Comparing with the arc-shaped bending, 3D bending is more complex. During the forming process, there exists a continuous change of bending planes along the bending line. Existing research has also shown that the torque can also influence the bending curvatures [17,18]. To precisely form a spatial tube product, an extensive relationship among mobile die motion, clamp motion, and geometrical characters of forming products should be proposed.
However, due to the irregular shape of these variablecurvature tubes, their geometric and physical non-linearity is stronger than those ordinary constant-curvature tubes. Therefore, the theoretical analytical models' insufficiency and test measuring data should be short and on the other hand be not accurate. Moreover, the control of its forming accuracy is more difficult, especially for the precise forming of variable-curvature tubes.
In this paper, a piece-wise function is adopted to simplify the curvature distribution along the bending line. Then, a series of control parameters is extracted, and the relation between controlling parameters and mobile tool motion is investigated. Based on this investigation, a new strategy based on the extended Kalman filter (EKF) algorithm is suggested for the forming control.

Fundamental theories about spatial curves
According to the theories of spatial geometry, a typical spatial curve can be expressed as Fig. 3. Assuming that the natural parametric function of a spatial curve is r = r(s) and point O is on the spatial curve, the active coordinate system at point O can be marked by tangent vector t(s) , normal vector n(s) , and bi-normal vector b(s) in Fig. 3. The triangular pyramid constituted by t(s) , n(s), and b(s) is known as the Frenet mark. The Frenet mark can reflect the variation tendency of curvature and torsion on the spatial curve.
According to the theories of spatial geometry, the relationship of t(s) , n(s) , and b(s) can be expressed as Eq. (1).
where k(s) and (s) refer to the curvature and torsion separately. k(s) and (s) can be calculated by Eq. (2).

Principle fundamental
The main features of the EKF are the Taylor series expansion for the non-linear system, and if just one order term is considered, this non-linear system would be treated as linear. After this approximation, the state evaluation should be done through the EKE algorithm framework. Commonly, a nonlinear system can be described as Eqs. (3)  where Eq. (3) includes the control function u k−1 and noise w k−1 ; function f will establish the mapping between the current state at time t and the state at time t − 1. In Eq. (4), function h represents the non-linear relation between the measurement variable z k and state variable x k .

Evaluation part
State evaluating: Co-variables predicting: 2. adjustment part Matrix of transmission gaining: Filtering equation: Up-dated co-variables: The above five equations from Eqs. (5) to (9) are the key elements for EKF, and the Jacobian matrix through the differential of function f with x is as And in the same way, the Jacobian matrix through the differential of function h with x is shown as: After the linear processing using the one order Taylor series, we have the following formulas.
State formula: Observation formula: The key characteristic of the EKF algorithm consists in the Kalman gaining K k , and the Jacobin matrix H can transfer useful observing information. Put simply, the framework of the EKF algorithm is presented in Fig. 4.
The application of EKF in 3D tube

The simulation of the free bending process
In this paper, a numerical model for push bending (shown in

Strategies for selecting the axis of the tube
When using the numerical method for tube springback analysis, a key procedure is to extract the tube  axis from the FE model. In this section, a strategy for extracting tube axis from the FE model considering the influence of wall thinning and cross-section distortion is introduced.
As shown in Fig. 7, nodes on the tube outer surface are defined as the intersections of U-curves and V-curves (U-curves are defined as curves that are approximately parallel to the tube axis and V-curves are defined as curves that are approximately vertical to the tube axis). The nodes on the same V-curve are chosen and noted as x i , y i , z i . To determine the best cross-section plane, Eq. (14) is defined as an optimization function.
where n refers to the number of nodes; A , B , C, and D are parameters that need to be determined.
Equation (14) can be solved. When A , B , C , and D are determined, the best cross-section plane can be described as Eq. (15) in the Cartesian coordinate system.
If the nodes x i , y i , z i are projected to the cross-section plane above, Fig. 8 can be obtained. Node coordinates after projection are noted as Fig. 12, (x, y, z) is defined as the center of the ideal surface, of which the diameter is equal to that of the tube blank. Equation (16) is defined as a function to optimize the center coordinates. The solution of Eq. (16) is considered the center of the tube cross-section.
where R refers to the outer radius of the tube. Repeating the strategies above several times, a series of tube cross-sections and center points can be determined. If we connect the center points using a B-spline, the tube axis can be geometrically built as Fig. 9. After extracting the tube axis, the curvature distribution along the bending line can be analyzed in the CREO modeling system as Fig. 10.

The forecasting and optimization of formed tube trajectory
A new strategy for acquiring the forming parameters is proposed in this paper, which is based on segmented arcs (or segmented spirals in 3D) approximation. In this way, the motion model, which would be fundamental for the EKF solution together with the measurement model, will be more effective. And the forecasting and optimization of formed tube trajectory would be dealt with efficiently.

Initial processing parameters
In the case of the 2D tube in Fig. 11, eleven uniformly distributed points are selected up along its axis from a designed CAD model. Making interpolation through these points using the cubic B-spline, we can acquire the spline equation, and through differential operation, the curvature k(s) can be calculated, shown in Fig. 12. The fluctuation phenomenon in Fig. 12 came from the curvature calculating algorithm that is much more sensitive to the point number and its distribution. Therefore, the "lowness" method in MATLAB is used to wipe off this noise and smooth the curvature curve. Selecting parameters t(0 ≤ t ≤ 1) corresponding with eleven uniformly distributed points, the curvature k and associated radii are as in Table 1.
Then, making an approximation of tube axis through segmented arcs to guarantee the co-alignment of tangent and normal directions at the starting point on the second arc and the final point on the first arc, respectively. After finishing the approximation section-by-section in the same way mentioned above, the approximation's accuracy could be evaluated by Eq. (17), and a better approximation is taken in Table 2. where Δl i is the distance between the objective point and the correspondent arc endpoint; is the total averaged tolerance; is the maximum tolerance. After approximating, some transformation from the design model to the processing model should be done. For simplicity, one section arc is taken as an example. As shown in Fig. 13, if Generation of the tube axis using a series of center points Grant being free of jam, the bending die should not only rotate around fixed die indexed by 2 but also around the axis itself indexed as 3 − 2 . In this situation, if the radii R is available, the angle 2 , and 3 − 2 can be calculated from Eq. (18). When the arc length s, feeding speed v, and feeding time is given, we have the following simple relations where is the center angle of wanted tube arc.

Tube motion model
Correctly selecting state variables could make the system model simple, and the EKE algorithm is more efficient. In general, the system state variables evolve with time step-bystep. In this way, differential equations, which are formulated by one order Taylor series with time, are used to represent this evolution. Through differential regarding state variables, the  The schematic diagram is shown in Fig. 14. Firstly, when the forming process begins, the feeding speed varies from 0 to v, then the bending die rotates hypothetically to a position with the arc center angle 1 in a short time. Continually feeding, the arc no. 1 is formed, and 1 value has a change of v∕R from time k to time k + 1. After finishing arc no. 1 forming, the die moves to a new position with arc center angle 2 , and then repeating this step, all processing progress can be finished.
The motion equation for a certain arc is

Through one order Taylor series approximation, then we have
The state transformation equation is Making differential to state variables, we obtain the following Jacobian matrix, the state transformation matrix.

Observation model
Monitoring the formed tube instantly, shown in Fig. 15, we can get the observation model that we need. Because the EKE algorithm is adopted, the observation model should be non-linear. With an eye on point A at a certain time, the distance r from point B to point A and the angle Δ between axis x and line AB are selected as observing variables, which can be calculated from monitored coordinates x and y.
Observation variables The relationship between observed variables and state variables are The observation matrix H is As shown in Fig. 16, at first, monitoring the trajectory of point A from the time t = 0 , if t = t 1 ∕4 , while we think a quarter of the arc no. 1 was formed, we stop monitoring point A and continue our monitoring on point B. By the same way, if t = t 1 ∕2 , we stop the monitoring on point B. Through this process, the formed tube part ABC, which is just half part of arc no. 1 section, is available and this obtained half part will be regarded as an observation model. Based on a comparison between the design shape and this half part, further calculation relating to compensation could be considered on the second half tube part. As for the second half part monitoring is finished, just repeat the above steps for the next section.

EKE application
For calculating efficiently and easily, the diagonal noise co-variable matrix could be used in this study with smaller initial values P 0 = diag[0.001,0.001]. The matrix A and Q could be calculated by multiple experiments (Piippo A). Because the tube axis calculated through the FEM data does have some tolerances, the observing noise covariable matrix is assumed as R = diag[0.09,0.09] in this situation. According to reference (Piippo A.), we choose Q = diag[0.01,0.01,0.01,0]. Assuming the coordinates x i,c , y i,c are objective point i, and the point coordinates x i,e , y i,e on the trajectory are calculated by EKF, therefore, the deviation between these two points is The maximum and mean deviations are In this study, the maximum deviation limit is set to be max = 2 mm, and the mean value is ave = 1.5 mm. Running the EKF algorithm section by section, the 1st section arc is performed by processing parameters 1 = 0.1745 rad, v 1 = 30 mm/s, and t 1 = 2 s, simulating and monitoring the first arc part, then finishing the filtering operation both on the established motion model, and observing model taken from FEM simulation. As shown in Fig. 17, the true half tube part trajectory can be calculated by the EKF using motion and observation models, and since the deviation between the true shape and objective one is within the tolerance, so no adjusting operation is needed.
For the second half part, we choose the processing parameters as 2 = 0.3054 rad, v 2 = 30 mm/s, and t 2 = 2 s. After the EKF calculating and deviation analysis, shown in Fig. 18, the deviation is larger than the mean value ave and the projected true trajectory goes away from the objective one; therefore, processing parameters should be adjusted to make some compensation. Fig. 15 Schematic of monitoring the formed tube According to the deviation measured from the comparison between predicted trajectory and objective, we can determine a compensated direction. Based on careful analysis, the projected bending radii are too small and the bending angle is too larger; therefore, a new adjusted processing parameter bending angle ′ 2 = 0.2618 rad is used for the second half tube section, shown in Fig. 19.
And after this compensation, the mean deviation is well below the limit E ave , which means the true trajectory is redressed to a reasonable position closer to the designed shape.
A set of initial processing parameters 3 = 0.4538 rad, v 3 = 30 mm/s, and t 3 = 1 s is adopted for the third arc part of the tube; the deviation concerning the first half is within the limit value, shown in Fig. 19, so nothing would be done for the second half. After all, through the EKF calculating and compensating, a better result is shown in Fig. 20b, contrasting with a result without compensation in Fig. 20a. Figure 21 shows us the normalized processing parameters varying with time, the red dot line representing the adjusted die rotating values, and the black dot line the initial values.

Spiral based approximation for 3D tube
A spiral line can be drawn out by moving a point both along a circle with radii a and an axis vertical to the circle plane at some constant speeds. As shown in Fig. 22, radii a and  where the value a is the base radii of the column, b is pitch, and = t , where is the angle speed. If we choose the arc length s as the parameter, then we have Eq. (33). And the arc length can be calculated using Eq. (34).
The Frenet frame at starting point is The Frenet frame at finishing points is The normal vector always points to the curvature center. And curvature can reflect the bending intensity of the curve. The torsion of the curve can express the deviation tendency of this curve to the osculating plane. In other words, curvature and torsion are used to describe the bending and twisting intensity of the spatial curve.
When a spatial product is available, we firstly should define the representation and connection method of the tube elements so that the forming parameters ( , ) can be obtained. Generally, if the tube elements are small enough, they can be considered a series of spirals. Figure 23 shows the representation of the spatial tube by spirals. Two adjacent elements share the same tangent vector t(s) , normal vector n(s) , and binormal vector b(s) at the common joint point. For each spiral element, curvature and torsion can be determined. So, a series of parameters ( , ) can be obtained.

Geometric parameters
For a case of the 3D tube in Fig. 24, eleven uniformly distributed points are selected up along its axis. Making interpolation using these points based on the triple spline, we can get the spline equation, and through differential operation, we can calculate the curvature and torsion shown in Fig. 25. The curvature and torsion at each point are in Table 3.
According to this analysis, three sections of spiral lines with different features, shown in Table 4, were selected to construct an approximation curve that will be used to coordinate the processing parameters.

Processing parameters
For 3D tube forming, the rotating speed of the clamp should be considered another processing parameter. In the situation of the spiral, the base radii a and pitch b will be coordinated with processing parameters: die rotating angle, feeding speed, and clamp rotating speed.
In the simulation using ABAQUS, the die rotating angle, feeding speed, and clamp rotating speed will be set as boundary conditions, similar to 2D simulation. When the die rotating angle is = 0.14 rad, feeding speed v = 100 mm/s, and clamp rotating speed is set as 0.1 rad/s, 0.2 rad/s, 0.3 rad/s, and 0.4 rad/s respectively, the values a and b of formed spiral tubes are shown in Fig. 26.
It is clear that the base radii change a little and the pitch b increases much more with the rotating speed . And for a point at the spiral line, the curvature k and torsion could be calculated from Eq. (38).
The clamp rotating speed is After the above analyses, the relationship between geometrical parameters and initial processing parameters has been established.

Motion model
In the case of 3D tube forming, we choose status variables such as object position coordinates x, y, z, rotating angle , feeding speed v, and twisting speed .
The schematic diagram is shown in Fig. 27. Firstly, when the forming process starts, the feeding speed varies from 0 to v, then the bending die rotates to a position with the arc center angle 1 in a short time and simultaneously clamp twists at a speed of 1 . Continually feeding, the spiral no. 1 is formed, and 1 value has changed to 2 and (40) Through one order Taylor series approximation, then we have Making differential to state variables, we obtain the following Jacobian matrix that is the state transformation matrix.

Observation model
Monitoring the formed tube, shown in Fig. 28. For point A at a certain time, the distance r from point B to point A, the angle between axis z and line AB, the angle between  The maximum deviation is set as max = 3 mm, and the mean deviation E ave is ave = 2 mm. The 1st section arc is performed by processing parameters 1 = 0.1745 rad,  Fig. 29, the deviation between the true shape and objective one is within the tolerance, so no adjusting operation is needed. For the second part of the spiral, the initial processing parameters are set as 2 = 0.4494 rad/s, v 2 = 40 mm/s, and t 2 = 2 s.
Results from EKF, related to the first half of the spiral, are presented in Fig. 30a, together with objective and observation curves. The deviation between objective and EKF trajectory is becoming somewhat larger and larger along with the forming process. In this situation, an adjustment should be invoked. According to the base radii a and pitch b calculated from fitting the half EKF curve through spiral formula, it is clear that a set of smaller values of radii a and pitch b makes this kind of deviation. Through analysis, the bending angle ′ 2 = 0.6435 rad and twisting speed ′ 2 = 0.4639 rad/s are carried out after several errors and tries. Figure 30b shows a better solution for the second half part of this section. A set of initial processing parameters 3 = 0.3657 rad, 3 = 0.1838 rad/s, v 3 = 40 mm/s, and t 3 = 2 s is adopted for the third spiral part of the tube; the deviation concerning the first half is not within the limit value, shown in Fig. 31a; therefore, an adjustment has to be done for the second half.
On the analysis similar to the second section, the base radii a should be somewhat larger. Therefore, new processing parameters ′ 3 = 0.4049 rad and ′ 3 = 0.2320 rad/s are tried out, and a reasonable result is presented in Fig. 31b.
Finally, the deviation cloud chart, obtained from the adjusted as well as not adjusted processing parameters, is shown in Fig. 32 respectively. Apparently, it is inspiring for us to consider this proposed processing parameter adjusting strategy. Figure 33 shows us the normalized processing parameters varying with time, the red dot line representing the adjusted die rotating and twisting values, and the black dot line the initial values.

Discussions and conclusions
1. For complex tubes, according to the basic principles of free bending forming, the forming process of free bending is analyzed and the cubic B-spline curve is used to pick up the shape characteristic feature. This work lays a theoretical foundation for the subsequent correlation of geometric parameters and processing parameters. 2. The EKF algorithms were applied to analyze and calculate a typical nonlinear system with a concern of solving the problem of predicting the forming trajectory of the variable curvature tubes. 3. A finite element analysis model for the free tube bending forming is established, and the axis information of the formed tube is extracted. For typical plane curvature variable tubes, arc interpolation is used to obtain its curvature change, and its geometric parameters are obtained by segmented arc approximation; for typical spatial curvature variable tubes, its geometric parameters are obtained by segmented spiral approximation. A kinematic model and an observation model (numerical simulation) were established for the curvature variable tubes, and the EKF algorithm was used for the related calculation and process optimization, which verified the effectiveness of the algorithm in predicting the forming trajectory of the curvature variable tubes.