Abstract
In this study, the dynamic model for the herringbone planetary gear transmission system is established by the lumped parameter method based on the system dynamics and the Lagrange equation, and the impact of the support stiffness and the torsional stiffness on dynamic characteristics is studied. The research results have a guiding significance for the design of the herringbone gear transmission system. In this model, the herringbone gear is treated as a special gear coupled by 2 opposite helical gears, where the stagger angle, comprehensive meshing error, support stiffness, support damping, and load inertia are considered in the analysis of dynamics. Moreover, the dynamic characteristic of the carrier is considered as well. By calculating the meshing force curve of the transmission system, the impact of the stagger angle, supporting stiffness, and the torsional stiffness on meshing force and load sharing coefficient is analyzed. The results show that the stagger angle has an obvious impact on load sharing coefficient while it has little impact on maximum meshing force. And the support stiffness has a more obvious impact on the dynamic characteristics of the system. The recommendary support stiffness of the system is that all of the support stiffness of the sun gear, planetary gear, ring gear, and carrier is 107 N/m. The torsional stiffness has little impact on the dynamic characteristics of transmission system, except the torsional stiffness of planetary gear, and carrier has an obvious impact on load sharing coefficient. The commercial software ADAMS carried out dynamics analysis of the transmission system to verify the necessity validity of the theoretical analysis.
1. Introduction
Herringbone gear planetary transmission system, which is one of the most important parts of wind power generator, is wildly used in ship, aviation, and other transmission fields because of its distinctive advantages, such as small volume, light weight, large transmission ratio, strong bearing capacity, high transmission efficiency, and so on. The load distribution of planetary gears in the transmission directly affects the working life, stability, and reliability of gears. It is of great practical and engineering significance to carry out the dynamic load characteristics of herringbone gear transmission system.
Many scholars around the world have done many theoretical and experimental studies on the load sharing characteristics of planetary transmission. SHI Wan-kai et al. [1] studied the transmission dynamic characteristic of wind speed increase gearbox to study the effect on wind turbine gearbox dynamic load characteristics caused by stagger angle and left-right coupling torsional stiffness. ZHANG Lin-lin et al. [2] researched load sharing characteristic of herringbone gear transmission system influenced by meshing phase and conducted an experiment on different meshing phase planetary transmission system to study load sharing on single-stage herringbone planetary gear test stand. MO Shuai et al. [3–5] proposed a new dynamic model of herringbone star gear transmission system to study meshing error, tooth thickness error, base circle error, and assembly error influencing the load sharing characteristic of GTF aeroengine. Ming LI [6] established a model to predict the reliability of helicopter planetary gear train under the condition with partial load and analyzed structure and part load state of the transmission system in detail based on the study of load sharing characteristics of a helicopter gearbox. The impact of assembly error on load sharing characteristic of closed differential herringbone gear transmission system was researched by ZHU Zeng-bao [7]. LU Jun-hua [8] established a calculation model of 2K-H transmission system to analyze load sharing characteristic of planetary transmission system from the point of view of dynamics. XU Xiang-yang [9] studied impact of stiffness of flexible pin shaft on load sharing of star transmission system and verified the effectiveness of theoretical result. Kahraman A et al. [10–14] researched the impact of gear indexing error on dynamic response of a typical multigrid gear system and compared subharmonic resonance of response with numerical integration result directly to prove the accuracy of prediction and predicted other nonlinear phenomena by nonlinear model and studied the variation of the position errors of the planetary gear and the accompanying star wheel under the condition of light load in the small module planetary gear system as well. Howard I [15] proposed a new dynamic model of flexible supporting gear train to calculate the meshing stiffness with varied gear center distance and coupled it to dynamic model by iteration. Wang J et al. [16] developed a three free degrees torsional-vibration model of locomotive spur gear transmission system, where nonlinear gap, static transfer error, and time varying meshing stiffness and other factors were considered, and the influence of speed and stiffness of pinion on the dynamic characteristics of the transmission system is studied by numerical method. Fernandez-Del-Rincon A et al. [17] proposed a model which can simultaneously consider the internal excitation caused by variable meshing stiffness and the variable compliance results of bearing. The model combines with hybrid formula for calculating meshing force and nonlinear variable flexibility of bearing to prove a detail described for meshing force. Liang X et al. [18] focus on modeling, detection, and diagnosis of gearbox failure based on dynamics. The research result was limited on some key and basic aspects by literature review: gear meshing stiffness evaluation, gearbox damage modeling and fault diagnosis technology, gearbox transmission path modeling, and method verification. TAN Yuan-qiang [19] researched dynamic and load sharing characteristics of closed differential planetary gearbox by comparative trial, which proves that the average load performance of the closed stage is better than that of the differential stage, and the greater the load, the better the load performance of the two stages. The transmission performance of the planetary system with the three kinds of load sharing methods is better than the single load sharing system. Singh [20] established an analysis model of planetary transmission load sharing characteristic to physically explain basic mechanism of the phenomenon which causes inhomogeneous load. Wink et al. [21] researched the calculation method of load distribution of tooth surface and proved three different methods to solve load distribution of gear teeth in calculation of static transmission error of gear pair.
The herringbone gear planetary transmission system of wind power train was taken as the research object in this paper. The support stiffness, support damping, meshing stiffness, meshing damping, and other factors of herringbone gear in the condition of stagger angle which is 0 are overall considered to study their impact on dynamic characteristic. Compared with other papers, this paper especially researched the impact of carrier on dynamic characteristic and deduced the formula of carrier as one part of system dynamics formula. The dynamic model of herringbone gear planetary transmission system is established to analyze the dynamic characteristics of the system. The theoretical dynamic characteristic of system is verified through simulation.
2. Dynamics Theories and Equations
2.1. Dynamic Model of Planetary Transmission System
The three-dimensional model of herringbone planetary transmission system is shown in Figure 1.
(a) 3D model with carrier of system
(b) 3D model without carrier of system
The impact of the input and output shafts is not considered in this paper. Thus, the input torque is directly acting on carrier. Because each single-helical gear, respectively, has the translational motion along the direction of x, y, and z axes and turning motion in the direction of the z axis, each herringbone gear, which is coupled by 2 helical gear with opposite direction, has 8 free degrees. The carrier has the translational motion along the direction of x, y, and z axes and turning motion in the direction of the z axis. Thus, the carrier has 4 free degrees. Because the herringbone gear is divided into left and right parts, subscript L means the left side gear while subscript R means the right-side gear. The transmission system has free degrees, and the generalized coordinate is
where x, y, z, respectively, are lateral, longitudinal, and axial displacement of a component and θ is the angle displacement of the component. is the ith planetary gear.
According to the 3D model of herringbone planetary transmission system, the dynamic model of herringbone planetary transmission system is established, as shown in Figure 2.
OXY is the fixed coordinate system, and oxy is the moving reference system.
To simplify the angle displacement as the equivalent line displacement on the corresponding base circle, where means the left side and right side of the component. means the jth planetary gear. rbs, rbpj, rbr, respectively, are the radius of base circle of sun gear, the jth planetary gear, and ring gear. rc is the distance from the center of planetary gear to gyration center of carrier,where rs, rp, respectively, are the radius of pitch circle of sun gear and planetary gear.
2.2. Equivalent Displacement
Figure 3 shows the projection relation of each component of the planetary gear system on the gear face.
According to Figure 3, projection of relative displacement of each component in the meshing line or coordinate system can be calculated. The projection of relative displacement of sun gear and planetary gear along the direction of meshing line isThe projection of relative displacement of planetary gear and ring gear along the direction of meshing line isAnd the projection of relative displacement of planetary gear and carrier along the directions of xc, yc, zc, and uc iswhere Φpi is the position angle formed by the ith planetary gear relative to the first planetary gear, , αsp is the end face meshing angle of sun-planetary gear meshing, and αrp is the end face meshing angle of planetary-ring gear meshing. esp is the comprehensive meshing error in the sun-planetary gear meshing line, erp is the comprehensive meshing error in the planetary-ring gear meshing line, L is the left side, and R is the right side.
2.3. Equilibrium Differential Equation of System Dynamics
According to the Newton second law and system dynamics, the differential equation of the dynamics of sun gear can be written as (6) and (7). Equation (6) is the dynamics equation of the left side of sun gear and (7) is the dynamics equation of the right side of sun gear.where ms is the equivalent mass of sun gear, Ks is the radial support stiffness of sun gear, Cs is the radial support damping of sun gear, Ksa is the axial support stiffness, Csa is the axial support damping, Ks12 is the axial stiffness between two helical gears, Cs12 is the axial damping between two helical gears, K1t is the left-right coupling torsional stiffness, C1t is the left-right coupling torsional damping, and Tin is the input torque.
The axial support stiffness and support damping of sun gear are as follows:
The differential equation of the dynamics of planetary gear can be written as (10) and (11). Equation (10) is the dynamics equation of the left side of planetary gear and (11) is the dynamics equation of the right side of planetary gear.where mp is the equivalent mass of planetary gear, Kp is the radial support stiffness of planetary gear, Cp is the radial support damping of planetary gear, Kpa is the axial support stiffness, Cpa is the axial support damping, Kp12 is the axial stiffness between two helical gears, Cp12 is the axial damping between two helical gears, K2t is the left-right coupling torsional stiffness, and C2t is the left-right coupling torsional damping.
The axial support stiffness and support damping of planetary gear are as follows:
The differential equation of the dynamics of ring gear can be written as (13) and (14). Equation (13) is the dynamics equation of the left side of ring gear and (14) is the dynamics equation of the right side of ring gear.where mr is the equivalent mass of planetary gear, Kr is the radial support stiffness of planetary gear, Cr is the radial support damping of planetary gear, Kra is the axial support stiffness, Cra is the axial support damping, Kr12 is the axial stiffness between two helical gears, Cr12 is the axial damping between two helical gears, K3t is the left-right coupling torsional stiffness, and C3t is the left-right coupling torsional damping.
The axial support stiffness and support damping of ring gear are as follows:
Figure 4 is the dynamics model of carrier, which combined with the Newton second law can deduct the differential equation of the carrier dynamics, as shown in where mc is the equivalent mass of carrier, Kc is the radial support stiffness of carrier, Cc is the radial support damping of carrier, Kca is the axial support stiffness, Cca is the axial support damping, Kct is the torsional stiffness, and Cct is the torsional damping.
The axial support stiffness and support damping of carrier is as follows:
All above equations can be arranged as matrix form as (18):where M is the generalized mass matrix, X is the generalized coordinate displacement matrix, C is the damping matrix, K is the stiffness matrix, and T is the external excitation column matrix.
2.4. Meshing Error Excitation
Short-period error is mainly composed by base circle error and involute tooth profile error. In the filed of gear transmission, the product of the rotation frequency of the shaft and the number of teeth, which is named meshing frequency, is noted as .where ns, np, nr, respectively, are the revolution speed of sun gear, planetary gear, and ring gear and zs, zp, zr, respectively, are the teeth number of sun gear, planetary gear, and ring gear.
The circular frequency corresponding to the meshing frequency is .
There will be errors in the manufacture and installation of gears due to objective factors, which is the main factor to make vibration of gear system. The variation rule of short period error can be expressed by Fourier series, which is based on the fundamental frequency of the wave. There is a total of 4N periodic excitation.
Error excitation of inner meshing isAnd error excitation of outer meshing iswhere is the 1st order amplitude of the transfer error between the sun gear and the planetary gear, is the initial phase of the transfer error between the sun gear and the planetary gear, is the 1st order amplitude of the transfer error between the sun gear and the planetary gear, is the initial phase of the transfer error between the sun gear and the planetary gear, is the phase difference between the meshing of ring gear and planetary gear and the meshing of sun gear and planetary gear, and is the stagger angle of herringbone gear.
3. Impact of Stagger Angle on Dynamic Characteristic of System
There is a wind power gearbox whose main parameters are shown in Table 1. And the initial amplitude and phase of gear transmission error can be calculated in Kisssoft. The numerical value of first 3-order harmonics amplitude, shown in Table 2, is needed because the numerical value of higher order harmonics amplitude is very small.
Because of the special tooth shape of herringbone gear, there will be formed unavoidable stagger angle in the processing of gear manufacturing. As Figure 5 shown, left side is the plan sketch of stagger angle while the right side is the 3D model of it.
(a) Sketch of stagger angle is 0
(b) Sketch of stagger angle is
While stagger angle, respectively, is 0, π/4, π/2, 3π/4, and π, the corresponding maximum meshing force can be calculated. Figure 6 shows the maximum meshing force curve.
(a) Inner meshing
(b) Outer meshing
According to the equation of load sharing coefficient, the load sharing coefficient of a single gear pair can be defined aswhere n is the period number of meshing frequency, N is the number of planetary gear, and F is the corresponding meshing force of the meshing cycle.
Likewise, the load sharing coefficient of herringbone gear inner and outer meshing is
The load sharing coefficient of inner and outer meshing can be calculated by (23), (24), and Figure 6, as shown in Figure 7.
(a) Inner meshing
(b) Outer meshing
In summary, the change of stagger angle has obvious influence on the dynamic characteristics of the left side of the herringbone gear but has little effect on the right side. In the processing of stagger angle increase from 0 to π, both the left and right side maximum meshing force have an unobvious change. But the load sharing coefficient of left side decreases at first and increases at last, while the right side decreases all the time. The load sharing coefficient of inner meshing meets the minimum when stagger angle is π/2, and the outer meshing meets the minimum when stagger angle is 3π/4. Overall factors which are comprehensively considered, the best stagger angle of this system is 3π/4.
4. Impact of Support Stiffness on Dynamic Characteristic of System
4.1. Theoretical Dynamics Analysis on Initial Parameter
To analyze the impact of support stiffness on dynamic characteristics of the herringbone gear planetary transmission system from the view of practical application, the parameters researched in this paper were gotten from the wind power gearbox described. In the transmission system, all types of the initial support stiffness of sun gear, planetary gear, and ring gear are 1×107N/m, and the stagger angle is 0. All of these parameters are substituted in the dynamics program; the meshing force is calculated and shown in Figure 8.
(a) Left-side meshing force of inner mesh
(b) Right-side meshing force of inner mesh
(c) Left-side meshing force of outer mesh
(d) Right-side meshing force of outer mesh
4.2. Theoretical Dynamics Analysis on Changing Support Stiffness
Keep the stiffness of planetary gear and ring gear be a constant value which is 107N/m, and change the stiffness of sun gear from 106N/m to 108N/m. The stiffness, respectively, is 106, 5×106, 107, 5×107, and 108N/m, and the corresponding meshing force can be gotten and be coupled as 3D sketch as shown in Figures 9–14.
(a) 3D curve of left-side meshing force
(b) 3D curve of right-side meshing force
(a) Left-side meshing force 3D curve
(b) Right-side meshing force 3D curve
(a) 3D curve of left-side meshing force
(b) 3D curve of right-side meshing force
(a) 3D curve of left-side meshing force
(b) 3D curve of right-side meshing force
(a) 3D curve of left-side meshing force
(b) 3D curve of right-side meshing force
(a) 3D curve of left-side meshing force
(b) 3D curve of right-side meshing force
The maximum meshing force of corresponding support stiffness of sun gear can be coupled as shown in Figure 15. Therefore, the load sharing coefficient curve is shown in Figure 16.
(a) Inner meshing maximum meshing force curve
(b) Outer meshing maximum meshing force curve
(a) Inner meshing load sharing coefficient curve
(b) Outer meshing load sharing coefficient curve
Likewise, keep the stiffness of sun gear and ring gear be a constant value which are 107N/m, and change the stiffness of planetary gear from 106N/m to 108N/m. The stiffness is 106, 5×106, 107, 5×107, and 108N/m, respectively. And the corresponding maximum meshing force and load sharing coefficient can be gotten as shown in Figures 17 and 18. Keep the stiffness of sun gear and planetary gear be a constant value which are 107N/m, and change the stiffness of ring gear from 106N/m to 108N/m. The corresponding maximum meshing force and load sharing coefficient can be shown in Figures 19 and 20.
(a) Inner meshing maximum meshing force curve
(b) Outer meshing maximum meshing force curve
(a) Inner meshing load sharing coefficient curve
(b) Outer meshing load sharing coefficient curve
(a) Inner meshing maximum meshing force curve
(b) Outer meshing maximum meshing force curve
(a) Inner meshing load sharing coefficient curve
(b) Outer meshing load sharing coefficient curve
And change the support stiffness of carrier while the support stiffness of gears is a constant value, as shown in Figures 21 and 22. According to the research results, commendatory support stiffness of the studied system is as follows: support stiffness of sun gear 107N/m, support stiffness of planetary gear 107N/m, support stiffness of ring gear 107N/m, and support stiffness of carrier 107N/m. Curves of meshing fore and load sharing coefficient with commendatory support stiffness are shown in Figures 23 and 24, where load sharing coefficient is written as LSC.
(a) Inner meshing maximum meshing force curve
(b) Outer meshing maximum meshing force curve
(a) Inner meshing load sharing coefficient curve
(b) Outer meshing load sharing coefficient curve
(a) The left-side meshing force of inner mesh
(b) The right-side meshing force of inner mesh
(c) The left-side meshing force of outer mesh
(d) The right-side meshing force of outer mesh
(a) The left-side LSC of inner mesh
(b) The right-side LSC of inner mesh
(c) The left-side LSC of outer mesh
(d) The right-side LSC of outer mesh
In summary, firstly, with the increase of the support stiffness of sun gear, both of the maximum meshing force and load sharing coefficient of inner meshing increase while these of outer meshing decrease. Secondly, with the increase of support stiffness of planetary gear, both of the maximum meshing force and load sharing coefficient of inner meshing increase; the load sharing coefficient of outer meshing increases while the maximum meshing force of outer meshing fluctuates near the support stiffness of 107N/m. Thirdly, with the increase of support stiffness of ring gear increase, both of the maximum meshing force and load sharing coefficient of inner meshing decrease while these of outer meshing increase. Finally, there is mere change of the maximum meshing force and load sharing coefficient of inner meshing and outer meshing, no matter how much the support stiffness of carrier changes.
5. Impact of Torsional Stiffness on Dynamic Characteristic of System
5.1. Impact of Coupling Torsional Stiffness of Gears on the Dynamic Characteristics
Suppose the support stiffness of system is the commendatory value, the torsional stiffness of left-right helical gear is changed to study the impact of torsional stiffness on dynamic characteristic of system. The torsional stiffness of planetary gear and ring gear is a constant value that is 108 Nm/rad, while changing it of sun gear which is 107, 5×107, 108, 5×108, and 109 Nm/rad, respectively. The inner and outer meshing force and load sharing coefficient of system are shown in Figures 25 and 26.
(a) The maximum meshing force of inner meshing
(b) The maximum meshing force of outer meshing
(a) The load sharing coefficient of inner meshing
(b) The load sharing coefficient of outer meshing
Likewise, when the torsional stiffness of sun gear and ring gear is a constant value, the curves of maximum meshing force and load sharing coefficient with the changing of torsional stiffness of planetary gear from 107 Nm/rad to 109 Nm/rad are shown in Figures 27 and 28. Change the torsional stiffness of ring gear; the curves of maximum meshing force and load sharing coefficient with the changing of torsional stiffness of ring gear from 107 Nm/rad to 109 Nm/rad are shown in Figures 29 and 30.
(a) The maximum meshing force of inner meshing
(b) The maximum meshing force of outer meshing
(a) The load sharing coefficient of inner meshing
(b) The load sharing coefficient of outer meshing
(a) The maximum meshing force of inner meshing
(b) The maximum meshing force of outer meshing
(a) The load sharing coefficient of inner meshing
(b) The load sharing coefficient of outer meshing
5.2. Impact of Torsional Stiffness of Carrier on the Dynamic Characteristics
For researching impact of torsional stiffness of carrier on dynamic characteristic of system, the torsional stiffness of carrier is changed from 107 Nm/rad to 109Nm/rad, while the left-right coupling torsional stiffness of sun gear, planetary gear, and ring gear is a constant value. The curves of maximum meshing force and load sharing coefficient are shown in Figures 31 and 32.
(a) The maximum meshing force of inner meshing
(b) The maximum meshing force of outer meshing
(a) The load sharing coefficient of inner meshing
(b) The load sharing coefficient of outer meshing
In summary, the left-right coupling torsional stiffness of sun gear and ring gear has seldom impact on dynamic characteristic when it changes in a reasonable range, and the left-right coupling torsional stiffness of planetary gear has obvious impact on load sharing coefficient but has little impact on maximum meshing force. Likewise, the torsional stiffness of carrier has obvious impact on load sharing coefficient but has little impact on maximum meshing force.
6. Dynamic Simulation Verifications
To verify the validity of theoretical calculation, 3D model of the researched system is built to be simulation analyzed in commercial software ADAMS. The similar simulation curves are fitted together, as shown in Figure 33, where simulation curve of meshing force is written as SCMF.
(a) SCMF of inner meshing of left-side
(b) SCMF of inner meshing of right-side
(c) SCMF of outer meshing of left-side
(d) SCMF of outer meshing of right-side
Meanwhile, according to the simulation meshing force, (23) and (24), load sharing coefficient can be calculated. The curve of load sharing coefficient is shown in Figure 34, where LSC means load sharing coefficient and ROM, RIM, LIM, and LOM mean right-side outer meshing, right-side inner meshing, left-side inner meshing, and lest-side outer meshing, respectively.
(a) The Maximum LSC of LIM
(b) The Maximum LSC of RIM
(c) The Maximum LSC of LOM
(d) The Maximum LSC of ROM
In summary, the tendency and phase distribution of meshing force of simulation curves are similar to these curves of theoretical calculation. And the curves of load sharing coefficient of simulation are relatively consistent with these curves of theoretical calculation. Both of them verify the correctness of theoretical dynamic analysis method.
7. Conclusion
(1) Considering the impact of support stiffness and left-right coupling torsional stiffness on dynamics characteristic of herringbone gear transmission and the impact of carrier as one part of system, whose dynamic characteristic will influence the system dynamics, the dynamic equations of herringbone gear transmission system were established by lumped parameter method based on system dynamics and Lagrange equation.
(2) Study the impact of stagger angle of gear on dynamic characteristic of system. Stagger angle has obvious impact on load sharing coefficient but little on maximum meshing force. In the processing of stagger angle increase from 0 to π, both the left and right-side maximum meshing force have an unobvious change. But the load sharing coefficient of left side decreases at first and increases at last, while the right side decreases all the time. The load sharing coefficient of inner meshing meets the minimum when stagger angle is π/2, and the outer meshing meets the minimum when stagger angle is 3π/4, and the best stagger angle of this system is 3π/4.
(3) In the herringbone gear transmission system, the changing of support stiffness will have impact on dynamic characteristics. With the increase of the support stiffness of sun gear, both of the maximum meshing force and load sharing coefficient of inner meshing increase while these of outer meshing decrease. With the increase of support stiffness of planetary gear, both of the maximum meshing force and load sharing coefficient of inner meshing increase; the load sharing coefficient of outer meshing increases while the maximum meshing force of outer meshing fluctuates near the support stiffness of 107N/m. With the increase of support stiffness of ring gear increase, both of the maximum meshing force and load sharing coefficient of inner meshing decrease while these of outer meshing increase. And, there is mere change of the maximum meshing force and load sharing coefficient of inner meshing and outer meshing, no matter how much the support stiffness of carrier changes. The commendatory support stiffness of this paper is as follows: all of support stiffness of sun gear, planetary gear, ring gear, and carrier is 107N/m.
(4) The impact of torsional stiffness on system characteristic was considered. The meshing force and load sharing coefficient seldomly changed while the changing of left-right coupling torsional stiffness changed from 107N m/rad to 109Nm/rad except the change of coupling torsional stiffness of planetary gear. When the coupling torsional stiffness changed from 107N m/rad to 109Nm/rad, the maximum meshing force is nearly a constant value while the load sharing coefficient increases slowly. Moreover, with the change of torsional stiffness of carrier, the meshing force is merely changed and the load sharing coefficient changes a lot.
(5) The dynamic characteristic was simulated by commercial software ADAMS. And curves of meshing force and load sharing coefficient gotten from ADAMS simulation were compared with these of theoretical calculation. The result is shown that both the meshing force curves and the load sharing curves are in good agreement with the theoretical curves, which validates the necessity and correctness of theoretical analysis method.
Data Availability
All data generated or analyzed during this study are included in this published article, and other pieces of information are available from the corresponding author on reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research is supported by National Natural Science Foundation of China (no. 51805368), Natural Science Foundation of Tianjin (no. 17JCQNJC04300), Open Funding of the State Key Laboratory of Materials Processing and Die & Mould Technology-Huazhong University of Science and Technology (no. P2019-022), Fundamental Research Funds for the Tianjin Universities (no. 2017KJ083), Applied Basic Research Project of China Textile Industry Association (no. J201806), and the Program for Innovative Research Team in University of Tianjin (no. TD13-5037). Thanks also go to Professor Haruo Houjoh, Shigeki Matsumura of Tokyo Institute of Technology, Professor Izhak Bucher of Israel Institute of Technology, and Professor Yidu Zhang of Beihang University for their constant assistances.