Dynamic Modeling and Vibration Characteristics Analysis of Transmission Process for Dual-Motor Coupling Drive System

The dual-motor coupling drive system (DCDS), which is widely used in electric vehicles, has attracted increasing attention due to its high transmission efficiency and economical fuel consumption. Current research has mainly focused on the control scheme of dual motors and has ignored the dynamic characteristics of the asymmetrical transmission structure. This paper presents an investigation of a dynamic model and an analysis method of the transmission process for the DCDS. The entire dynamic model of the DCDS was established by considering the planetary gear, differential bevel gear, and drive shaft with the transfer matrix method (TMM). Then, a detailed theoretical analysis was developed to study the influence of meshing stiffness and excitation source on the dynamic characteristics. Finally, the DCDS experimental platform was utilized to validate the effectiveness of the proposed dynamic model. For susceptibility to low-frequency vibrations, the first four natural frequencies and vibration modes of the DCDS were analyzed through the processing and analysis of acceleration signals. The experimental dynamic responses were generally consistent with the numerically computed results, which demonstrates the effectiveness of the proposed dynamic model with TMM. Furthermore, the proposed dynamic analysis method may be helpful for developing effective control algorithms to suppress vibrations and achieving highly smooth motions for electric vehicles.


Introduction
With the increasing demand for a green economy and environmental protection, modern vehicles are required to produce zero emissions, much less pollution, and little noise [1]. Electric vehicles, which use a drive motor and a battery instead of a traditional internal combustion engine, can achieve higher driving efficiency and control precision. Especially for the characteristics of constant torque at low speeds and constant power at high speeds, the power transmission system of electric vehicles can be tremendously simplified so as to improve the dynamic performance [2,3].
The dual-motor coupling drive system (DCDS) is a popular research topic in the field of electric vehicles [4]. Compared with a normal single-motor driving axle system, the DCDS, which adopts two motors rather than a high-peak-power motor, can significantly improve a vehicle's dynamic performance and energy efficiency [5]. The typical applications of the single-motor driving axle system and the DCDS are shown in Figures 1 and 2, respectively. In Figure 2, two brushless direct current motors (BDCM) are utilized with particular gear transmission equipment, such as a planetary gear, a differential bevel gear, and so on. The drive power of the dual-motor system is synthesized gear, a differential bevel gear, and so on. The drive power of the dual-motor system is synthesized and transferred by the gear transmission system. During the driving process, the main motor is operated at a rated speed to improve energy efficiency. The assistant motor is controlled to adjust the rotational speed. To improve the smoothness of stepless speed regulation, the dynamic and vibration characteristics of the DCDS need to be analyzed.  Implementation of a DCDS firstly requires high manufacturing and assembling precision to ensure motion transmission smoothness. Moreover, achieving dual-motor coordinate motion is restricted by the complicated power transmission chain from the dual motors to the wheels. The operation of a DCDS, which usually adopts a planetary gear and differential bevel gear set, can suffer from meshing deflection and torsional vibration. The actual driving dynamics are time varying and commonly affected by the continuous change of the meshing stiffness. Hence, it is crucial to study the dynamic response and natural characteristics of the transmission process for the DCDS.
At present, research on the DCDS has mainly focused on two aspects: optimizing the control mode and improving energy efficiency. The dynamic characteristics of the transmission process for the DCDS itself are ignored. Although many scholars have conducted extensive research on the dynamics of the gear set, research on the dynamic performance of the DCDS in electric vehicles has attracted less attention. A coupling structure diagram of a DCDS was studied in [6]. A dynamic model of the planetary gear and vehicle was presented by considering the transmission ratio, which is an oversimplification of the transmission system. Some studies have provided a seven-degrees-offreedom vehicle dynamic model. Several model-based vibration suppression algorithms have been proposed and tested by road experiments [7,8]. Considering external disturbances and resistance, Kim et al. [9] proposed a load torque estimation method for a parallel hybrid vehicle. The drive shafts from motor to wheel suffer from the transmission power, which is amplified by the reducing gearbox. The power transmission shaft also undertakes a load during the transmission process [10]. Xie et al. and Berriri et al. merged the rotational inertia of the engine and gearbox in an equivalent inertia [11,12]. The transmission system was simplified as a drive shaft model, and the total stiffness of the shaft was obtained by considering the axial stiffness and gearbox stiffness. Symmetry 2020, 12, x FOR PEER REVIEW 2 of 21 gear, a differential bevel gear, and so on. The drive power of the dual-motor system is synthesized and transferred by the gear transmission system. During the driving process, the main motor is operated at a rated speed to improve energy efficiency. The assistant motor is controlled to adjust the rotational speed. To improve the smoothness of stepless speed regulation, the dynamic and vibration characteristics of the DCDS need to be analyzed.  Implementation of a DCDS firstly requires high manufacturing and assembling precision to ensure motion transmission smoothness. Moreover, achieving dual-motor coordinate motion is restricted by the complicated power transmission chain from the dual motors to the wheels. The operation of a DCDS, which usually adopts a planetary gear and differential bevel gear set, can suffer from meshing deflection and torsional vibration. The actual driving dynamics are time varying and commonly affected by the continuous change of the meshing stiffness. Hence, it is crucial to study the dynamic response and natural characteristics of the transmission process for the DCDS.
At present, research on the DCDS has mainly focused on two aspects: optimizing the control mode and improving energy efficiency. The dynamic characteristics of the transmission process for the DCDS itself are ignored. Although many scholars have conducted extensive research on the dynamics of the gear set, research on the dynamic performance of the DCDS in electric vehicles has attracted less attention. A coupling structure diagram of a DCDS was studied in [6]. A dynamic model of the planetary gear and vehicle was presented by considering the transmission ratio, which is an oversimplification of the transmission system. Some studies have provided a seven-degrees-offreedom vehicle dynamic model. Several model-based vibration suppression algorithms have been proposed and tested by road experiments [7,8]. Considering external disturbances and resistance, Kim et al. [9] proposed a load torque estimation method for a parallel hybrid vehicle. The drive shafts from motor to wheel suffer from the transmission power, which is amplified by the reducing gearbox. The power transmission shaft also undertakes a load during the transmission process [10]. Xie et al. and Berriri et al. merged the rotational inertia of the engine and gearbox in an equivalent inertia [11,12]. The transmission system was simplified as a drive shaft model, and the total stiffness of the shaft was obtained by considering the axial stiffness and gearbox stiffness. Implementation of a DCDS firstly requires high manufacturing and assembling precision to ensure motion transmission smoothness. Moreover, achieving dual-motor coordinate motion is restricted by the complicated power transmission chain from the dual motors to the wheels. The operation of a DCDS, which usually adopts a planetary gear and differential bevel gear set, can suffer from meshing deflection and torsional vibration. The actual driving dynamics are time varying and commonly affected by the continuous change of the meshing stiffness. Hence, it is crucial to study the dynamic response and natural characteristics of the transmission process for the DCDS.
At present, research on the DCDS has mainly focused on two aspects: optimizing the control mode and improving energy efficiency. The dynamic characteristics of the transmission process for the DCDS itself are ignored. Although many scholars have conducted extensive research on the dynamics of the gear set, research on the dynamic performance of the DCDS in electric vehicles has attracted less attention. A coupling structure diagram of a DCDS was studied in [6]. A dynamic model of the planetary gear and vehicle was presented by considering the transmission ratio, which is an oversimplification of the transmission system. Some studies have provided a seven-degrees-of-freedom vehicle dynamic model. Several model-based vibration suppression algorithms have been proposed and tested by road experiments [7,8]. Considering external disturbances and resistance, Kim et al. [9] proposed a load torque estimation method for a parallel hybrid vehicle. The drive shafts from motor to wheel suffer from the transmission power, which is amplified by the reducing gearbox. The power transmission shaft also undertakes a load during the transmission process [10]. Xie et al. and Berriri et al. merged the rotational inertia of the engine and gearbox in an equivalent inertia [11,12]. The transmission system was simplified as a drive shaft model, and the total stiffness of the shaft was obtained by considering the axial stiffness and gearbox stiffness. The dynamic characteristics of the DCDS are critically influenced by the gear transmission system. The driving power is transferred from the dual motors to actuators by complicated planetary and differential gears. Previous research found that the meshing stiffness, which is varied periodically during the engagement process, can be the main excitation source. Wang and Sinha [13] established the lumped parameter model of the planetary gear to study its natural characteristics. Parker and colleagues used the finite element method to analyze the nonlinear phenomenon in the gear transmission and presented the phase regulation criterion to eliminate specific modes [14,15]. Several scholars studied the influence of errors and stiffness on gear contact load characteristics [16,17]. The dynamic model and response of the rotor system, which contains the planetary gear, was provided and the system stability was analyzed in [18,19]. Moreover, the dynamics of the dual-motor motion mechanism with constraints can be further applied in the asymptotical stability control [20]. Hence, how the meshing stiffness and structural features affect the natural characteristics and dynamic response of the DCDS needs to be further analyzed.
This paper is organized as follows: Section 2 presents the dynamic modeling and analysis of a DCDS with the transfer matrix method (TMM). The dynamic model of the DCDS is established by considering the transmission process of the planetary gear, differential bevel gear, and drive shafts. The effects of meshing stiffness and structural features on the time-frequency dynamic response of the DCDS are analyzed in Section 3. Reports on the dynamic characteristics test experiment that was carried out, which adopted the DCDS platform and an NI vibration analyzer. The natural characteristics and dynamic response were detected and analyzed by discrete Fourier transform (DFT). Further, the effectiveness of the proposed dynamic model for the DCDS was verified by comparing the results of numerical computing to the experimental test. Conclusions are drawn in Section 4.

Dynamic Model of DCDS with TMM
In this section, the dynamic model of the DCDS by considering the asymmetrical structure of the DCDS is proposed with TMM. Then, the generating principles, geometries, and dynamics of the DCDS are investigated. Based on this, the dynamic analysis is then carried out in Section 3.

Establishment of the Model for the Whole DCDS
The proposed DCDS consists of dual motors, a planetary gear, a differential gear, and drive shafts, as shown in Figure 3. To provide high motion stability and accuracy, the driving process is accomplished by two-set power sources. The main motor actuator which takes charge of the rotation of sun gear Z s is mounted on the frame. The assistant motor actuator which is connected to the internal gear ring Z r is used to adjust the motion velocity. Due to the differential speed of the dual motors and the gear engagement law, both the revolution and rotation of planet gear Z p can be generated. The planet carrier Z c which supports the planets is rotated around the center line. Furthermore, the differential bevel gear is promoted by the assisted gear Z 3 , which is engaged with the planet carrier. During the driving process, the sun gear Z s is promoted by motor 1 at a rated speed to improve the energy efficiency. The gear ring Z r is coordinately rotated by assistant motor 2 to adjust the rotational speed. The electric vehicle can be driven by the differential, which is composed of two pairs of bevel gears.  The generalized dynamic model of the DCDS is presented in Figure 4. To study the dynamic characteristics of the transmission process, the mass (inertia) of the dual motors and the main body of the gears are supposed to be a lumped mass (inertia). The shaft segment between each lump mass is considered to be a massless torsional spring. Considering the torsional vibrations of the DCDS, the state vector of the lumped parameter shaft contains five shaft segments in the system, as shown in Figure 3. Shafts I and II are the input shafts, which are connected to the dual motors; Shafts III and IV are the transmission shaft; and Shaft V is the output shaft, which drives the wheels directly. The overall state vector of the DCDS can be defined as The generalized dynamic model of the DCDS is presented in Figure 4. To study the dynamic characteristics of the transmission process, the mass (inertia) of the dual motors and the main body of the gears are supposed to be a lumped mass (inertia). The shaft segment between each lump mass is considered to be a massless torsional spring. The generalized dynamic model of the DCDS is presented in Figure 4. To study the dynamic characteristics of the transmission process, the mass (inertia) of the dual motors and the main body of the gears are supposed to be a lumped mass (inertia). The shaft segment between each lump mass is considered to be a massless torsional spring. Considering the torsional vibrations of the DCDS, the state vector of the lumped parameter shaft contains five shaft segments in the system, as shown in Figure 3. Shafts I and II are the input shafts, which are connected to the dual motors; Shafts III and IV are the transmission shaft; and Shaft V is the output shaft, which drives the wheels directly. The overall state vector of the DCDS can be defined as Considering the torsional vibrations of the DCDS, the state vector of the lumped parameter shaft contains five shaft segments in the system, as shown in Figure 3. Shafts I and II are the input shafts, Symmetry 2020, 12, 1171 5 of 20 which are connected to the dual motors; Shafts III and IV are the transmission shaft; and Shaft V is the output shaft, which drives the wheels directly. The overall state vector of the DCDS can be defined as where Z i (i = I, II, III, IV, and V) are the state vectors of each shaft, and the state equation is presented in the following formula: where Z jL and Z jD are the input drive torque and output load torque of the jth substructure, and j = 1, 2, . . . , 9. T j is the transmission matrix between the two adjacent substructures. Then, the whole transmission matrix of the DCDS based on TMM can be expressed as The following principles need to be followed when the whole transfer matrix T W is constructed, owing to the difference of the location for each substructure and the coupling relationship between two transmission components. The states of each substructure are transmitted simultaneously with the overall state vector Z.
(1) If the state of the substructure is transmitted independently in a certain period, the transfer matrix is established by the pure rotational matrix according to the corresponding position in the whole transfer matrix; in addition, the coupling transmission matrix is 0.
(2) If the state of the substructure is not transmitted, the relative positions can be expressed by the unit matrix E (the same order as the transfer matrix of individual pair of gears).
(3) If the state of the substructure is transmitted with both coupling and rotation motions, the coupling unit of the whole transfer matrix is not 0. The transfer matrix needs to be established by the coordination relation of specific torque balance and meshing deflection.
Hence, the diagram of the substructure division by TMM is shown in Figure 4, and the whole transfer matrix of the DSDC can be established as where E 1 , E 2 , E 3 , and E 4 are 5 order, 10 order, 15 order, and 20 order unit matrices, respectively; T I , T II , T III , T IV , and T V represent the transfer matrix of Shafts I-V at certain transmission processes; and T C n 11 , T C n 12 , T C n 21 , and T C n 22 are the elements of the coupling transfer matrix for the engagement of the individual gear and planetary gears.

Modeling of Meshing Process for an Individual Pair of Gears
The meshed model of an individual pair of gears is shown in Figure 5. The meshing dynamic model of can be expressed as where θ 1 and θ 2 are the rotational angles of the driving gear and the driven gear, respectively. θ D and θ L are the rotational angles caused by input drive torque and output load torque, respectively. Similarly, T 1D and T 2D are the input drive torques of the driving gear and the driven gear, respectively. T 1L and T 2L are the output load torques of the driving gear and the driven gear, respectively. R b1 and R b2 are the base radii. e(t) is the geometric transmission error along with the meshing line. k m is the meshing stiffness, and c is the damp.
Symmetry 2020, 12, x FOR PEER REVIEW 6 of 21 where 1  and 2  are the rotational angles of the driving gear and the driven gear, respectively. D  and L  are the rotational angles caused by input drive torque and output load torque, respectively.
Similarly, T1D and T2D are the input drive torques of the driving gear and the driven gear, respectively. T1L and T2L are the output load torques of the driving gear and the driven gear, respectively. Rb1 and Rb2 are the base radii. e(t) is the geometric transmission error along with the meshing line. km is the meshing stiffness, and c is the damp. Meshing error is caused by manufacturing and assembly errors. The basic frequency of the meshing error is equal to the engagement frequency, which has obvious periodicity. Then, the meshing error can be obtained by Fourier transform: where m e is the average error.
Then, the meshing model of an individual pair of gears can be expressed as The coefficient matrix in the right side of (7) is the transfer matrix, which can be rewritten as Meshing error is caused by manufacturing and assembly errors. The basic frequency of the meshing error is equal to the engagement frequency, which has obvious periodicity. Then, the meshing error can be obtained by Fourier transform: where e m is the average error.
Then, the meshing model of an individual pair of gears can be expressed as Symmetry 2020, 12, 1171 7 of 20 The coefficient matrix in the right side of (7) is the transfer matrix, which can be rewritten as The second matrix on the right side of (8) denotes the pure rotational motion of each gear. If the gear is considered to be a rigid body, then the meshing deflection of the gear is zero. In this case, the transfer process can be expressed as a pure lumped inertia matrix. The first matrix on the right side of (8) points out the mechanical coupling between the pair of gears, which is affected by the meshing stiffness k m , damp c, and base radius.

Modeling of the Transmission Process for the Planetary Gear
Since the planetary gear set has the significant feature of mass concentration, the lumped parameter model is proposed in this subsection. The lumped parameter modeling method contains two aspects: a pure torsional model and a translational-rotational model. In the DCDS of an electric vehicle, since the ratio of the support stiffness of the gears and the meshing stiffness at its engagement surface is greater than 10, the natural characteristics computed by these two models are quite close [21]. Assume that the three planet gears are uniformly distributed along the circumference and the support stiffness and inertia of each planet gear are equal. The torsional model of the transmission process for the planetary gear is shown in Figure 6.
The second matrix on the right side of (8) denotes the pure rotational motion of each gear. If the gear is considered to be a rigid body, then the meshing deflection of the gear is zero. In this case, the transfer process can be expressed as a pure lumped inertia matrix. The first matrix on the right side of (8) points out the mechanical coupling between the pair of gears, which is affected by the meshing stiffness m k , damp c, and base radius.

Modeling of the Transmission Process for the Planetary Gear
Since the planetary gear set has the significant feature of mass concentration, the lumped parameter model is proposed in this subsection. The lumped parameter modeling method contains two aspects: a pure torsional model and a translational-rotational model. In the DCDS of an electric vehicle, since the ratio of the support stiffness of the gears and the meshing stiffness at its engagement surface is greater than 10, the natural characteristics computed by these two models are quite close [21]. Assume that the three planet gears are uniformly distributed along the circumference and the support stiffness and inertia of each planet gear are equal. The torsional model of the transmission process for the planetary gear is shown in Figure 6.   The sun gear is designated as s, the planet gear as p i (i = 1, 2, 3), the carrier as c, and the gear ring as r. The rotational angles of the sun gear, planet gear, and carrier are labeled as θ s , θ pi , and θ c . The meshing stiffness between the sun gear and the ith planet gear is k spi . The base radii of the sun gear, planet gear, and gear ring are R bs , R bp , and R br . Several sets of coordinates can be used to describe the geometrical relationship of the DCDS. One set is given by the absolute coordinate, which is established based on the fixed center O(x, y, z). The second set is given by the actual coordinate of the sun gear, which is designated as {S s -X s Y s Z s }. The planet gears are given by the equivalent coordinates {S p -X p Y p Z p }, which orbit around the sun gear with a variable angle velocity w z .
Suppose that, under the input drive torque, the rotational direction of the sun gear and carrier is positive, and the rotational direction of the planet gear is the opposite. The input drive torque and output torque are T D and T L . Then, the vibration differential equation can be established as Considering the transmission balance between the input drive torque and output load torque, the meshing model of the sun gear can be expressed as Similarly, the meshing model of the internal gear ring can be expressed as The meshing model of the planet gear can be expressed as The meshing model of the carrier can be expressed as where I s = m s R 2 bs , I pi = m pi R 2 bpi , and I c = m c R 2 bc are the equivalent masses of the sun gear, planet gear, and carrier. F D = T D /R bs and F L = T L /R bs are the equivalent forces.
During the engagement process, the meshing force is determined by the relative displacement along with the meshing line and the meshing stiffness of each component. The relative displacement along with the meshing line is composed of the rotational displacement and the eccentric error of each gear. The eccentric errors which contain the manufacturing and vibration errors of the sun gear, planet gear, and gear ring are labeled E s , E pi , and E r . ϕ s is the angle between E s and the external meshing line of the ith planet gear. δ s is the angle between E s and the x axis. γ pi and ϕ pi are the angles between E pi and the internal and external meshing lines, respectively. δ pi is the angle between E pi and the x pi axis. ϕ r is the angle between E r and the internal meshing line of the ith planet gear. δ r is the angle between E r and the x axis. From Figures 7 and 8, the following geometric relationship can be obtained: where w s , w p , and w c are the rotational speeds of the sun gear, planet gear, and gear ring. ψ i is the initial angle between the ith planet gear and the x axis. The projection of E s and E pi on the external meshing line can be expressed as The projection of E pi and E r on the internal meshing line can be expressed as The comprehensive error along with the meshing line can be obtained by combining (6), (15), and (16): For the ith planet gear, the relative displacements along with the meshing line can be deduced as where x s = R bs θ s and x pi = R bpi θ bpi are the rotational displacements of the sun gear and planet gear, respectively. Then, the drive force and damp of the ith gear can be expressed as where c pi = 2ξ k pi 1/m s +1/m pi is the coefficient of the meshing damp. Hence, the dynamic model of the planetary gear shown in (10)- (13) can be rewritten as follows: The coefficient matrix on the right side of (21) is the transfer matrix of the planetary gear, which can be defined as   The sun gear is designated as s, the planet gear as pi (i = 1, 2, 3), the carrier as c, and the gear ring as r. The rotational angles of the sun gear, planet gear, and carrier are labeled as s  , pi  , and c  .
The meshing stiffness between the sun gear and the ith planet gear is kspi. The base radii of the sun gear, planet gear, and gear ring are Rbs, Rbp, and Rbr. Several sets of coordinates can be used to describe the geometrical relationship of the DCDS. One set is given by the absolute coordinate, which is established based on the fixed center   , , O x y z . The second set is given by the actual coordinate of the sun gear, which is designated as {Ss-XsYsZs}. The planet gears are given by the equivalent coordinates {Sp-XpYpZp}, which orbit around the sun gear with a variable angle velocity wz.
Suppose that, under the input drive torque, the rotational direction of the sun gear and carrier is positive, and the rotational direction of the planet gear is the opposite. The input drive torque and output torque are TD and TL. Then, the vibration differential equation can be established as   The sun gear is designated as s, the planet gear as pi (i = 1, 2, 3), the carrier as c, and the gear ring as r. The rotational angles of the sun gear, planet gear, and carrier are labeled as s  , pi  , and c  .
The meshing stiffness between the sun gear and the ith planet gear is kspi. The base radii of the sun gear, planet gear, and gear ring are Rbs, Rbp, and Rbr. Several sets of coordinates can be used to describe the geometrical relationship of the DCDS. One set is given by the absolute coordinate, which is established based on the fixed center   , , O x y z . The second set is given by the actual coordinate of the sun gear, which is designated as {Ss-XsYsZs}. The planet gears are given by the equivalent coordinates {Sp-XpYpZp}, which orbit around the sun gear with a variable angle velocity wz.
Suppose that, under the input drive torque, the rotational direction of the sun gear and carrier is positive, and the rotational direction of the planet gear is the opposite. The input drive torque and output torque are TD and TL. Then, the vibration differential equation can be established as

Modeling of the Transmission Process for the Differential Bevel Gear
The differential bevel gears can be considered to be composed of one driving gear and two driven gears. The meshing dynamic model of the differential bevel gears can be expressed as where θ cL , θ 1L , and θ 2L are the output load rotational angle of the driving gear and driven gear, and T cL , T c1 , and T c2 are the output torque of the driving gear and driven gear, respectively. The first matrix on the right side of (28) is the transfer matrix of the differential bevel gear without the rotational motion of lumped inertia.
The shaft segments I-V are considered to be a massless torsional spring with a uniform cross section. Considering the effect of torsional deformation, the transfer matrix of each shaft segment can be obtained as where G is the shear modulus, L i and I si are the length and rotational inertia of the ith shaft segment, and k s is the torsional stiffness.

Dynamic Characteristics Analysis and Test of the Transmission Process for DCDS
The transfer matrix can be obtained by combining (8), (22), (28), and (29) sequentially from substructure 1 to substructure 9. The dynamic response can be computed by the characteristics equations of w. After the natural frequency computing, the certain order of frequency can be substituted in the whole transfer matrix. Then, one can obtain the linear state equation and the state variables. The eigenvalues represent the natural frequency, and the state vectors are the vibration mode shapes.

Dynamic Characteristics Analysis
Solve the characteristic equation w i (i = 1 · · · 5), which shows that the natural frequencies can be obtained. The rotational speeds of the dual motors were set to 0-3000 r/min. The diameters of shaft segments I-V were 0.020, 0.018, 0.025, 0.030, and 0.030 m, respectively, and the modulus of elasticity E = 2.06 × 10 11 Pa. The parameters of the planetary gear and differential bevel gear for the computational analysis are shown in Table 1. Since the main exciting sources of the DCDS are the alternating resistance torques, the torsional and radial vibrations were the main modes which were modeled in Section 2, and the axial motion was ignored. To study the effect of axial motion on the dynamic characteristics, the dynamic analysis of the planetary gear, which is the most complex component in the DCDS, is presented. In practice, the driving force which is transferred by the contact of the teeth's surface is composed of elastic force and damping force. In the numerical computation and simulation, the impact function method (IFM) was adopted to express the contact force: where K is the stiffness coefficient, q 0 is the initial distance of the contact surface, q is the actual motion displacement, and STEP is the step function. The simulated results are shown in Figures 9-11. As the rotational speed increased from 1 to 50 r/s, the torsional, radial, and axial vibrations became more severe. The torsional and radial vibrations were almost the same when considering axial motion or not, which means that the axial motion made no difference to the dynamic characteristics. Moreover, Figure 11 shows that the difference of the axial vibration amplitudes in the z axis was several times lower than the torsional and radial vibrations. Each pair of meshing gears had a symmetrical structure, and the dynamic characteristics had periodic behavior. Figure 10 shows that the radial vibrations were very close in the x and y axes. Hence, the axial vibration of the gears could be neglected in the dynamic analysis of the DCDS, especially at low frequencies.
where K is the stiffness coefficient, q0 is the initial distance of the contact surface, q is the actual motion displacement, and STEP is the step function. The simulated results are shown in Figures 9-11. As the rotational speed increased from 1 to 50 r/s, the torsional, radial, and axial vibrations became more severe. The torsional and radial vibrations were almost the same when considering axial motion or not, which means that the axial motion made no difference to the dynamic characteristics. Moreover, Figure 11 shows that the difference of the axial vibration amplitudes in the z axis was several times lower than the torsional and radial vibrations. Each pair of meshing gears had a symmetrical structure, and the dynamic characteristics had periodic behavior. Figure 10 shows that the radial vibrations were very close in the x and y axes. Hence, the axial vibration of the gears could be neglected in the dynamic analysis of the DCDS, especially at low frequencies.   Then, the time-frequency domain dynamic characteristics of the different substructures in the transmission process were computed and analyzed by the dynamic model. Figure 12 shows the dynamic response of the DCDS with the driving of the single main motor 1 at the speed of 600 r/min. At this time, the rotational speeds of the gears related to the assistant motor 2 were 0, and the ring gear could be considered to be fixed. Figure 13 shows the dynamic response of the DCDS with the driving of the single assistant motor 2 at the same rotational speed. At this time, the rotational speeds of the sun gear and related gears were 0. From the time domain curve, the varying trend of each substructure driven by the main motor 1 and assistant motor 2 was consistent. For the asymmetrical transmission structure, the amplitude of the dynamic response driven by the main motor 1 was relatively larger than that driven by the assistant motor 2. As can be seen in the frequency domain  Then, the time-frequency domain dynamic characteristics of the different substructures in the transmission process were computed and analyzed by the dynamic model. Figure 12 shows the dynamic response of the DCDS with the driving of the single main motor 1 at the speed of 600 r/min. At this time, the rotational speeds of the gears related to the assistant motor 2 were 0, and the ring gear could be considered to be fixed. Figure 13 shows the dynamic response of the DCDS with the driving of the single assistant motor 2 at the same rotational speed. At this time, the rotational speeds of the sun gear and related gears were 0. From the time domain curve, the varying trend of each substructure driven by the main motor 1 and assistant motor 2 was consistent. For the asymmetrical transmission structure, the amplitude of the dynamic response driven by the main motor 1 was relatively larger than that driven by the assistant motor 2. As can be seen in the frequency domain curve, in the initial state, the dynamic response of each substructure mainly consisted of first, second, Then, the time-frequency domain dynamic characteristics of the different substructures in the transmission process were computed and analyzed by the dynamic model. Figure 12 shows the dynamic response of the DCDS with the driving of the single main motor 1 at the speed of 600 r/min. At this time, the rotational speeds of the gears related to the assistant motor 2 were 0, and the ring gear could be considered to be fixed. Figure 13 shows the dynamic response of the DCDS with the driving of the single assistant motor 2 at the same rotational speed. At this time, the rotational speeds of the sun gear and related gears were 0. From the time domain curve, the varying trend of each substructure driven by the main motor 1 and assistant motor 2 was consistent. For the asymmetrical transmission structure, the amplitude of the dynamic response driven by the main motor 1 was relatively larger than that driven by the assistant motor 2. As can be seen in the frequency domain curve, in the initial state, the dynamic response of each substructure mainly consisted of first, second, and third harmonic frequencies. The amplitude of each harmonic frequency increased with the increasing of the rotational speed, and the stabilizing time became longer as well.  The natural characteristics are related to the meshing stiffness of the gear. The variation of the natural frequency along with the increase of the meshing stiffness is shown in Figure 14. The first and second vibration modes were barely affected by the increasing of the meshing stiffness; the first two vibration modes were mainly caused by the mass wheel at the end of the shaft. The variation trend of the third and fourth order vibrations indicates that improving meshing stiffness can increase the natural frequency. The third vibration mode increased as the meshing stiffness improved from 0 to 10 N/m; the third order vibration was mainly generated by the bevel gear set, and the amplitude was obviously affected by the meshing stiffness. Meanwhile, the fourth vibration mode increased from 1167.1 to 1574.5 rad/s. The fourth order vibration mode reached a stable value more slowly than the first three vibration modes. When the meshing stiffness was set to 8 6.5 10  N/m, which is approximate to the actual value, the first four natural frequencies could be calculated as 219.604, 500.219, 1083.379, and 1389.928 rad/s, respectively.  The natural characteristics are related to the meshing stiffness of the gear. The variation of the natural frequency along with the increase of the meshing stiffness is shown in Figure 14. The first and second vibration modes were barely affected by the increasing of the meshing stiffness; the first two vibration modes were mainly caused by the mass wheel at the end of the shaft. The variation trend of the third and fourth order vibrations indicates that improving meshing stiffness can increase the natural frequency. The third vibration mode increased as the meshing stiffness improved from 0 to 10 N/m; the third order vibration was mainly generated by the bevel gear set, and the amplitude was obviously affected by the meshing stiffness. Meanwhile, the fourth vibration mode increased from 1167.1 to 1574.5 rad/s. The fourth order vibration mode reached a stable value more slowly than the first three vibration modes. When the meshing stiffness was set to 8 6.5 10  N/m, which is approximate to the actual value, the first four natural frequencies could be calculated as 219.604, 500.219, 1083.379, and 1389.928 rad/s, respectively. The natural characteristics are related to the meshing stiffness of the gear. The variation of the natural frequency along with the increase of the meshing stiffness is shown in Figure 14. The first and second vibration modes were barely affected by the increasing of the meshing stiffness; the first two vibration modes were mainly caused by the mass wheel at the end of the shaft. The variation trend of the third and fourth order vibrations indicates that improving meshing stiffness can increase the natural frequency. The third vibration mode increased as the meshing stiffness improved from 0 to 10 N/m; the third order vibration was mainly generated by the bevel gear set, and the amplitude was obviously affected by the meshing stiffness. Meanwhile, the fourth vibration mode increased from 1167.1 to 1574.5 rad/s. The fourth order vibration mode reached a stable value more slowly than the first three vibration modes. When the meshing stiffness was set to 6.5 × 10 8 N/m, which is approximate to the actual value, the first four natural frequencies could be calculated as 219.604, 500.219, 1083.379, and 1389.928 rad/s, respectively.
Symmetry 2020, 12, x FOR PEER REVIEW 16 of 21 Figure 14. The variation of natural frequency along with the increase of meshing stiffness.
The DCDS can be divided into an input shaft (Shafts I and II), an intermediate shaft (Shaft III), and an output shaft (Shafts IV and V). The normalized vibration modal shapes of the first four natural frequencies from the dual motors to the output shaft are shown in Figure 15. The meshing stiffness were set to 2, 5, and 8 MN/m. The relative deflection of each part could be obtained with Gauss elimination and the recurrence formula of the differential motion equation on the basis of the calculated natural frequency. Figure 15a,b indicate that the first two vibration modes were mainly caused by the dual motors, which were mounted at the input end of the equivalent shafting. Figure  15a shows that the increase of the meshing stiffness could obviously reduce the vibration amplitude of the front end of the DCDS. Figure 15b,d indicate that the change of the meshing stiffness had a distinct effect on the second and fourth vibration modes of the DCDS. Figure 15c shows that the third vibration mode shape with 8 MN/m stiffness was the opposite of the other two stiffness situations, while the variation rules were similar. With the increase of the meshing stiffness, the amplitudes of the vibration mode at the front end of the DCDS became lower. Overall, the variation trend of the third and fourth vibration modes had an increasing trend with the improving of the stiffness. Hence, increasing the meshing stiffness will suppress the vibration amplitude of the low-order vibration, especially the front and tail ends of the DCDS. Nevertheless, it may amplify the vibration of the highorder modes of the DCDS. Figure 15. The normalized model of the first four vibration modes with different meshing stiffness: (a) The first mode shape varying from the axial length; (b) the second mode shape varying from the axial length; (c) the third mode shape varying from the axial length; (d) the fourth mode shape varying from the axial length. The DCDS can be divided into an input shaft (Shafts I and II), an intermediate shaft (Shaft III), and an output shaft (Shafts IV and V). The normalized vibration modal shapes of the first four natural frequencies from the dual motors to the output shaft are shown in Figure 15. The meshing stiffness were set to 2, 5, and 8 MN/m. The relative deflection of each part could be obtained with Gauss elimination and the recurrence formula of the differential motion equation on the basis of the calculated natural frequency. Figure 15a,b indicate that the first two vibration modes were mainly caused by the dual motors, which were mounted at the input end of the equivalent shafting. Figure 15a shows that the increase of the meshing stiffness could obviously reduce the vibration amplitude of the front end of the DCDS. Figure 15b,d indicate that the change of the meshing stiffness had a distinct effect on the second and fourth vibration modes of the DCDS. Figure 15c shows that the third vibration mode shape with 8 MN/m stiffness was the opposite of the other two stiffness situations, while the variation rules were similar. With the increase of the meshing stiffness, the amplitudes of the vibration mode at the front end of the DCDS became lower. Overall, the variation trend of the third and fourth vibration modes had an increasing trend with the improving of the stiffness. Hence, increasing the meshing stiffness will suppress the vibration amplitude of the low-order vibration, especially the front and tail ends of the DCDS. Nevertheless, it may amplify the vibration of the high-order modes of the DCDS.

Experiment Setup and Test
To verify the proposed comprehensive dynamic model of the DCDS, the natural frequency and dynamic response were tested by experiment. The experimental platform of the DCDS is shown in Figure 16. The variation of the acceleration vibration response of the DCDS was measured by an NI vibration analyzer. Triaxial acceleration sensors were selected to collect vibration signals, and the dynamic performance of the DCDS was tested by eddy current sensors and an NI vibration monitoring analyzer, as shown in Figure 17, so that the torsional and radial vibration of the gear set and power shaft could be measured.
while the variation rules were similar. With the increase of the meshing stiffness, the amplitudes of the vibration mode at the front end of the DCDS became lower. Overall, the variation trend of the third and fourth vibration modes had an increasing trend with the improving of the stiffness. Hence, increasing the meshing stiffness will suppress the vibration amplitude of the low-order vibration, especially the front and tail ends of the DCDS. Nevertheless, it may amplify the vibration of the highorder modes of the DCDS.

Experiment Setup and Test
To verify the proposed comprehensive dynamic model of the DCDS, the natural frequency and dynamic response were tested by experiment. The experimental platform of the DCDS is shown in Figure 16. The variation of the acceleration vibration response of the DCDS was measured by an NI vibration analyzer. Triaxial acceleration sensors were selected to collect vibration signals, and the dynamic performance of the DCDS was tested by eddy current sensors and an NI vibration monitoring analyzer, as shown in Figure 17, so that the torsional and radial vibration of the gear set and power shaft could be measured.   Figure 17 shows the filtered average curve after the separation of the trend term and the noise signal by an IIR digital filter. In order to obtain displacement signals, the acceleration signal needed to be integrated twice through the Finally, the frequency characteristic curve

Experiment Setup and Test
To verify the proposed comprehensive dynamic model of the DCDS, the natural frequency and dynamic response were tested by experiment. The experimental platform of the DCDS is shown in Figure 16. The variation of the acceleration vibration response of the DCDS was measured by an NI vibration analyzer. Triaxial acceleration sensors were selected to collect vibration signals, and the dynamic performance of the DCDS was tested by eddy current sensors and an NI vibration monitoring analyzer, as shown in Figure 17, so that the torsional and radial vibration of the gear set and power shaft could be measured.   Figure 17 shows the filtered average curve after the separation of the trend term and the noise signal by an IIR digital filter. In order to obtain displacement signals, the acceleration signal needed to be integrated twice through the Finally, the frequency characteristic curve The rotational speed of the dual motors could be set from 500 to 3000 r/min. According to sampling theorem f s > 2 f m , the sampling frequency was 1000 Hz in order to ensure the reliability of the experimental data. The vibration voltage signal a(k) (k = 0, 1, 2, . . . , N, where N is the number of samples) collected by the acceleration sensor was converted into an acceleration signal with unit conversion; then, the time domain signal could be obtained. Figure 17 shows the filtered average curve after the separation of the trend term and the noise signal by an IIR digital filter. In order to obtain displacement signals, the acceleration signal needed to be integrated twice through the trapezoidal integral formula y(k) = ∆t k i=1 a (i−1) + a (i) /2. Finally, the frequency characteristic curve could be obtained by DFT: Y(k) = N−1 r=0 y(r)e − j2πkr/N . Figure 18 shows the original time domain characteristic curve and the first four natural frequencies by DFT of the DCDS. The comparison of the first four natural frequencies between the dynamic model and the experimental test is shown in Table 2. On the whole, the first four natural frequencies of the experimental test corresponded closely to the computed result of the dynamic model. The relative errors, which ranged from 2.61% to 7.75%, were within allowable error. Due to the presence of damping, the actual experimental results were relatively higher than the computed frequencies. The comparison results confirmed that the computed results were approximate to the experimental natural frequencies, which illustrates the efficiency of the proposed dynamic model. The experimental results of the dynamic response for the DCDS are shown in Figure 19. The variation trend of the experimental dynamic responses was consistent with that of the computed results. Furthermore, the first four harmonic frequencies were also approximate, which illustrates that the dynamic model can reflect the actual operation process of the DCDS. Figure 19a-d show that the first and second harmonic frequencies were significantly stimulated in substructures 2, 4, 6, and 8. However, the amplitudes of the second harmonic frequency of substructures 6 and 8 were much lower than substructures 2 and 4, which indicates that the second harmonic component was mainly caused by Shafts I and II connected with the dual motors. The amplitude of the second harmonic frequency became lower as the distance away from the dual motors increased. Moreover, the third The comparison of the first four natural frequencies between the dynamic model and the experimental test is shown in Table 2. On the whole, the first four natural frequencies of the experimental test corresponded closely to the computed result of the dynamic model. The relative errors, which ranged from 2.61% to 7.75%, were within allowable error. Due to the presence of damping, the actual experimental results were relatively higher than the computed frequencies. The comparison results confirmed that the computed results were approximate to the experimental natural frequencies, which illustrates the efficiency of the proposed dynamic model. The experimental results of the dynamic response for the DCDS are shown in Figure 19. The variation trend of the experimental dynamic responses was consistent with that of the computed results. Furthermore, the first four harmonic frequencies were also approximate, which illustrates that the dynamic model can reflect the actual operation process of the DCDS. Figure 19a-d show that the first and second harmonic frequencies were significantly stimulated in substructures 2, 4, 6, and 8.
However, the amplitudes of the second harmonic frequency of substructures 6 and 8 were much lower than substructures 2 and 4, which indicates that the second harmonic component was mainly caused by Shafts I and II connected with the dual motors. The amplitude of the second harmonic frequency became lower as the distance away from the dual motors increased. Moreover, the third and fourth harmonic frequencies were stimulated, which was unapparent in the computed results. That was perhaps caused by nonlinear factors such as meshing backlash and friction, which can be a future research topic.

Conclusions
A theoretical and experimental study of the dynamic characteristics of the DCDS by TMM was proposed. The main conclusions of this study are as follows: (1) Firstly, the transmission dynamics model for the DCDS was developed, and the whole system was divided into several substructures by considering the asymmetrical transmission structure of the dual motors. The dynamic model of the transmission process for the planetary gear and differential bevel gear was proposed by considering the torsional vibration, transmission error, and meshing stiffness, which can be excited in the practical driving process. The dynamic characteristics and vibration shapes were analyzed with different driving mode and time-varying stiffness.
(2) Then, a detailed analysis of the variation trend of the natural characteristics and dynamic responses was numerically computed and presented. The effect of the meshing force was affected by the rotational speed of the dual motors. Compared with the torsional vibrations, the axial vibration can be neglected in the dynamic analysis of a DCDS, especially in low frequencies.
(3) The natural characteristics were analyzed by considering the time-varying meshing stiffness of the gear. The first two vibration modes are mainly caused by the input of the equivalent shafting, and increasing the meshing stiffness can obviously reduce the vibration amplitude. The main excitation sources are the input front end and output tail end of the transmission chain. Increasing the meshing stiffness will suppress the vibration amplitude of the low-order vibration. Nevertheless, it may amplify the vibration of the high-order modes of a DCDS.

Conclusions
A theoretical and experimental study of the dynamic characteristics of the DCDS by TMM was proposed. The main conclusions of this study are as follows: (1) Firstly, the transmission dynamics model for the DCDS was developed, and the whole system was divided into several substructures by considering the asymmetrical transmission structure of the dual motors. The dynamic model of the transmission process for the planetary gear and differential bevel gear was proposed by considering the torsional vibration, transmission error, and meshing stiffness, which can be excited in the practical driving process. The dynamic characteristics and vibration shapes were analyzed with different driving mode and time-varying stiffness.
(2) Then, a detailed analysis of the variation trend of the natural characteristics and dynamic responses was numerically computed and presented. The effect of the meshing force was affected by the rotational speed of the dual motors. Compared with the torsional vibrations, the axial vibration can be neglected in the dynamic analysis of a DCDS, especially in low frequencies.
(3) The natural characteristics were analyzed by considering the time-varying meshing stiffness of the gear. The first two vibration modes are mainly caused by the input of the equivalent shafting, and increasing the meshing stiffness can obviously reduce the vibration amplitude. The main excitation sources are the input front end and output tail end of the transmission chain. Increasing the meshing stiffness will suppress the vibration amplitude of the low-order vibration. Nevertheless, it may amplify the vibration of the high-order modes of a DCDS.
(4) Finally, the dynamic characteristics experiment for the DCDS was carried out. The natural characteristics and dynamic responses were detected and analyzed by a vibration monitoring analyzer. The first four natural frequencies and variation trend of the experimental test were generally consistent with the numerically computed results. The relative errors ranged from 2.61% to 7.75%, which illustrates the efficiency of the proposed dynamic model. Furthermore, the experimental results of the dynamic responses indicate that the first and second harmonics of meshing frequencies are mainly caused by the shafts connected to the dual motors, which is consistent with the analysis of the dynamic model. The third and fourth harmonics are mainly caused by the planetary gear and bevel gear, respectively. In the practical driving process, the damping and control scheme can be specifically designed.
The main implication of this study is the dynamic characteristics analysis of the transmission process for the DCDS in electric vehicles. The modeling and experimental test method, which has been validated by the vibration detection system, can be used in the dynamic analysis of a DCDS. The proposed analysis and test method can also help develop effective control algorithms to create a proper driving control mode so as to suppress vibration and improve driving smoothness.

Conflicts of Interest:
The authors declare no conflict of interest.