Research on the Method of Solving Kinematics Parameters of Three-axis Dynamic Centrifuge


 With the development of aerospace technology, the maneuverability of various types of aircraft continues to improve, and these aircraft experience an overload environment with high acceleration and high acceleration rate. Due to the special influence brought by the high acceleration rate, the dynamic centrifugal test technology, which is different from the traditional steady-state centrifugal test technology, came into being. The steady-state centrifuge has only one rotor, while the dynamic centrifuge has multiple rotors; so, the relationship between the acceleration of the control point and the motion parameters of the rotors is more complicated. Therefore, a key issue of the dynamic centrifugal test technology is the inverse kinematics of the dynamic centrifuge, which is to calculate the kinematic parameters of the dynamic centrifuge according to the expected acceleration environment that needs to be simulated on the centrifuge. After the kinematic parameters is calculated, the control target of each rotor of the dynamic centrifuge could be known, then the expected acceleration environment could be produced. In this paper, 1) on the basis of the predecessors, the equation for solving the angular velocity of the main arm of the centrifuge is improved; 2) and then a time step adaptive method is proposed, which takes into account the calculation accuracy and efficiency. As a result, an inverse kinematics algorithm that is more accurate and adaptable to various acceleration history curve is obtained. Finally, the inverse kinematics algorithm in this paper is verified through experiments and numerical simulations.


INTRODUCTION
With the development of aerospace technology, the maneuverability of many aircraft continues to improve [1]. These aircraft experience an overload environment with high acceleration and high acceleration rate [2][3][4][5][6]. This high acceleration rate will affect the performance of the electronic components on the aircraft, and change the dynamic characteristics of the lightweight and thin plate shell structure on the aircraft [7][8][9][10]. At the same time, it also places higher requirements on the pilot's anti-overload capability [11][12][13]. Therefore, the dynamic centrifugal test technology, which is different from the traditional steady-state centrifugal test technology [14][15][16][17][18][19], came into being.
The steady-state centrifuge has only one rotor. In the steady-state centrifugal test, it is stipulated that the loading and unloading rate of the centrifuge must be slow enough. Therefore, the tangential acceleration is ignored in the steady-state centrifugal test, and the test acceleration a is considered to be Rω1 2 . According to the test acceleration a and the centrifuge arm length R, the angular velocity of the centrifuge can be quickly obtained as / aR   . For the dynamic centrifugal test [20], 1) in order to generate an overload environment with a high rate of change, the loading rate of the dynamic centrifuge may be very large, and 2) the dynamic centrifuge has multiple rotors (e.g., the GLJ-3R dynamic centrifuge shown in Fig. 1 has three rotors: main arm, outer frame and inner frame). Therefore, on the one hand, the axial and tangential acceleration must be both considered, and on the other hand, the coordinate transformation method of multi-body dynamics needs to be applied to obtain the acceleration at the control points on the dynamic centrifuge. As a result, the forward kinematics of the dynamic centrifuge is far more complicated than that of the steady-state centrifuge, and ultimately leads to the inverse kinematics of the dynamic centrifuge (that is, calculating the kinematic parameters of the dynamic centrifuge according to the size of the centrifuge and the expected acceleration environment that needs to be simulated on the centrifuge) is extremely complicated. Only by solving the inverse kinematics problem of the dynamic centrifuge can the control target of each rotor of the dynamic centrifuge be obtained, and then the expected acceleration environment could be generated. Therefore, the inverse kinematics problem of the dynamic centrifuge has received many attentions. In [21], a new type of three-channel centrifuge is proposed, and the analytical solution of the kinematic parameters of centrifuge the under the condition of uniform rotation is obtained. In [22], the biaxial acceleration based on the simplified kinematics model is obtained. In [23], the angular velocity of the main arm during the loading phase is solved. However, the solution in the unloading phase is more difficult. To solve this problem, in [24] Jacobi elliptic function is introduced to solve it, but the error of this method is a little large, reaching 0.49g (g is the acceleration of gravity). In [25], the two-dimensional interpolation is used to estimate the motion parameters of the unloading section, but it is necessary to establish a database before using this method. In [2], the method of constructing symmetrical loads is presented to solve this problem, the calculation accuracy is high, and the process is simple. However, the size of the time step has a great influence on the accuracy of the method, especially during the acceleration platform stage (that is, the stage where the acceleration is constant). In this paper, the equation for solving the angular velocity of the main arm of the centrifuge is derived again, and then a method which automatically selects the size of time step according to the accuracy is proposed for solving the equation, so both the calculation accuracy and efficiency are took into account in this paper. Finally, the algorithm in this paper is verified through experiments and numerical simulations, and the results shows that our inverse kinematics algorithm is accurate and adaptive to various kinds of acceleration history curve.
The contents of this paper are organized as follows. The forward and especially inverse kinematics of GLJ-3R dynamic centrifuge is introduced in Section 2; then the experiment and numerical simulations are shown in Section 3 to verify and determine the accuracy of our inverse kinematics algorithm; at last, the conclusion is in Section 4.

Forward Kinematics of GLJ-3R Dynamic Centrifuge
The GLJ-3R dynamic centrifuge in Fig. 1 contains three rotors: main arm, outer frame and inner frame, and its schematic diagram can be simplified as shown in Fig. 2. The three rotors rotate around the Axes 1, 2, and 3 respectively. The rotation angle of the main arm relative to the ground is q1, the rotation angle of the outer frame relative to the main arm is q2, and the rotation angle of the inner frame relative to the outer frame is q3; the coordinate system O0-x0y0z0 is an inertial system, O-x1y1z1 is the body fixed frame of the main arm, O-x2y2z2 is the body fixed frame of the outer frame, and O-x3y3z3 is the body fixed frame of the inner frame. This article stipulates that when q1 is zero, x1 and x0 coincide; when q2 and q3 are both zero, O-x1y1z1, O-x2y2z2 and O-x3y3z3 coincide.   (1) where (1) a represents the coordinates of a in O-x1y1z1, R represents the length of O0O, the superscript "T" represents matrix transpose. The coordinates of a in O-x3y3z3 can be obtained through coordinate transformation as In the dynamic centrifugal test, the test object is fixed at point O on the inner frame, so the acceleration of the object in its body-fixed frame O-x3y3z3 is exactly (3) a . Therefore, the control objective of the dynamic centrifugal test is to ensure that (3) a is equal to the expected acceleration of the test object. From Eq. (2), it can be seen that (3) a is completely determined by the parameters R, 1 q &, 1 q && , 2 q , and 3 q ; here R is a constant. The goal of this paper is to calculate the kinematic parameters 1 q &, 1 q && , 2 q , and 3 q of the dynamic centrifuge according to the expected acceleration (3) a of the object.

Inverse Kinematics of GLJ-3R Dynamic Centrifuge
Calculating the kinematic parameters 1 q &, 1 q && , 2 q , and 3 q of the dynamic centrifuge according to the expected acceleration (3) a of the object is an inverse kinematics problem. There are some nonlinear terms (e.g., 2 1 q & , 2 sin q , 3 sin q ) and the first two-order derivatives of q1 in Eq. (2), so it is rather difficult to solve this equation directly. In order to simplify the solving, the magnitude a of the expected acceleration a is firstly analyzed.
It can be seen that a is completely determined by the kinematic parameter q1 of the main arm. According to Eq.
(3), q1 can be calculated firstly, so as to avoid the direct solving of Eq. (2). Eq. (3) is a high-order nonlinear differential equation about q1, which is not easy to solve directly. The following solution is based on the idea of the literature [2]. A small time period t  at the current moment t is firstly took. Assuming that 1 q && changes linearly during this period, then the 1 q & at the next moment can be calculated as is calculated by the rectangle integration formula, while the trapezoidal integration formula is used in this paper to refine the accuracy. The following equation can be obtained by substituting Eq. (4) into Eq. (3). 4 3 It can be seen that Eq. (5)  t  respectively. While in this paper, these two terms are both considered to refine the accuracy of Eq. (5). Besides, it is not too difficult to solve Eq. (5) even if these two terms are both considered, since this equation can be solved by the general solution formula of unary quartic equation. It is worth noting that in the unloading phase (that is, the phase when the acceleration magnitude a drops), there may be a situation where Eq. (5) has no real roots. To solve this problem, the method of constructing symmetrical loads is presented in [2] to solve the problem, the calculation accuracy is high, and the process is simple. However, the size of the time step has a great influence on the accuracy of the method, especially during the acceleration platform stage (that is, the stage where the acceleration is constant). In this paper, a method which automatically selects the size of time step according to the accuracy is proposed for solving the Eq. (5), so both the calculation accuracy and efficiency are took into account in this paper.

Method of adaptively selecting the size of time step
It can be known from Eq. (3) that The acceleration history curve in Fig. 3 Then it can be obtained that the solution will be as follows: 1 1 oscillates at a high frequency of 1 / 2 t  during the acceleration platform stage (this will be shown in Section 0). In fact, during the acceleration platform stage, the angular velocity of the main arm will be constant to ensure that the acceleration keeps unchanged; therefore, during the acceleration platform stage, 1 q && will tend to 0, 1 q & and will tend to / aR . Therefore, this high-frequency oscillation solution is a defect of the above algorithm. When using the algorithm above to calculate the control parameters of the dynamic centrifuge, the algorithm should be improved. The improved program flow chart is shown in Fig. 4. In the case of acceleration history as shown in Fig. 3, ω1 has not reached its stable value during the acceleration platform stage at time t1, and there is still room for its rise, so α1 > 0 is a true solution, here ω1 = 1 q & and α1 = 1 q && . The tangential and normal accelerations of point O on the main arm are Rα1 and -Rω1 2 respectively, so abs(α1/ω1 2 ) represents the relative magnitude of the tangential acceleration to the normal acceleration, and the abs(α1/ω1 2 ) ≤ tol in the program means that within a certain accuracy range, α1 = 0. At this moment, since angular acceleration is nearly zero, the angular velocity of the main arm has nearly reached the stable value during the acceleration platform stage. Therefore, it is stipulated in this program that as long as abs(α1/ω1 2 ) ≤ tol or α1 > 0, the calculation will end, otherwise the time step dT will continue to be reduced.

Rotation angle of the inner and outer frame
After 1 q && at the next moment is solved according to Eq.
arctan , , Then the q3 at the next moment can be solved according to Eq. (2) as (3) (3) where q3,0 is the q3 at the current moment, , so the error of solving q2 and q3 caused by the magnitude of a will be eliminated.

SIMULATION
In this section, the correctness of the inverse kinematics algorithm in this paper is verified through an experiment firstly, and then the accuracy of this algorithm is analyzed through several numerical simulations.

Simulation 1: Expected Acceleration with Two Non-Zero Components
The expected acceleration history curve shown in Fig. 5 is considered here, and the inverse kinematics algorithm in this paper is used to calculate kinematic parameters of the GLJ-3R dynamic centrifuge. Then the calculated parameters are used as the control target of the GLJ-3R dynamic centrifuge, so that the expected acceleration environment could be simulated on the centrifuge. The correctness of the algorithm can be verified by comparing the expected acceleration history curve with the acceleration history curve generated by the experiment.

Fig. 5 Expected acceleration history curve
In this simulation, the initial value of Δt is chosen as 0.01s. The kinematic parameters of the centrifuge obtained by the inverse kinematics algorithm of this paper are shown in Fig. 6, where the solid line represents the measured data of the kinematic parameters of the centrifuge, and the dotted line represents the calculation result, which is used as the control target of the centrifuge. It can be seen that the measured data (control result) basically coincides with the target curve (calculation), the angular velocity error is within 0.6 r/min, and the angle control error is within 0.5 deg. Since there is no measured data of the angular acceleration of the Axis 1, the calculation results of α1 are not given here.  Fig. 7. The solid line in the figure represents the expected acceleration, and the dotted line represents the measured acceleration. It can be seen that if the inverse kinematics solution in this paper is used as the control target of the dynamic centrifuge, the expected acceleration history curve will be simulated on the centrifuge, and the maximum error is about 0.2 g (g is the acceleration of gravity, taking 9.81m/s 2 ). The above-mentioned error consists of two parts: 1) the acceleration error of the inverse kinematics solution, this error causes the deviation between the control target and its expected value since the control target is calculated by the inverse kinematics algorithm, 2) the control error, that is, the error caused by the deviation between the control result and the control target. For this article, the main concern is the former error. The error is shown in Fig. 8. It can be seen that the acceleration error of the inverse kinematics solution is less than 0.0004 g, which is much less than the total error of 0.2 g. Therefore, the acceleration error in the experiment mainly comes from the control error.

Simulation 2: Expected Trapezoidal Acceleration
The expected trapezoidal acceleration history curve shown in Fig. 9 is considered in this simulation, and its math expression is where g = 9.81m/s 2 . In order to use the inverse kinematics algorithm in this paper, firstly, the acceleration curve needs to be divided into non-decreasing segment and non-increasing segment, and then these segments should be solved in forward and reverse directions respectively. The dividing point is the * in Fig. 9. In this simulation, the initial value of Δt is chosen as 0.01s. The kinematic parameters of the centrifuge obtained by the inverse kinematics algorithm is shown in Fig. 10. It can be seen that the angular velocity ω1 of the Axis 1 has a sudden change near t = 2s. In order to analyze the cause of this sudden change, the expected acceleration curve during the timespan of [0s, 2s] is observed firstly. It can be seen that before t = 1s, a increases monotonously and then remains at the level of 1g. Then the angular velocity curve should be checked, it can be seen that ω1 has not reached a stable value at t = 1s. It is easy to know that ω1 will continue to increase when a is rising; so when a just reaches the platform stage, ω1 must be less than its stable value due to the existence of the angular acceleration. Therefore, when a just reaches the platform stage, ω1 will continue to increase, and the angular acceleration α1 will gradually reduce to zero from a positive value, and when α1 reduce to zero (at about t = 1.29s in this simulation) ω1 just reaches its stable value. The result in Fig. 10 also confirms this. Therefore, the reason for this sudden change in the curve of ω1 could be analyzed as follows.
1) When solving the kinematic parameters corresponding to the last segment (2s < t ≤ 3s) of the acceleration curve in the reverse direction, since the dividing point is at t = 2s, the reverse calculation also ends at t = 2s, and ω1 has not yet reached the stable value at t = 2 + s, here 2 + represents a number that is a little more than 2.
2) However, for the acceleration segment of 0s ≤ t ≤ 2s, the calculation direction is forward, and ω1 remains at the stable value in the time interval of [1.29s, 2s], namely ω1 reaches the stable value when t = 2s. Therefore, the angular velocity curve is indeed discontinuous at t = 2s.
The history curve of angular acceleration α1 can be derived from the angular velocity curve, and then the history curve of acceleration a at point O on the dynamic centrifuge can be calculated by bringing the obtained ω1, α1, q2 and q3 into Eq. (2). The acceleration error of the inverse kinematics solution in this paper can be obtained by calculating the difference between the expected acceleration history curve and the calculated acceleration history curve. The error is shown in Fig. 11. Since the x and y components of the acceleration error are much smaller than the z component, only the z component result is given here. It can be seen from the figure that the inverse kinematics solution has a large error at t = 2s, which is caused by the sudden change of angular velocity. Sudden changes in angular velocity will cause the angular acceleration at t = 2s to be very large, thereby introducing an obvious tangential acceleration component, which will eventually cause this large error at t = 2s. Fig. 11 Acceleration error of the inverse kinematics solution After the analysis (in the penultimate paragraph above) of the cause of the sudden change of the obtained angular velocity, it can be seen that if the second dividing point of acceleration curve is changed to the midpoint of the overload platform stage, the problem of the sudden change of the obtained angular velocity could be solved. The dividing point after this change is shown in Fig. 12.

Simulation 3: Expected Acceleration from [2]
The following expected trapezoidal acceleration history curve shown in Fig. 15 is considered in this simulation, and its math expression is expressed in Eq. (11). As shown in Fig. 15, the dividing point are also represented by the mark "*", and here the dividing points at the end of the acceleration platform stage are all avoided.  Fig. 15 Expected trapezoidal acceleration history curve in [2] In this simulation, the initial value of Δt is chosen as 0.01s. The kinematic parameters of the centrifuge obtained by the inverse kinematics algorithm is shown in Fig. 16, and the corresponding acceleration error of the inverse kinematics solution is shown in Fig. 17. Similarly, since the x and y components of the acceleration error are much smaller than the z component, only the z component result is given here. As shown in Fig. 17, if the initial value of Δt is chosen as 0.005s, the corresponding acceleration error could be further reduced to 0.01g, while in [2], the error is 0.06g. The method of adaptively selecting the size of time step has been given in Section 2.2.1. If this method is not available, namely the time step is fixed; the angular acceleration curve corresponding to this example is shown in Fig. 18. It can be seen from the figure that there is high frequency oscillation in the angular acceleration curve. The theory in Section 2.2.1 has predicted this high-frequency oscillation, and the simulation results here verify the theory.

Simulation 4: Expected Acceleration with Three Non-Zero Components
The following history curve of acceleration with three non-zero components shown in Fig. 19 is considered in this simulation. Its math expression is shown in Eq. (12). 10g sin( ) 5gcos( ) 7.5g cos (2 ) x y z at at at         Fig. 19 History curve of expected acceleration with three non-zero components In this simulation, the initial value of Δt is chosen as 0.01s. The kinematic parameters of the centrifuge obtained by the inverse kinematics algorithm is shown in Fig. 20, and the corresponding acceleration error of the inverse kinematics solution is shown in Fig. 21. By comparing Figs. 8, 14, 17 and 21, it can be seen that, the error of this algorithm is relatively large (reach 0.02g) at the corner of the acceleration curve. If the acceleration curve is smooth and without corner, the magnitude of the error will be rather small (about 10 -4 g).

CONCLUSIONS
Firstly, the correctness of the inverse kinematics algorithm in this paper is verified through an experiment, and then a series of numerical simulations show that: 1) it is more suitable to choose the dividing point of the acceleration curve at the midpoint of the overload platform stage, instead of at the corner of the acceleration curve; 2) the method of adaptively selecting the size of time step in the Section 2.2.1 is effective and avoids the high-frequency oscillation type error of the kinematic parameter solution; 3) the algorithm has a large error at the corner of the acceleration curve, which can be up to 0.02g, while if the acceleration curve is smooth and has no corner, the magnitude of the error is about 10 -4 g; 4) this algorithm has strong applicability and is suitable for expected acceleration curves with one, two or three non-zero components, or acceleration curves containing corners; 5) in this paper the algorithm in [2] is improved, and the calculation accuracy is also refined.