Dynamic characteristic analysis of a multi-power coupling transmission system considering non-torque loads

1. Introduction

The multi-power coupling transmission system (MPCTS) is an electromechanical coupling system, which typically uses two or more induction motors to provide its driving torque. It is used in various industrial applications, such as cutter header driving systems of tunnel boring machinery, and yaw systems of wind turbines, because of its high dynamic performance, robustness, low cost, reliability, and efficiency [1, 2]. In practical engineering applications, MPCTS systems are used to provide high driving torque. However, the effects of loading from both the torque load and the non-torque, axial and radial loads, can be complex. During MPCTS operation, the non-torque load can lead to unequal loading between gears, and an increase in the dynamic meshing force among gear pairs, which negatively affect the dynamic characteristics and load-sharing properties of the system. This indicates that non-torque load is a major contributing factor in the safety and reliability of MPCTS operation, because it contributes to premature gear transmission system failures, and internal permanent damage [3]. Therefore, it is necessary to investigate the effects of non-torque load on MPCTS systems, in order to provide better theoretical guidance for system design. Researchers have recently begun to pay attention to the influence of non-torque loads on the dynamic characteristics of gear transmission systems. Tan et al. created a dynamic model of a wind turbine to investigate the effects of these loads on the performance of its drivetrain, bearing in mind the effects of non-torque loading caused by platform motion. The results suggest that the non-torque loads not only introduce additional excitation frequencies, but also increase the vibration amplitude of the drivetrain [4]. Likewise, Li et al. constructed a coupling model of a wind turbine drivetrain to analyze its dynamic characteristics under non-torque loads, caused by blade gravity, wind shear, tower shadow, and yawed inflow. Those results indicate that non-torque loads can lead to unequal planet load sharing, which are commonly carried by carrier 1 bearings [5]. Zhao et al. created a dynamic model of a helical gear system, to investigate the dynamic behaviors influenced by gear geometric eccentricity. It was found that there are sideband frequencies around the meshing frequency of gear pairs with geometric eccentricities [6]. Cao et al. built a dynamic model of a planetary gear set with four planets to analyze dynamic performance under non-torque loads caused by the gear eccentricities. This study points out that the non-torque loads caused by gear eccentricity error not only affect the dynamic response of the planetary gear set, but also deteriorate the load-sharing performance of the planetary gear set [7]. Zhai et al. studied the dynamic characteristics of one stage helical planetary gear system, under the influence of non-torque loading caused by errors in carrier assembly. The results show that those errors tend to increase the dynamic meshing forces of the gear pairs [8]. Park et al. developed two analysis models of a wind turbine gearbox, consisting of two different load configurations: one with both torque and non-torque loading, and the another with torque only loading. After investigating the influence of wind fluctuation on the load distribution of gear tooth flanks, along with the load sharing of planets, it was concluded that the influence of non-torque forces are very important to consider, to accurately determine the design load for a wind turbine gearbox [9]. Subsequently, Park et al. built a 3D analysis model of a three-point suspension gearbox for a 2-MW class wind turbine, to study the effects of axial force, radial force, and bending moment, on the load sharing of planet gears. The results show that the load sharing of planet gears is mainly affected by the non-torque factors of moment and radial force [10]. Choudhury and Tandon created a theoretical analysis model of a rotor bearing system to study its dynamic response under radial load conditions. They concluded that vibration amplitudes in both healthy and defective bearings increase with increasing speed but that those amplitudes are greater in defective bearings owing to radial load [11].

Qiu et al. developed a rotational-translational-axial dynamic model of a planetary gear transmission system, to investigate the effects of gravity, ring support stiffness, and bedplate tilt angle on the load-sharing characteristics of a planetary gear transmission system [12]. Zhu et al. created a dynamic model of a multi-floating gear system using the lumped parameter method, and analyzed the effects of a flexible pin on gear natural frequency and dynamic meshing force [13]. Wei et al. proposed a dynamic modeling method for a multistage planetary gear system by combining the lumped parameter and finite element methods (FEM). The dynamic response calculated by the proposed method was compared with the individual results calculated the individual methods. The results show that the accuracy of the proposed combination method is higher than that of the lumped parameter method, and lower than that of FEM, although the calculation efficiency is higher than that of FEM [14]. Frabco et al. developed a dynamic model of a planetary gearbox by combining the lumped parameter method and FEM, yielding very similar modal property results from simulation versus experimentation [15]. Zeng et al. built a nonlinear dynamic model of a gear-shaft-housing system by using FEM, to study the influence of the meshing and bearing support stiffness on vibration response. It was concluded that if the bearing support and meshing stiffness increase, the vibration acceleration and velocity of the housing also increase. However, the influence from the bearing stiffness is greater than that of the meshing stiffness [16]. Wang et al. created a coupled dynamic model of a spur, gear-shaft-bearing system, to investigate the effects of the shaft flexibility on bearing response and radiated noise [17].

In the aforementioned research studies, the focus is mainly on either dynamic characteristic of a system affected by non-torque loads, or a dynamic modeling method which only considers the characteristics of the gear transmission system itself. Subsequently, the effects of the time-varying excitations caused by the induction motor are ignored, although the mechanical and electric assemblies of the MPCTS cannot be divided into two non-correlated parts. Since the vibration caused by the time-varying meshing stiffness of the gear transmission system can be transmitted to the induction motor through the coupling, and the vibration caused by the speed or torque fluctuations of the motor can to be transmitted to the gear transmission system through the coupling. Therefore, in this study, an electromechanical coupling dynamic model of MPCTS with a gear transmission and induction motor, was used to analyze dynamic characteristic under non-torque loads. Section 2 presents the development of the MPCTS model used, while Section 3 details its resulting dynamic characteristics under a non-torque load condition. Section 4 presents the conclusions drawn on the basis of the analysis of the results.

2. Electromechanical coupling dynamic model for MPCTS

The MPCTS, as shown in Fig. 1, is composed of a single-stage planetary gear system and a parallel shaft gear system. In the parallel shaft gear system, the three pinions are driven separately by three induction motors, and the wheel is driven by the three pinions, which drive the wheel. The sun shaft is connected to the wheel via a spline, and the power from the wheel is transmitted to the single-stage planetary gear system through the sun shaft.

Fig. 1Three-dimensional-model of MPCTS

2.1. Dynamic model of Timoshenko beam

The dynamic model of the rotating shafts and gears was developed using the Timoshenko beam theory proposed by Nelson [18], and is shown in Fig. 2. The subscripts $pii$ and $piiw$ represent the $i$th pinion and ith gear pair formed by the $i$th pinion and wheel, respectively. The subscript $w$ represents the wheel. The mathematical model of the two-node Timoshenko beam can be written as:

1

$M_{j,j+1}{\ddot{X}}_{j,j+1}+C_{j,j+1}{\dot{X}}_{j,j+1}+K_{j,j+1}{X}_{j,j+1}=F_{j,j+1}.$

where $M_{j,j+1}$ denotes the mass matrix of the Timoshenko beam. $C_{j,j+1}$ and $K_{j,j+1}$ denote the damping matrix and stiffness matrix of the Timoshenko beam, and their detailed matrix forms are shown in Ref. [19]. $X_{j,j+1}$ and $F_{j,j+1}$ denote the displacement and exciting force column vectors of the Timoshenko beam.

Fig. 2Timsoshenko beam

2.2. Dynamic model of parallel shaft gear system

In the parallel shaft gear system, all the gears and shafts are equivalent to the Timoshenko beam. The basic parameters of the parallel shaft gear system are listed in Table 1. The dynamic model of the parallel shaft gear system can be obtained using the node finite element method, which is formed from the bearing, meshing gear pair, and shaft, as shown in Fig. 3(a). Nodes $i$ and $j$ are assumed to be the meshing nodes of the pinion-wheel gear pair. Nodes $i-2$, $i+2$, and $i+3$ are the bearing supporting points of the pinion and wheel, respectively. As shown in Eq. (2), the dynamic model of the supporting bearing can then be expressed in matrix form, using $i-2$ as the bearing supporting point in this example:

2

$C_{i-2}{\dot{X}}_{i-2}+K_{i-2}{X}_{i-2}=0,$

where $C_{i-2}$ and $K_{i-2}$ denote the damping matrix and stiffness matrix of the bearing, respectively. $X_{i-2}$ denotes the displacement vector of bearing supporting point $i-2$. The detailed expressions of the matrices can be found in Ref. [20].

Fig. 3a) Finite element node diagram of parallel shaft gear system; b) meshing pair model of parallel shaft gear system

a)

b)

The meshing pair model of the parallel shaft gear system is shown in Fig. 3 (b), where the gear meshing pair is modeled as a stiffness-damping spring system, with the spring acting along the line of action between the two gears. The deformation along the line of action can be written as:

3

${\delta }_{piiw}=-x_{pii}\mathrm{s}\mathrm{i}\mathrm{n}(\alpha -{\varphi }_{piiw})\mathrm{c}\mathrm{o}\mathrm{s}\beta +y_{piiw}\mathrm{c}\mathrm{o}\mathrm{s}(\alpha -{\varphi }_{piiw})\mathrm{c}\mathrm{o}\mathrm{s}\beta +x_w\mathrm{s}\mathrm{i}\mathrm{n}(\alpha -{\varphi }_{piiw})\mathrm{c}\mathrm{o}\mathrm{s}\beta$
$-y_w\mathrm{c}\mathrm{o}\mathrm{s}\left(\alpha -{\varphi }_{piiw}\right)\mathrm{c}\mathrm{o}\mathrm{s}\beta +r_{pii}{\theta }_{pii}\mathrm{c}\mathrm{o}\mathrm{s}\beta -r_w{\theta }_w\mathrm{c}\mathrm{o}\mathrm{s}\beta +z_{pii}\mathrm{s}\mathrm{i}\mathrm{n}\beta -z_w\mathrm{s}\mathrm{i}\mathrm{n}\beta -e_{piiw},$

where $x_j$, $y_j$ and $z_j\mathrm{}(j=pii,w)$ denote the displacement of the pinion and wheel along the $x$-, $y$-, and $z$-axes, respectively. ${\theta }_j$ and $r_j\mathrm{}(j=pii,w)$ denote the angle displacement around the $z$-axis and the base circle radius of the pinion and wheel, respectively. $\alpha$ is the pressure angle, and ${\varphi }_{piiw}$ is the position angle of the pinion. $\beta$ denotes the helix angle of the gear. $e_{piiw}$ is the transmission error, which can be expressed as:

4

$e_{piiw}={\stackrel{-}{e}}_{piiw}+\sum _{l=1}^7{\widehat{e}}_{piiwl}\cdot \mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{l\cdot 2\pi \cdot r_{pii}\cdot {\theta }_{pii}}{p_{bpii}}+{\gamma }_i\right),$

where ${\widehat{e}}_{piiw}$ and ${\stackrel{-}{e}}_{piiw}$ are the $l$th harmonic term amplitude and mean transmission error, respectively, and ${\gamma }_i$ is the $i$th pinion phase angle.

Hence, the differential equation of the pinion-wheel pair can be deduced as:

5

$\left\{\begin{array}{l}{m}_{pii}{\ddot{x}}_{pii}-c_{piiw}{\dot{\delta }}_{piiw}\mathrm{s}\mathrm{i}\mathrm{n}(\alpha -{\varphi }_{piiw})\mathrm{c}\mathrm{o}\mathrm{s}\beta -k_{piiw}{\delta }_{piiw}\mathrm{s}\mathrm{i}\mathrm{n}(\alpha -{\varphi }_{piiw})\mathrm{c}\mathrm{o}\mathrm{s}\beta =0,\\ m_{pii}{\ddot{y}}_{pii}+c_{piiw}{\dot{\delta }}_{piiw}\mathrm{c}\mathrm{o}\mathrm{s}(\alpha -{\varphi }_{piiw})\mathrm{c}\mathrm{o}\mathrm{s}\beta +k_{piiw}{\delta }_{piiw}\mathrm{c}\mathrm{o}\mathrm{s}(\alpha -{\varphi }_{piiw})\mathrm{c}\mathrm{o}\mathrm{s}\beta =0,\\ m_{pii}{\ddot{z}}_{pii}+c_{piiw}{\dot{\delta }}_{piiw}\mathrm{s}\mathrm{i}\mathrm{n}\beta +k_{piiw}{\delta }_{piiw}\mathrm{s}\mathrm{i}\mathrm{n}\beta =0,\\ J_{pii}{\ddot{\theta }}_{pii}+c_{piiw}{\dot{\delta }}_{piiw}{r}_{pii}\mathrm{c}\mathrm{o}\mathrm{s}\beta +k_{piiw}{\delta }_{piiw}{r}_{pii}\mathrm{c}\mathrm{o}\mathrm{s}\beta =T_{mi},\\ m_w{\ddot{x}}_w+c_{piiw}{\dot{\delta }}_{piiw}\mathrm{s}\mathrm{i}\mathrm{n}(\alpha -{\varphi }_{piiw})\mathrm{c}\mathrm{o}\mathrm{s}\beta +k_{piiw}{\delta }_{piiw}\mathrm{s}\mathrm{i}\mathrm{n}(\alpha -{\varphi }_{piiw})\mathrm{c}\mathrm{o}\mathrm{s}\beta =0,\\ m_w{\ddot{y}}_w-c_{piiw}{\dot{\delta }}_{piiw}\mathrm{c}\mathrm{o}\mathrm{s}(\alpha -{\varphi }_{piiw})\mathrm{c}\mathrm{o}\mathrm{s}\beta -k_{piiw}{\delta }_{piiw}\mathrm{c}\mathrm{o}\mathrm{s}(\alpha -{\varphi }_{piiw})\mathrm{c}\mathrm{o}\mathrm{s}\beta =0,\\ m_w{\ddot{z}}_w-c_{piiw}{\dot{\delta }}_{piiw}\mathrm{s}\mathrm{i}\mathrm{n}\beta -k_{piiw}{\delta }_{piiw}\mathrm{s}\mathrm{i}\mathrm{n}\beta =0,\\ J_w{\ddot{x}}_w-c_{piiw}{\dot{\delta }}_{piiw}{r}_w\mathrm{c}\mathrm{o}\mathrm{s}\beta -k_{piiw}{\delta }_{piiw}{r}_w\mathrm{c}\mathrm{o}\mathrm{s}\beta =-T_l,\end{array}\right.$

where $m_j$ and $J_j\mathrm{}(j=pii,w)$ denote the mass and inertia of the pinion and wheel, respectively. $T_{mi}$ and $T_l$ denote the input torque and load torque, while $c_{piiw}$ and $k_{piiw}$ denote the meshing damping and stiffness of the pinion-wheel pair, respectively. They can be expressed as Eq. (6):

6

$\left\{\begin{array}{l}{k}_{piiw}={\stackrel{-}{k}}_{piiw}+\sum _{l=1}^7{\widehat{k}}_{piiwl}\cdot \mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{l\cdot 2\pi \cdot r_{pii}\cdot {\theta }_{pii}}{p_{bpii}}+{\gamma }_i\right),\\ c_{piiw}=2\zeta \sqrt{\frac{k_{piiw}}{\left(1/m_{eq,1}+1/m_{eq,2}\right)}},\end{array}\right.$

where ${\widehat{k}}_{piiw}$ and ${\stackrel{-}{k}}_{piiw}$ are the $l$th harmonic term amplitude and mean mesh stiffness, respectively. $\zeta$ is the meshing damping ratio, $m_{eq,i}=I_{xi}/r_i^2(i=pii,w)$ represents the equivalent quality. Therefore, the matrix form of the pinion-wheel pair dynamic model can be written as:

7

$\left\{\begin{array}{l}{M}_{piiw}{\ddot{X}}_{piiw}+C_{piiw}{V}_{piiw}^T{\dot{\delta }}_{piiw}+K_{piiw}{V}_{piiw}^T{\delta }_{piiw}=F_{piiw},\\ {\delta }_{piiw}=V_{piiw}{X}_{piiw}-e_{piiw},\end{array}\right.$

where, $V_{piiw}$ is the meshing vector of the pinion-wheel pair, which can be written as:

8

$V_{piiw}=\left[\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Phi },\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Phi },\mathrm{s}\mathrm{i}\mathrm{n}\beta ,r_{pii}\mathrm{c}\mathrm{o}\mathrm{s}\beta ,-\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Phi },\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Phi },-\mathrm{s}\mathrm{i}\mathrm{n}\beta ,-r_w\mathrm{c}\mathrm{o}\mathrm{s}\beta \right],$

where $M_{piiw}$, $C_{piiw}$, and $K_{piiw}$ denote the mass matrix, meshing damping matrix, and meshing stiffness matrix of the pinion-wheel pair, respectively. $X_{piiw}$ and $F_{piiw}$ denote the displacement and exciting force vectors, respectively. In Eq. (8), $\mathrm{\Phi }=\alpha -{\varphi }_{piiw}$.

Table 1Basic parameters of parallel shaft gear system

	Tooth number	Module/mm	Face width/mm	Pressure angle/°	Supporting stiffness/Nm-1	Diameter/mm
Pinion	19	4	146	20	1e8	D1 = 100; D2 = 110; D3 = 160; D4 = 55; D5 = 50
Wheel	99	4	146	20	1e8

2.3. Dynamic model of single-state planetary gear system

The dynamic model of a single-stage planetary gear system can also be obtained using the node finite element method, consisting of the bearing, sun-planet meshing pair, planet-ring meshing pair, pin, and carrier. The carrier, ring, planet, sun, and pin are modeled as different cross-section types of Timoshenko beams, and can be divided into different element nodes, as shown in Fig. 4(a). The basic parameters of the planetary stage are listed in Table 2. Herein, the dynamic model of the supporting bearing in a single-stage planetary gear system can be expressed using Eq. (2) as well.

Table 2Basic parameters of single-stage planetary gear system

	Tooth number	Module / mm	Face width / mm	Pressure angle /°	Helix angle /°	Supporting stiffness / Nm-1	Diameter/mm
Sun	18	6	146	20	25.5	1e8	D6 =260; D7 = 174.6; D8 = 268; D9 = 320
Planet	52	6	146	20	25.5	1e9
Ring	123	6	146	20	25.5	1e10

The meshing pair model of the single-stage planetary gear system is shown in Fig. 4(b). The subscripts $s$, $r$, $c$, and $p_i$ represent the sun, ring, carrier, and $i$th planet, respectively. There are two types of meshing pairs existing in the model: an external, sun-planet meshing pair, and an internal, planet-ring meshing pair. Without loss of generality, the external and internal meshing pairs are also modeled as a stiffness-damping spring system, with the spring acting along the line of action between the two gears.

For the sun-planet $i$ meshing pair, the deformation along the line of action can be written as:

9

${\delta }_{spi}=x_s\sin \left({\alpha }_t+{\varphi }_n\right)\mathrm{c}\mathrm{o}\mathrm{s}\beta -y_s\cos \left({\alpha }_t+{\varphi }_n\right)\mathrm{c}\mathrm{o}\mathrm{s}\beta -x_{pi}\sin \left({\alpha }_t\right)\mathrm{c}\mathrm{o}\mathrm{s}\beta$
$+y_{pi}\mathrm{c}\mathrm{o}\mathrm{s}\left({\alpha }_t\right)\mathrm{c}\mathrm{o}\mathrm{s}\beta -r_s{\theta }_s\mathrm{c}\mathrm{o}\mathrm{s}\beta -r_{pi}{\theta }_{pi}\mathrm{c}\mathrm{o}\mathrm{s}\beta +z_{pi}\mathrm{s}\mathrm{i}\mathrm{n}\beta -z_s\mathrm{s}\mathrm{i}\mathrm{n}\beta -e_{spi},$

where $x_j$, $y_j$ and $z_j\mathrm{}(j=s,pi)$ denote the displacement of the sun and planet $i$ along the $x$-, $y$-, and $z$-axes, respectively. ${\theta }_j$ and $r_j\mathrm{}(j=s,pi)$ denote the angle displacement around the $z$-axis and the base circle radius of the sun and planet $i$, respectively. ${\alpha }_t$ represents the pressure angle, and ${\varphi }_n$ represents the position angle of the planet. $e_{spi}$ is the transmission error of the sun-planet $i$ meshing pair, which can be expressed as Eq. (15).

The differential equation of the sun-planet $i$ meshing pair can be deduced as follows:

10

$\left\{\begin{array}{l}{m}_s{\ddot{x}}_s+c_{spi}{\dot{\delta }}_{spi}\mathrm{s}\mathrm{i}\mathrm{n}({\alpha }_t+{\varphi }_n)\mathrm{c}\mathrm{o}\mathrm{s}\beta +k_{spi}{\delta }_{spi}\mathrm{s}\mathrm{i}\mathrm{n}({\alpha }_t+{\varphi }_n)\mathrm{c}\mathrm{o}\mathrm{s}\beta =0,\\ m_s{\ddot{y}}_s-c_{spi}{\dot{\delta }}_{spi}\mathrm{c}\mathrm{o}\mathrm{s}({\alpha }_t+{\varphi }_n)\mathrm{c}\mathrm{o}\mathrm{s}\beta -k_{spi}{\delta }_{spi}\mathrm{c}\mathrm{o}\mathrm{s}({\alpha }_t+{\varphi }_n)\mathrm{c}\mathrm{o}\mathrm{s}\beta =0,\\ m_s{\ddot{z}}_s-c_{spi}{\dot{\delta }}_{spi}\mathrm{s}\mathrm{i}\mathrm{n}\beta -k_{spi}{\delta }_{spi}\mathrm{s}\mathrm{i}\mathrm{n}\beta =0,\\ J_s{\ddot{\theta }}_s-c_{spi}{\dot{\delta }}_{spi}{r}_s\mathrm{c}\mathrm{o}\mathrm{s}\beta -k_{spi}{\delta }_{spi}{r}_s\mathrm{c}\mathrm{o}\mathrm{s}\beta =T_{in},\\ m_{pi}{\ddot{x}}_{pi}-c_{spi}{\ddot{\delta }}_{spi}\mathrm{s}\mathrm{i}\mathrm{n}\left({\alpha }_t\right)\mathrm{c}\mathrm{o}\mathrm{s}\beta -k_{spi}{\delta }_{spi}\mathrm{s}\mathrm{i}\mathrm{n}\left({\alpha }_t\right)\mathrm{c}\mathrm{o}\mathrm{s}\beta =0,\\ m_{pi}{\ddot{y}}_{pi}+c_{spi}{\dot{\delta }}_{spi}\mathrm{c}\mathrm{o}\mathrm{s}\left({\alpha }_t\right)\mathrm{c}\mathrm{o}\mathrm{s}\beta +k_{spi}{\delta }_{spi}\mathrm{c}\mathrm{o}\mathrm{s}\left({\alpha }_t\right)\mathrm{c}\mathrm{o}\mathrm{s}\beta =0,\\ m_{pi}{\ddot{z}}_{pi}+c_{spi}{\dot{\delta }}_{spi}\mathrm{s}\mathrm{i}\mathrm{n}\beta +k_{spi}{\delta }_{spi}\mathrm{s}\mathrm{i}\mathrm{n}\beta =0,\\ J_{pi}{\ddot{\theta }}_{pi}-c_{spi}{\dot{\delta }}_{spi}{r}_{pi}\mathrm{c}\mathrm{o}\mathrm{s}\beta -k_{spi}{\delta }_{spi}{r}_{pi}\mathrm{c}\mathrm{o}\mathrm{s}\beta =0,\end{array}\right.$

where $m_j$ and $J_j\mathrm{}(j=s,pi)$ denote the mass and inertia of the sun and planet $i$, respectively. $c_{spi}$ and $k_{spi}$ denote the meshing damping and stiffness of the sun-planet $i$ meshing pair, and can be expressed as Eq. (15). $T_{in}$denotes the input torque of the sun. The matrix form of the dynamic model of the sun-planet $i$ meshing pair can be written as:

11

${\delta }_{spi}=V_{spi}{X}_{spi}-e_{spi}\mathrm{a}\mathrm{n}\mathrm{d}{M}_{spi}{\ddot{X}}_{spi}+C_{spi}{V}_{spi}^T{\dot{\delta }}_{spi}+K_{spi}{V}_{spi}^T{\delta }_{spi}=F_{spi},$

where $M_{spi}$, $C_{spi}$, and $K_{spi}$ denote the mass matrix, meshing damping matrix, and meshing stiffness matrix of the sun-planet $i$ pair, respectively. $X_{spi}$ and $F_{spi}$ denote the displacement and excitation force vectors, respectively. Moreover, $V_{spi}$ is the meshing vector of the sun-planet $i$ meshing pair, which can be written as:

12

$V_{spi}=\left[\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Psi },-\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Psi },-\mathrm{s}\mathrm{i}\mathrm{n}\beta ,-r_s\mathrm{c}\mathrm{o}\mathrm{s}\beta ,-\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_t,\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_t,\mathrm{s}\mathrm{i}\mathrm{n}\beta ,-r_{pi}\mathrm{c}\mathrm{o}\mathrm{s}\beta \right],$

where $\mathrm{\Psi }={\alpha }_t+{\varphi }_n$.

Fig. 4a) dynamic model of single-stage planetary gear system; b) meshing pair model of single-stage planetary gear system

a)

b)

For the planet $i$-ring meshing pair, the deformation along the line of action can be written as:

13

${\delta }_{pir}=x_{pi}\mathrm{s}\mathrm{i}\mathrm{n}\left({\alpha }_t\right)\mathrm{c}\mathrm{o}\mathrm{s}\beta +y_{pi}\mathrm{c}\mathrm{o}\mathrm{s}\left({\alpha }_t\right)\mathrm{c}\mathrm{o}\mathrm{s}\beta -x_r\mathrm{s}\mathrm{i}\mathrm{n}({\varphi }_n-{\alpha }_t)\mathrm{c}\mathrm{o}\mathrm{s}\beta$
$-y_r\mathrm{c}\mathrm{o}\mathrm{s}\left({\varphi }_n-{\alpha }_t\right)\mathrm{c}\mathrm{o}\mathrm{s}\beta +r_{pi}{\theta }_{pi}\mathrm{c}\mathrm{o}\mathrm{s}\beta -r_r{\theta }_r\mathrm{c}\mathrm{o}\mathrm{s}\beta +z_{pi}\mathrm{s}\mathrm{i}\mathrm{n}\beta -z_r\mathrm{s}\mathrm{i}\mathrm{n}\beta -e_{pir},$

where $x_r$, $y_r$ and $z_r$ denote the displacement of the ring along the $x$-, $y$-, and $z$-axes, respectively. ${\theta }_r$ and $r_r$ denote the angle displacement around the $z$-axis and the base circle radius of the ring, respectively. $e_{pir}$ is the transmission error of the planet $i$-ring meshing pair, which can be expressed as Eq. (15).

The differential equation of the planet $i$-ring meshing pair can be deduced as follows:

14

$\left\{\begin{array}{l}{m}_{pi}{\ddot{x}}_{pi}+c_{pir}{\dot{\delta }}_{pir}\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_t+k_{pir}{\delta }_{pir}\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_t=0,\\ m_{pi}{\ddot{y}}_{pi}+c_{pir}{\dot{\delta }}_{pir}\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_t+k_{pir}{\delta }_{pir}\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_t=0,\\ m_{pi}{\ddot{z}}_{pi}+c_{pir}{\dot{\delta }}_{pir}\mathrm{s}\mathrm{i}\mathrm{n}\beta +k_{pir}{\delta }_{pir}\mathrm{s}\mathrm{i}\mathrm{n}\beta =0,\\ J_{pi}{\ddot{\theta }}_{pi}+c_{pir}{\dot{\delta }}_{pir}{r}_{pi}\mathrm{c}\mathrm{o}\mathrm{s}\beta +k_{pir}{\delta }_{pir}{r}_{pi}\mathrm{c}\mathrm{o}\mathrm{s}\beta =0,\\ m_r{\ddot{x}}_r-c_{pir}{\dot{\delta }}_{pir}\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}({\varphi }_n-{\alpha }_t)-k_{pir}{\delta }_{pir}\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}({\varphi }_n-{\alpha }_t)=0,\\ m_r{\ddot{y}}_r-c_{pir}{\dot{\delta }}_{pir}\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}({\varphi }_n-{\alpha }_t)-k_{pir}{\delta }_{pir}\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}({\varphi }_n-{\alpha }_t)=0,\\ m_r{\ddot{z}}_r-c_{pir}{\dot{\delta }}_{pir}\mathrm{s}\mathrm{i}\mathrm{n}\beta -k_{pir}{\delta }_{pir}\mathrm{s}\mathrm{i}\mathrm{n}\beta =0,\\ J_r{\ddot{\theta }}_r-c_{pir}{\dot{\delta }}_{pir}{r}_r\mathrm{c}\mathrm{o}\mathrm{s}\beta -k_{pir}{\delta }_{pir}{r}_r\mathrm{c}\mathrm{o}\mathrm{s}\beta =0,\end{array}\right.$

where $m_r$ and $J_r$ denote the mass and inertia of the ring, respectively. $c_{pir}$ and $k_{pir}$ denote the meshing damping and stiffness of the planet $i$-ring meshing pair, respectively, and can be expressed as Eq. (15):

15

$\left\{\begin{array}{l}{k}_{spi}={\stackrel{-}{k}}_{spi}+\sum _{l=1}^7{\widehat{k}}_{spil}\cdot \mathrm{s}\mathrm{i}\mathrm{n}\left[\left(\frac{l\cdot 2\pi \cdot r_s\cdot {\theta }_s}{p_{bs}}+Z_s{\psi }_m+\mathrm{\Gamma }\right)+{\varphi }_{spi}\right],\\ e_{spi}={\stackrel{-}{e}}_{spi}+\sum _{l=1}^7{\widehat{e}}_{spil}\cdot \mathrm{s}\mathrm{i}\mathrm{n}\left[\left(\frac{l\cdot 2\pi \cdot r_s\cdot {\theta }_s}{p_{bs}}+Z_s{\psi }_m\right)+{\varphi }_{spi}\right],\\ k_{pir}={\stackrel{-}{k}}_{pir}+\sum _{l=1}^7{\widehat{k}}_{pirl}\cdot \mathrm{s}\mathrm{i}\mathrm{n}\left[\left(\frac{l\cdot 2\pi \cdot r_{pi}\cdot {\theta }_{pi}}{p_{bpi}}+Z_r{\psi }_n+\mathrm{\Gamma }+{\gamma }_{sr}\right)+{\varphi }_{pir}\right],\\ e_{pir}={\stackrel{-}{e}}_{pir}+\sum _{l=1}^7{\widehat{e}}_{pirl}\cdot \mathrm{s}\mathrm{i}\mathrm{n}\left[\left(\frac{l\cdot 2\pi \cdot r_{pi}\cdot {\theta }_{pi}}{p_{bpi}}+Z_r{\psi }_n+{\gamma }_{sr}\right)+{\varphi }_{pir}\right],\\ c_j=2\zeta \sqrt{\frac{k_j}{\left(1/m_{eq,1}+1/m_{eq,2}\right)}},j=spi,pir,\end{array}\right.$

${\widehat{e}}_i$ (or ${\widehat{k}}_i$, $i=spil$ and $pirl$) and ${\stackrel{-}{e}}_j$ (or ${\stackrel{-}{k}}_j$, $j=spil$ and $pirl$) represent the $l$th harmonic term amplitude and mean transmission error (or mean mesh stiffness), respectively. $Z_r$ and $Z_s$ represent the number of teeth in the ring and sun, respectively. ${\gamma }_{sr}$ represents the phase angle between planet $i$ (or ring) and sun (or planet $i$) meshing pair. $p_{bs}$ and $p_{bpi}$ are the base pitch of the sun and planet $i$, respectively. $\mathrm{\Gamma }$ represents the phase difference between $e_j$ and $k_j$. $z_r{\varphi }_n$ and $z_s{\psi }_m$ are the phase relationships between the meshing pairs. $m_{eq,i}=I_{xi}/r_i^2(i=s,pi,r)$ represents the equivalent quality. The matrix form of the dynamic model of the planet $i$-ring meshing pair can be written as follows:

16

${\delta }_{pir}=V_{pir}{X}_{pir}-e_{pir}\mathrm{a}\mathrm{n}\mathrm{d}{M}_{pir}{\ddot{X}}_{pir}+C_{pir}{V}_{pir}^T{\dot{\delta }}_{pir}+K_{pir}{V}_{pir}^T{\delta }_{pir}=F_{pir},$

where $M_{pir}$, $C_{pir}$, and $K_{pir}$ denote the mass matrix, meshing damping matrix, and meshing stiffness matrix of the planet $i$-ring meshing pair, respectively. $X_{pir}$ and $F_{pir}$ denote the displacement and excitation force vectors, respectively. Further, $V_{pir}$ is the meshing vector of the planet $i$-ring pair, which can be expressed as follows:

17

$V_{pir}=\left[\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_t,\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_t,\mathrm{s}\mathrm{i}\mathrm{n}\beta ,r_{pi}\mathrm{c}\mathrm{o}\mathrm{s}\beta ,-\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Theta },-\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Theta },-\mathrm{s}\mathrm{i}\mathrm{n}\beta ,-r_r\mathrm{c}\mathrm{o}\mathrm{s}\beta \right],$

where $\mathrm{\Theta }={\varphi }_n-{\alpha }_t$.

If it is assumed that the number of planets is $n$ in the planetary gear system (Fig. 5), it can be concluded that the number of pins and connections are also $n$, respectively, according to the structural features of the planetary gear system. The stiffness-damping spring system is used to represent the coupled relationship between the side plate, pin, and connecting plate; $c_{c1}$ and $k_{c1}$ represent the coupling damping and stiffness between the side plate and pin, while $c_{c2}$ and $k_{c2}$ represent the coupling damping and stiffness between the side plate and connecting plate. $c_i$ and $k_i$$(i=pin,cp)$ represent the damping and stiffness of the pin and connecting plate, respectively.

The matrix form of the dynamic model of the carrier pin can be expressed as:

18

$M_c{\ddot{X}}_c+C_c{\dot{X}}_c+K_c{X}_c=F_c,$

where $M_c$, $C_c$ and $K_c$ denote the mass matrix, damping matrix, and stiffness matrix of the carrier pin, respectively. $X_c$ and $F_c$ denote the displacement and exciting force vector of the carrier pin, respectively, and can be obtained using Eq. (13) and (14) in Ref. [19].

Fig. 5Dynamic model of carrier pin

2.5. Electromechanical coupling dynamic model

Using the node finite element method and the structural features of MPCTS, the system can be divided into 88 nodes with 352 degrees of freedom. The entire stiffness matrix of the MPCTS can be obtained using the coupled method described in Ref. [20], and the distribution of nodes illustrated in Fig. 6; which shows the stiffness of the Timoshenko beam supporting the bearing and gear meshing.

Fig. 6Dynamic model of MPCTS

The MPCTS dynamic model can be developed after obtaining the mass matrix, stiffness matrix, damping matrix, and exciting force vector of each sub-element. The components of each sub-element matrix can be superimposed on the corresponding position of the total matrix [20], based upon the total number of system nodes, and the relationship between each node. In addition, the input torques $T_{mi}(i=1,2,3)$ act on nodes 1, 11, and 21, respectively, as the driving torques of the system. The radial loads $F_X$, $F_Y$, and the torque load $T_L$ act on node 87. The master-slave control method was used to control the three induction motors, and the torque signal was selected as the control signal between the induction motors. Finally, the electromechanical coupling dynamic model of the MPCTS was be deduced as:

21

$M\ddot{X}+C\dot{X}+KX=F,$

where $M$, $C$, and $K$ denote the mass, damping, and stiffness matrices of the MPCTS, respectively. $X$ and $F$ denote the displacement and exciting force column vectors of the MPCTS, respectively.

3. Result and discussion

Two load conditions were selected to investigate the influence of non-torque load on the dynamic characteristics of MPCTS: (a) the torque load $T_L$, and (b) the non-torque, radial loads of $F_X$ and $F_Y$ in addition to the torque load $T_L$, which is shown in Fig. 7.

The motion trajectories of the pinions and the wheel are shown in Fig. 8. Geometrically, the motion trajectories of the pinions are linear, as shown in Fig. 8(a) and (b). Marked changes were observed in the angle ($\gamma$) between trajectories and the length ($L$) of the trajectories. The difference between the angles is small under the torque load condition, where ${\gamma }_1=$ 120.0020°, ${\gamma }_2=$ 120.0007°, and ${\gamma }_3=$ 119.9973°, and the maximum angle difference is between ${\gamma }_1$ and ${\gamma }_3$, namely, ${\gamma }_1-{\gamma }_3=$ 0.0047°. A similar phenomenon also occurred in the length of trajectory under torque load condition, where $L_1=$ 729.17 μm, $L_2=$ 732.25 μm and $L_3=$ 749.50 μm. The maximum length difference was between $L_1$ and $L_3$, namely, $L_1-L_3=$ 20.33 μm, as shown in Fig. 8(a). With regards to the effects of non-torque load, both the angle and the length vary significantly, but the length varied more, as shown in Fig. 8(b). Owing to the radial loads $F_X$ and $F_Y\text{,}$ the difference between the angles was further increased, where ${\gamma }_1=$ 120.7559°, ${\gamma }_2=$ 120.2356°, and ${\gamma }_3=$ 119.0085°. Notably, the maximum difference was still betwee n ${\gamma }_1$ and ${\gamma }_3$, which is ${\gamma }_1-{\gamma }_3=$ 1.7474°. Meanwhile, the lengths of trajectories increased to $L_1=$1000 μm, $L_2=$ 1200 μm and $L_3=$ 791.44 μm, and increased by 270.89 μm, 467.75 μm and 41.94 μm, respectively, as compared to the torque load effects. However, some difference occurred when the maximum length difference was between $L_2$ and $L_3\text{;}$ reaching $L_2-L_3=$ 408.56 μm.

Fig. 7Load: a) torque load; b) non-torque load

a)

b)

As shown in Fig. 8(c) and (d), the shape of the trajectory was close to a circle under either the torque load condition or the non-torque load condition. When non-torque load condition was considered, the vibration equilibrium positions of the wheel in $X$ and $Y$ directions went from –0.83 μm to 1.22 μm, and from 3.80 μm to 6.02 μm – an increase of 2.47 and 0.58 times, respectively. Moreover, the load-sharing characteristics of the pinion-wheel gear pairs also changed significantly because of the change in the vibration equilibrium position, as shown in Fig. 9.

Fig. 8Motion trajectories: a) and b) motion trajectories of the pinions under torque load and non-torque load (the blue solid line, red dotted line and green dash-dotted line are the motion trajectories of pinion 1, pinion 2 and pinion 3, respectively); c) and d) motion trajectories of the wheel under torque load and non-torque load, respectively

a)

b)

c)

d)

In this study, the maximum load-sharing coefficient (LSC) of the gear pair was selected to represent the LSC of the corresponding gear pair [23]. The LSC of the pinion-wheel gear pairs was 1.072 under torque load conditions, and 1.179 under non-torque load conditions – an increase of 9.98 %. In summary, non-torque load has a significant influence on both the vibration equilibrium positions of the pinions and the wheel, and the load-sharing characteristics of the pinion-wheel gear pairs.

The motion trajectory of the sun was analogous to a circle, as shown in Fig. 10(a) and (b). Significant changes occurred in the vibration equilibrium position of the sun. When the non-torque load wass considered, the vibration equilibrium positions of it in both $X$ and $Y$ directions were markedly offset, as compared with the torque load condition. They changed from –0.04 μm to 3.69 μm, and from 0.25 μm to 14.74 μm – increases of 93.25 and 57.96 times, respectively. In contrast, the motion trajectories of the planets were analogous to an ellipse under torque load condition, and their offsets of vibration equilibrium positions in $X$ and $Y$ directions were about 0.20 μm and –10.52 μm, respectively, as shown in Fig. 10(c). In Fig. 10(d), there was an obvious difference, as the shapes of the motion trajectories of the planets were analogous to a rectangle under non-torque load condition, and their vibration equilibrium positions in the $X$ direction were offset in two different directions. For planet 3, the vibration equilibrium position was shifted in the negative direction of the $X$-axis, with an offset of approximately $\mathrm{\Delta }{X}_3=$ 24.46 μm; increasing by 123.42 times. However, the vibration equilibrium positions for the other two planets were shifted in the positive direction of the $X$-axis, and their offsets were $\mathrm{\Delta }{X}_1=$ 20.33 μm and $\mathrm{\Delta }{X}_2=$ 11.25 μm – increasing by 63.78 and 56.78 times, respectively. In addition, the vibration equilibrium positions of the planets in the $Y$ direction were shifted in the negative direction of the $Y$-axis and change from approximately –10.52 μm to –163.66 μm – increasing by 14.56 times.

Fig. 9LSC: a) LSC of pinion-wheel gear pairs under torque load; b) LSC of pinion-wheel gear pairs under non-torque

a)

b)

Fig. 10Motion trajectories: a) and b) motion trajectories of the sun under torque load and non-torque load, respectively; c) and d) motion trajectories of the planets under torque load and non-torque load (the blue solid line, red dotted line and green dash-dotted line are the motion trajectories of planet 1, planet 2 and planet 3, respectively)

a)

b)

c)

d)

Fig. 11Motion trajectories: a) and b) motion trajectory of the carrier under torque load and non-torque load, respectively; c) and d) motion trajectory of the ring under torque load and non-torque load, respectively

a)

b)

c)

d)

When the non-torque load was considered, the same held true in that the motion trajectories of the carrier and ring, and their vibration equilibrium positions, were hardly affected by it, as shown in Fig. 11. In addition, their motion trajectories were also analogous to a circle. Given that the ring is a fixed part, its restraint stiffness is set to be relatively large. As a result, there was an obvious difference in the offsets of the vibration equilibrium positions of the carrier in the $X$ and $Y$ directions, as they were markedly larger than those of the ring. Compared with the torque load condition, the center of motion trajectory of the carrier changes from (0.009, 0.038) to (12.79, 20.52) when the non-torque load is considered. In other words, the offsets of their vibration equilibrium positions in the $X$ and $Y$ directions are increased by 1420 and 539 times, respectively, as shown in Fig. 11(a) and (b). However, compared with Fig. 10(c) and (d), the offsets of the vibration equilibrium positions of the ring in the $X$ and $Y$ directions are significantly reduced, and its center of motion trajectory goes from (–0.009, –0.008) to (–0.44, 0.36); increasing by 47.89 and 46 times.

The load-sharing coefficients of the sun-planet gear pairs and planet-ring gear pairs are shown in Fig. 12. S-P$i$ ($i=$ 1, 2, and 3) in Fig. 12(a) and (b) represent the external meshing pair formed by the sun and the $i$th planet, and P$i$-R ($i=$ 1, 2, and 3) in Fig. 12 (c) and (d) denote the internal meshing pair formed by the $i$th planet and ring. The vibration equilibrium positions of the sun, planets, ring, and carrier deviated significantly under the action of the radial loads $F_X$ and $F_Y$ resulting in a deterioration of the load-sharing characteristics of the sun-planet gear pairs and planet-ring gearpairs. The load-sharing coefficients of the sun-planet gear pair and the planet-ring gear pair increased from 1.087 to 1.237, and from 1.062 to 1.189, increased by 13.80 % and 11.96 %, respectively. Obviously, the influence of the non-torque load on the load-sharing characteristics of the sun-planet gear pairs is greater than that of the planet-ring gear pairs. This is mainly caused by the offset of the vibration equilibrium position of the sun being larger than that of the ring.

Fig. 12Load-sharing coefficient (LSC): a) and b) LSC of sun-planet gear pairs under torque load and non-torque load, respectively; c) and d) LSC of planet-ring gear pairs under torque load and non-torque load, respectively

a)

b)

c)

d)

The variations in the vibration equilibrium position of each central member are listed in Table 4. Obviously, in addition to the ring, the offsets of the other central members in the $Y$ direction are significantly greater than those in the $X$ direction, when the non-torque load is considered. This is because the radial load acting in the $X$ direction is significantly less than that in the $Y$ direction, as shown in Fig. 7. However, it is noteworthy to mention that the increased multiple of the offset of the vibration equilibrium position in the $X$ direction is obviously higher than that of the $Y$ direction. Moreover, the increased proportion of the offset of the vibration equilibrium position of the carrier is the largest, whether in the $X$ or $Y$ direction, closely followed by the sun and the wheel. This is because the carrier is closest to the point of the load application, followed by the sun and wheel. Similar variations can be seen in the LSC of each gear pair, as listed in Table 5. This indicates that the closer the central members are to the point of the load application, the more their vibration equilibrium positions and load-sharing characteristics are affected by it.

The offsets of the vibration equilibrium position of the central members at different non-torque loads are shown in Fig. 13. For convenience, the relationship between the resultant force $F_R$ and the offset, shown in the figure, where $F_R=\sqrt{F_X^2+F_Y^2}$. It is clear that the offsets of the vibration equilibrium position increase almost linearly with the increase in the resultant force $F_R$. Whether in the $X$ or $Y$ direction, the vibration equilibrium positions of the wheel, sun, and carrier are shifted in the positive direction of either the $X$-axis or $Y$-axis, with an increase in the resultant force $F_R$, as shown in Fig. 13(a), (b), and (c). However, with the increase in the resultant force $F_R$, the vibration equilibrium position of the ring in the $X$ direction is shifted in the negative direction of the $X$-axis, but the vibration equilibrium position of the ring in the $Y$ direction is shifted in the positive direction of the $Y$-axis, as shown in Fig. 13(d). As a result, the load-sharing characteristics of all gear pairs deteriorate further with the increase in the resultant force $F_R$, as shown in Fig. 14.

Fig. 13Offset of vibration equilibrium position at different non-torque loads: a) wheel; b) sun; c) carrier; and d) ring

a)

b)

c)

d)

Fig. 14Root mean square (RMS) of LSC at different non-torque loads: a) pinion-wheel gear pairs; b) sun-planet gear pairs; c) planet-ring gear pairs

a)

b)

c)

Fig. 15Offset of vibration equilibrium position at different speeds: a) wheel; b) sun; c) carrier; d) ring

a)

b)

c)

d)

Fig. 16RMS of LSC at different speeds: a) pinion-wheel gear pairs; b) sun-planet gear pairs; c) planet-ring gear pairs

a)

b)

c)

The offsets of the vibration equilibrium position of the central members at different speeds are shown in Fig. 15. When compared with Fig. 13, similarities can be observed. However, with the increase in the rotational speed, the offsets of the vibration equilibrium position of all the central members exhibit an approximate parabolic trend. When the rotational speed changes reaches the range of 800-1200 rpm, the offsets of the vibration equilibrium position of all central members increase slowly. The load-sharing characteristics of all gear pairs exhibit a very slow deterioration trend as well. Once the rotational speed exceeds 1200 rpm, the offsets of the vibration equilibrium position of the central members increase rapidly with the increase in speed, and the deterioration trend accelerates further, as shown in Fig. 16. In addition, the comparisons between Fig. 13 with Fig. 15, and Fig. 14 with Fig. 16, shows that the influence of non-torque load on the vibration equilibrium position and load-sharing characteristics is significantly greater than that of the rotational speed.

4. Conclusions

In this study, an electromechanical coupling dynamic model of the MPCTS, based upon the node finite element method, was established considering motor characteristics. The dynamic characteristics of the MPCTS are investigated under non-torque load conditions, and include a detailed study of the resulting motion trajectories and load-sharing characteristics. Based on the simulations and discussions, the following conclusions were drawn:

1) The motion trajectories of the pinions are linear. The angles and the lengths between the trajectories increase in different degrees when the non-torque load is considered. Here, the maximum angle deference and the maximum length deference increased from 0.0047° to 1.7474°, and from 20.33 μm to 408.56 μm, respectively. The motion trajectories of the wheel, sun, carrier, and ring are analogous to a circular shape. The motion trajectory of the planet is almost elliptical in shape, under the torque load condition, but is analogous to more rectangular under non-torque load condition. In addition, when the non-torque load is considered, the offsets of the vibration equilibrium position of the central members increase significantly, and the offsets in the $Y$ direction are noticeably greater than those in the $X$ direction.

2) For the central members, under the non-torque load condition, the increased multiples of their offsets of vibration equilibrium positions in the $X$ direction are clearly higher than those in the $Y$ direction. The closer the central members are to the applied load, the greater the increase in offset of the vibration equilibrium position of the carrier, closely followed by the sun and wheel. Similar characteristics were also observed in the load-sharing coefficient of each gear pair.

3) With an increase in the non-torque load, the offsets of the vibration equilibrium positions of the central members almost increased linearly, and the load-sharing characteristics of all gear pairs deteriorated further. However, with an increase in rotational speed, the offsets of the vibration equilibrium position of the central members exhibit a parabolic-like trend. When rotational speed is within the 800-1200 rpm range, the offsets increase slowly. However, once the rotational speed exceeds 1200 rpm, the offsets increase rapidly. Similarly, when the rotational speed is within the 800-1200 rpm range, the load-sharing characteristics of all gear pairs exhibit a very slow deterioration trend but accelerates further when the speed exceeds 1200 rpm. Moreover, the influence of the non-torque load on the vibration equilibrium position and load-sharing characteristics is significantly greater than that of increased rotational speed.