A mathematical model of the fractional differential method for structural design dynamics simulation of lower limb force movement step structure based on Sanda movement

In the above discussion, we have only considered the improvement of prediction methods, thinking from a methodological point of view, we should also consider the improvement of correction, and we can hope that the appropriate improvement will bring the return of increased accuracy. The following recommended solutions are commonly used [7]: (6) yk+10=yk−1m+2hf(tk,ykm) y_{k + 1}^0 = y_{k - 1}^m + 2hf\left( {{t_k},y_k^m} \right) (7) yk+1m=ykm+h2[ f(tk,ykm)+f(tk+1,yk+1i−1) ]i=1,2,⋯,m \matrix{ {y_{k + 1}^m = y_k^m + {h \over 2}\left[ {f\left( {{t_k},y_k^m} \right) + f\left( {{t_{k + 1}},y_{k + 1}^{i - 1}} \right)} \right]} \cr {i = 1,2, \cdots ,m} \cr } where yk+1m y_{k + 1}^m is the final calculation result at time t_k+1, and the number of corrections m is determined according to the following formula (8) |εa|=|yk+1i−yk+1i−1yk+1i|×100%≤ε \left| {{\varepsilon _a}} \right| = \left| {{{y_{k + 1}^i - y_{k + 1}^{i - 1}} \over {y_{k + 1}^i}}} \right| \times 100\% \le \varepsilon where is the pre-specified error accuracy.

Estimated with errors (8)–(11) (12) E=yk+10−yk+1m5=−112h3f‴(tk,ξ) E = {{y_{k + 1}^0 - y_{k + 1}^m} \over 5} = - {1 \over {12}}{h^3}f'''\left( {{t_k},\xi } \right) If |E| is greater than the previously specified error accuracy, then decrease the step size.

In integral formula, (13) ∫yk−nyk+1y′dy=y(tk+1)−y(tk−n)=∫tk−ntk+1f(t,y(t))dt⇒yk+1=yk−n+∫tk−ntk+1f(t,y(t))dt \matrix{ {\int_{{y_{k - n}}}^{{y_{k + 1}}} y'dy = y\left( {{t_{k + 1}}} \right) - y\left( {{t_{k - n}}} \right) = \int_{{t_{k - n}}}^{{t_{k + 1}}} f\left( {t,y\left( t \right)} \right)dt} \hfill \cr { \Rightarrow {y_{k + 1}} = {y_{k - n}} + \int_{{t_{k - n}}}^{{t_{k + 1}}} f\left( {t,y\left( t \right)} \right)dt} \hfill \cr } By replacing the integrand f(t, y(t)) with the n-degree interpolation polynomial P_n (t) of n + 1 nodes (t_i, f_i (t_i, y_i)) , i = k − n, k − n − 1, ··· , k + 1 the general linear multi-step method is obtained. As when n = 1, (14) yk+1=yk−1+2hf(tk,yk) {y_{k + 1}} = {y_{k - 1}} + 2hf\left( {{t_k},{y_k}} \right) As when n = 2.

In particular, by replacing the integrand f(t, y(t)) with the cubic interpolation polynomial P₃(t) at four nodes (t_i, f_i (t_i, y_i)), i = k − 3, k − 2, k − 1, k we get the following: (15) yk+1=yk+h24(55fk−59fk−1+37fk−2−9fk−3) {y_{k + 1}} = {y_k} + {h \over {24}}\left( {55{f_k} - 59{f_{k - 1}} + 37{f_{k - 2}} - 9{f_{k - 3}}} \right)

The four-node (t_i, f_i (t_i, y_i)) , i = k − 2, k − 1, k, k + 1 and cubic interpolation polynomial P₃ (t) are substituted and we get the following [8]: (16) yk+1=yk+h24(9fk+1+19fk−5fk−1−fk−2) {y_{k + 1}} = {y_k} + {h \over {24}}(9{f_{k + 1}} + 19{f_k} - 5{f_{k - 1}} - {f_{k - 2}}) The errors are 251720h5y(5)(ξ) {{251} \over {720}}{h^5}{y^{(5)}}(\xi ) and −19720h5y(5)(ξ) - {{19} \over {720}}{h^5}{y^{\left( 5 \right)}}\left( \xi \right) respectively. The accuracy of formula (16) is higher than that of formula (15).

Adams Forecast-Correction Formula: (17) {y¯n+1=yn+h24(55fn−59fn−1+37fn−2−9fn−3)yn+1=yn+h24(9f¯n+1+19fn−5fn−1+fn−2) \left\{ {\matrix{ {{{\bar y}_{n + 1}} = {y_n} + {h \over {24}}\left( {55{f_n} - 59{f_{n - 1}} + 37{f_{n - 2}} - 9{f_{n - 3}}} \right)} \hfill \cr {{y_{n + 1}} = {y_n} + {h \over {24}}\left( {9{{\bar f}_{n + 1}} + 19{f_n} - 5{f_{n - 1}} + {f_{n - 2}}} \right)}\hfill \cr } } \right.

Among them, f_n (t_n, y_n), f¯n+1(tn+1,y¯n+1) {\bar f_{n + 1}}\left( {{t_{n + 1}},{{\bar y}_{n + 1}}} \right) .

Note. Formula (18) cannot provide the three initial values required by itself and requires the cooperation of Runge-Kutta method.

Pro/E was used to model the mechanical structure of the new Sanda lower limb powering robot. Figure 3 shows the parts model of the training robot, and Figure 4 shows the assembly model in Pro/E.

Fig. 3

Fig. 4

The Pro/E model is imported into Adams, and the Adams model is generated through interface software of mechpro2005 of Adams and Pro/E [10], and the simulation analysis is performed in Adams. Model transformation includes three steps: defining rigid body, defining constraints and model transformation. Add necessary constraints and forces to the system in the Adams environment and improve the imported Adams model. The necessary forces and drives are added to define the reference quantities of the material density, stiffness, elastic deformation and acceleration of gravity of each rigid body. Finally, the interaction forces between various rigid bodies are defined, including the contact force between the hub and the ground, the acceleration of gravity and the motion drive. The built Adams dynamic model is shown in Figure 5.

Fig. 5

The specific statistics of each component, rigid body and constraint involved in the modelling process are shown in Table 1.

Given an omnidirectional wheel angular velocity of π6 rad/s, the curve of omnidirectional wheel angle change in the rotating state is shown in Figures 6, and 7 shows the robot trajectory curve and rotation angular velocity curve. Figure 6 samples the angular velocity of each wheel. Since the simulation purpose is to analyse the relationship between the angular velocity of each wheel and the robot trajectory, this sampling is used as a known input. Figure 7 reflects the change of the centre of gravity of the robot in space during the rotation. On the horizontal plane, the centre of gravity of the robot exhibits a sinusoidal change on the x and z axes, which meets the system control requirements. There is a 1 mm disturbance on the y-axis change, indicating that there is a vibration phenomenon during the rotation of the robot waiting for buffering [11]. Judging from the maximum partial velocity in the three directions of X, Y and Z, each joint has the largest X direction, the Z direction is second and the Y direction is the smallest. For example, for the toe, the maximum speed in the X direction is 13.02 m/s, the Z direction is 11.55 m/s and the Y direction is only 7.09 m/s. Another example is the knee joint, which is also the largest in the X direction, with a maximum average speed of 6.43 m/s, followed by the Z direction with 5.09 m/s and Y direction with only 4.19 m/s.

Fig. 6

Fig. 7

Given the angular velocity of the omnidirectional wheel as π6 rad/s, the curve of the omnidirectional wheel angle in the forward running state, that is, the input curve, is shown in Figure 8. The trajectory and speed curve of the robot in the forward state, that is, the change curve of the centre of gravity of the robot corresponding to Figure 7, are obtained under a rotating state. It can be seen from the figure that among the velocity components, only the Y direction (lateral or striking direction) has the highest velocity value now on striking, and the maximum values of other components are before the striking point. Obviously, the Y direction is the effective strike direction. Now of the strike, it is necessary to strike at the maximum speed value in order to produce a crisp striking sound or to achieve the purpose of damaging the opponent, and the other directions are ineffective speed directions. From the perspective of kinematics, the magnitude of other velocity components is only meaningful for shortening the swing time but not significant for the effectiveness of the strike. Therefore, in order to improve the effectiveness of the strike, the Y direction (side leg direction) and the maximum speed should be in the direction of the athlete's efforts.

Fig. 8

Figure 9 reflects the disturbance in the x-axis direction during the starting process. The maximum peak-to-peak value of the disturbance is 4 mm, and it becomes gentle after 15 s. According to the analysis, the cause of the disturbance is that the centre of gravity of the robot and its geometric centre are not in the same vertical direction, and the forces on the wheels are uneven when starting; it reflects the 1 mm disturbance on the y-axis change; the z-axis travels at a constant speed. Since the disturbance peaks sampled in Figure 9 are relatively small for the entire system, they have little effect on the overall system effect. Now of impact, the three velocity components of the hip and knee joints are relatively small, all <2 m/s. Among the X, Y and Z speed components of the ankle joint point, the speed in the Y direction is the smallest when it hits 4.03 m/s, while the speed component values in the X and Z directions have no significant difference (P> 0.05). The average is about 5.4 m/s. At the toes, there was no significant difference between the three speeds of X, Y and Z (P> 0.05), as now of on impact, with an average of about 7 m/s.

Fig. 9

According to the dynamic simulation results obtained in Adams, it can be analysed that the centre of gravity of the robot is not at its geometric centre, which will certainly cause some interference to the control of the robot. Improving the design structure to make it more symmetrical can reduce this tendency of interference. From the perspective of the dynamic mathematical model analysis, Adams simulation results are consistent with the dynamic mathematical model. When the given input is the same, the dynamic simulation of the design is consistent with the mathematical model results, indicating that the design structure meets the training requirements.