Inﬂuence of Small Teeth on Vibration for Dual-Redundancy Permanent Magnet Synchronous Motor

: The dual-redundancy permanent magnet synchronous motor (DRPMSM) has high reliability in its applications. Compared with the traditional permanent magnet synchronous motor (PMSM), small teeth were added to the stator in DRPMSM and its vibration has new characteristics. This paper derived the analytic expression of the radial electromagnetic force for PMSM in the air gap based on the Maxwell stress equation and summarized its sources. The simpliﬁed stator teeth model of the traditional PMSM and DRPMSM were established in 12-slot/10-pole (12s10p), and the analytical expression and Fourier expansion of radial electromagnetic force on the teeth were given. The electromagnetic and structural ﬁnite element models of the motor were established, and analyses of modal and harmonic response were carried out. The simulation results show that the stator vibration amplitude of DRPMSM is smaller than that of traditional PMSM at no-load and is also bigger under rated load. Furthermore, in order to guarantee the low amplitude of stator vibration under rated load for DRPMSM, the width of the small teeth should be smaller.


Introduction
The permanent magnet synchronous motor (PMSM) has been widely applied in various fields due to its simple structure, high power density, and reliable operation [1][2][3][4][5][6][7][8][9][10]. In this paper, a dual-redundancy permanent magnet synchronous motor (DRPMSM) with no electromagnetic coupling and low heat coupling between phase windings was evolved from PMSM of 12-slot/10-pole (12s10p). By adding six small teeth in the common slot center of the two adjacent phase windings, the leakage inductance between the different phase windings was basically eliminated and the thermal coupling between the phase windings was reduced by means of appending the heat insulation plate on both sides of the small teeth [11][12][13]. Two sets of three-phase symmetrical windings with star connection and the same phase of electromotive force were arranged on the DRPMSM. Under normal conditions, both of them work at the same time in double redundancy. When one of the three-phase windings malfunctions, its power supply is stopped and another normal three-phase winding continues to work.
The application of DRPMSM as the drive in the military field requires low vibration. Electromagnetic vibration is a main source of motor vibration and scholars worldwide have studied it in depth. Considering the influence of permanent magnet shape on motor vibration, by optimizing the shape of the permanent magnet, the amplitude of motor vibration can be reduced [14]. In reference [15], the finite element model was used to analyze the effect of stator teeth, slot wedge, and armature wedge, and armature winding on the natural frequency of the stator. For the special purpose of PMSM in the operation of variable frequency, the amplitude of vibration can be reduced by eliminating the radial electromagnetic force of specific frequencies [16]. Resonance occurs when the radial electromagnetic force is consistent with the modal frequency of the stator and will lead to great vibration and noise [17]. It is interesting to note that the vibration affecting permanent magnets on the motor is higher than the rated load in PMSM [18]. It has been recognized in [18][19][20][21][22][23][24] that the vibration and noise of a motor can be reduced by changing the pole-slot match to raise the minimum mode of the radial electromagnetic force.
The different parameters affecting the vibration of the motor have been studied in depth through the method of theoretical derivation [14][15][16][17][18][19][20][21][22][23][24] and are also demonstrated in Section 3 but with more details on a different object. On the basis of the abovementioned research, the radial electromagnetic force on the yoke equivalent to the superposition of different order force waves and the influence of extra small teeth on vibration in DRPMSM are analyzed in this paper.

The Structure and the Parameters of DRPMSM
The stator profile are shown in Figure 1. Both of them have the same rotor structure, and the main parameters are shown in Table 1.  Since the slot leakage flux of coils get closed through the small teeth in the common slot center of the two adjacent phase windings, the leakage inductance between the different phase windings are basically eliminated and the insulation panel placed on both sides of the small teeth can weaken the thermal coupling among different phase windings. Assuming that the mechanical angle of the  Since the slot leakage flux of coils get closed through the small teeth in the common slot center of the two adjacent phase windings, the leakage inductance between the different phase windings are basically eliminated and the insulation panel placed on both sides of the small teeth can weaken the thermal coupling among different phase windings. Assuming that the mechanical angle of the small

Derivation of the Magnetic Field in Air Gap
Assuming that the inner surface of the stator is smooth, the permanent magnetic field contains only the odd harmonics. The pole-pairs of the permanent magnetic field R ν are given by where P is the number of the permanent pole when  When the fundamental frequency of the stator current is f , the fundamental angular frequency is ω . The magnetomotive force of the permanent magnetic field in air gap R F (unit: A) is described as where R Rm ν F (unit: A) is the amplitude of R F ; θ (unit: rad) is the angular position of the rotor, measured in mechanical radian from the center of one of the permanent magnet; and the permanent magnetic field frequency R ν f of the pole-pairs R ν can be described by The time harmonic order of the stator current μ are given by The pole-pairs of the magnetic field produced by the armature reaction S ν are given by The magnetomotive force by the armature reaction in air gap S F (unit: A) is described as

Derivation of the Magnetic Field in Air Gap
Assuming that the inner surface of the stator is smooth, the permanent magnetic field contains only the odd harmonics. The pole-pairs of the permanent magnetic field ν R are given by where P is the number of the permanent pole when k R = 0, ν R = P. When the fundamental frequency of the stator current is f , the fundamental angular frequency is ω. The magnetomotive force of the permanent magnetic field in air gap F R (unit: A) is described as where F ν R Rm (unit: A) is the amplitude of F R ; θ (unit: rad) is the angular position of the rotor, measured in mechanical radian from the center of one of the permanent magnet; and the permanent magnetic field frequency f ν R of the pole-pairs ν R can be described by The time harmonic order of the stator current µ are given by The pole-pairs of the magnetic field produced by the armature reaction ν S are given by The magnetomotive force by the armature reaction in air gap F S (unit: A) is described as where F µ,ν S Sm (unit: A) is the amplitude of F S , θ µ,ν S 0 (unit: rad) is the phase angle of F S corresponding to with µ and ν S , and θ (unit: rad) is the angular position of stator, measured in mechanical radian from the center of two big teeth.
Considering the opening of the stator slot, the coordinate origin is the intersection between the line of the big (or normal) teeth center and the stator inner diameter, and the magnetic conductance at air gap λ δ (unit: H −1 ) are given by where λ 0 is the average conductance of the magnetic field in air gap, λ k Z is the amplitude of the k Z th order harmonic of air gap magnetic conductance, and Z is the number of the stator slot. According to Equations (2) and (7), the permanent magnetic field in air gap B Rδ (unit: T) can be described by where According to Equations (6) and (7), the magnetic field by the armature reaction in air gap B Rδ (unit: T) can be described by

The Theoretical Derivation of the Radial Electromagnetic Force
According to the Maxwell stress equations, the radial electromagnetic force in air gap p r can be expressed as where B r and B t are the radial and tangential components of the magnetic field in the air gap, Since the value of |B t | is much smaller than |B r |, p r can be reduced to For the formulas for B Rδ and B Sδ considering only the radial components above, p r can be described by where the radial electromagnetic forces in air gap B 2 Rδ /(2µ 0 ) are produced by the permanent magnetic field alone, B 2 Sδ /(2µ 0 ) by the armature reaction magnetic field alone, and (B Rδ B Sδ )/µ 0 by the interaction between the permanent magnetic field and the armature reaction magnetic field.
According to Equation (12), when the fundamental frequency of the stator current is f and the harmonic of air gap magnetic conductance k Z is equal to 1, the exciting mode number and frequency of p r can be obtained, as listed in Table 2. Table 2. The exciting mode number and frequency of p r .

Magnetic Field in Air Gap of 12s10p
When the fundamental current frequency of 12s10p is f and the harmonic of air gap magnetic conductance k Z is equal to 1, according to Equation (12), the exciting mode number and the frequency times of the radial electromagnetic force produced by permanent magnetic field B 2 Rδ /(2µ 0 ), armature reaction magnetic field B 2 Sδ /(2µ 0 ), and the interaction between permanent magnetic field and armature reaction magnetic field (B Rδ B Sδ )/µ 0 can be obtained, as listed in Tables 3-5.
According to Equation (8), the pole-pair numbers of radial electromagnetic force consist of ν R and Table 3. When ν R =5, 15, 25, 35, ν R + Z =17, 27, 37, 47, and ν R − Z =−7, 3,13, 23.   The mechanical angle of the stator slot opening b was 3.76 • . For Cases 1-4, the mechanical angle of the small teeth a were 0 • , 2 • , 3.76 • , and 5.63 . The amplitude of the radial magnetic field by the permanent magnet and armature reaction alone at rated load obtained by finite element analysis are shown in Tables 6 and 7. The radial electromagnetic force in the air gap will affect the radial electromagnetic force on the stator teeth directly. The mode number of the radial electromagnetic force l is equal to 0, 2, 4, 6, 8... only, according to Tables 3-5. When l is equal to 0, the excitation frequency is far away from the resonance frequency, so it will not cause big vibration noise, and for the 2nd, 4th, 6th, 8th, 10th, 12th and 14th order, its amplitude curve of the radial electromagnetic force in the air gap with the load times can be obtained by the 2D finite element, as shown in Figure 3 (2D diagrams) and 4 (3D diagrams). The four different colors in the diagram correspond to Cases 1-4. In the four cases of the same mode, the amplitudes of the radial forces in the air gap display few differences, as shown in Figures 3 and 4.  It is known from Tables 3-7 that the amplitude of the magnetic field which can produce 2nd-, 10th-, and 12th-order radial electromagnetic forces is larger than that of 4th-, 6th-, and 8th-order forces, so the amplitude of the 2nd-, 10th-, and 12th-order radial forces is higher than that of the 4th, 6th, and 8th order, as shown in Figures 3 and 4. According to Tables 6 and 7, the amplitude of the armature reaction magnetic field is far less than that of the permanent magnetic field, so the amplitude of the 10th-order radial force does not change following the load times (Figure 3e). In Table 7, the magnetic field amplitudes of   With the increase of the load times, the amplitude of the magnetic field in the air gap generated by the armature reaction increases, so the amplitude of the radial electromagnetic force become larger by the armature reaction magnetic field and by the interaction between the permanent magnetic field and the armature reaction magnetic field. Hence, the amplitude of the radial force in the air gap on l = 2, 4, 6, 8, 12 and 14 increases.
The amplitudes of the 2nd-, 12th-, and 14th-order radial forces are higher than that of the 4th-, 6th-, and 8th-order radial forces, and the amplitude of the 10th-order radial force is basically unchangeable.
In the four cases of the same mode, the amplitudes of the radial forces in the air gap display few differences, as shown in Figures 3 and 4.  It is known from Tables 3-7 that the amplitude of the magnetic field which can produce 2nd-, 10th-, and 12th-order radial electromagnetic forces is larger than that of 4th-, 6th-, and 8th-order forces, so the amplitude of the 2nd-, 10th-, and 12th-order radial forces is higher than that of the 4th, 6th, and 8th order, as shown in Figures 3 and 4. According to Tables 6 and 7, the amplitude of the armature reaction magnetic field is far less than that of the permanent magnetic field, so the amplitude of the 10th-order radial force does not change following the load times (Figure 3e). In Table 7, the magnetic field amplitudes of   S 7 ν are as big as that of  S 5 ν , so the amplitudes of the 14th-order radial force are higher than that of the 4th-, 6th-, and 8th-order, as shown in Figures  3 and 4. It is known from Tables 3-7 that the amplitude of the magnetic field which can produce 2nd-, 10th-, and 12th-order radial electromagnetic forces is larger than that of 4th-, 6th-, and 8th-order forces, so the amplitude of the 2nd-, 10th-, and 12th-order radial forces is higher than that of the 4th, 6th, and 8th order, as shown in Figures 3 and 4. According to Tables 6 and 7, the amplitude of the armature reaction magnetic field is far less than that of the permanent magnetic field, so the amplitude of the Assuming that the amplitude of radial electromagnetic force on the stator teeth surface is identical, the lth-order radial force in air gap rl p is described as     r cos l l p p κωt lθ (13) where the amplitude of the lth-order radial electromagnetic force in the air gap is l p , and its frequency times of the lth-order radial force is κ . Because the stator silicon-steel was rigid, the radial electromagnetic force on the stator teeth was calculated by the average value as where the average value of the th l order radial electromagnetic force on teeth number i (in Figure 5) is Til p (unit: N/m 2 ) and Ti θ (unit: rad) is its angular position.
The average teeth forces are given by The amplitudes of the average force on the stator teeth in four cases are shown in Table 8.
The radial electromagnetic force on the stator teeth is transmitted to the stator yoke through Assuming that the amplitude of radial electromagnetic force on the stator teeth surface is identical, the lth-order radial force in air gap p rl is described as p rl = p l cos(κωt + lθ) (13) where the amplitude of the lth-order radial electromagnetic force in the air gap is p l , and its frequency times of the lth-order radial force is κ.
Because the stator silicon-steel was rigid, the radial electromagnetic force on the stator teeth was calculated by the average value as where the average value of the lth order radial electromagnetic force on teeth number i (in Figure 5) is p Til (unit: N/m 2 ) and θ Ti (unit: rad) is its angular position. The average teeth forces are given by where ∆θ n1 = 2π/Z, ∆θ n2 = 2π/Z − b, ∆θ d1 = 2π/Z + a/2 + b/2, ∆θ d2 = 2π/Z − a/2 − 3b/2. The amplitudes of the average force on the stator teeth in four cases are shown in Table 8. The radial electromagnetic force on the stator teeth is transmitted to the stator yoke through the stator teeth. The amplitude of the radial electromagnetic force on the stator yoke with teeth is unequal to 0 and equal to 0 on the stator yoke without teeth, so the radial electromagnetic force acting on the stator yoke p Yl is periodic, continuous, and bounded during the period and conforms to the condition of Fourier expansion.
It can be equivalent to the superposition of different order force waves applied on the inner surface of the stator yoke. p Yl is described as where the ratio of stator teeth width to stator tip width is η, and its value is variable for the different stator teeth. a lk Y is the k Y th order Fourier expansion coefficient of the p Yl . a lk Y is the amplitude of a lk Y , so p lk Y = a lk Y is the amplitude of the radial electromagnetic force on the inner surface of the stator yoke. a lk Y can be obtained by The relation between the mode number and amplitude of deformation for m ≥ 2 can be expressed as [21] where Y ls is the amplitude of the deformation caused by mode m of the radial electromagnetic force on the stator yoke, K s is the coefficient determined by dimensions and structural properties of the stator, and p is the amplitude of mode m in radial electromagnetic force density distribution. In fact, this basic formula cannot be used for accurate calculation and a detailed finite element structural analysis is needed. However, for a rough comparison, it is beneficial to use Equation (22) to study the resultant amplitude of the deformation caused by the radial electromagnetic force. When K s is constant, as shown in Equation (22), Y ls decreases with m 4 .
In this paper, p is considered to be equal to p lk Y and m is considered to be equal to k Y . According to Equations (14)- (22), the Fourier expansion of radial electromagnetic force on the inner surface of stator yoke with different widths can be obtained, as shown in Table 9. In Table 9, with the same structure of the stator, small teeth (yes) stands for the radial electromagnetic force acting on big teeth and small teeth together, and small teeth (no) stands for that on big teeth alone. Table 9 shows that, for the radial electromagnetic force on the stator yoke p Yl , when the order l is equal to 2, 10, 14, the amplitude of |a l2 | obtained by Fourier expansion is higher, so the amplitude of the k Y th (k Y = 2) order radial electromagnetic force on the inner surface of the stator yoke is higher.
When the order l is equal to 4, 8, the amplitude of |a l4 | is higher, so the amplitude of the k Y th (k Y = 4) order radial electromagnetic force on the inner surface of the stator yoke is higher. When the order l is equal to 6, 12, the amplitude of |a l6 | is higher, so the amplitude of the k Y th (k Y = 6) order radial electromagnetic force on the inner surface of the stator yoke is higher. According to Figures 3 and 4 and Table 9, the amplitude of |a l2 |/ 2 4 (l = 2, 10) is much bigger than that of |a l2 |/ 2 4 (l = 4, 6, 8, 12, 14), |a l4 |/ 4 4 (l = 2, 4, 6, 8, 10,12,14), and |a l6 |/ 6 4 (l = 2, 4, 6, 8, 10, 12, 14), so |a l2 | (l = 2, 10) can be considered only below. For p Yl , the amplitude of |a l2 | (l = 2, 10) increases with the small teeth width increasing on condition that the radial electromagnetic force on the big teeth only work. When l = 2, the radial electromagnetic force on the small teeth raises the amplitude of |a l2 |, and the larger the small teeth width, the larger the amplitude of |a l2 |, so the amplitude in Y ls increases. When l = 10, the radial electromagnetic force on the small teeth decreases the amplitude of |a l2 |, and the larger the small teeth width, the smaller the amplitude of |a l2 |, so the amplitude in Y ls decreases. For p Yl (l = 2, 10), its excited radial electromagnetic force frequencies are mainly 2f (600 Hz), as shown in Tables 3-5.

Modal Analysis of Stator
It was important that the modal analysis of the stator structure is completed. This determined whether the resonance occurred on the stator when the frequencies of the radial harmonic components of electromagnetic force coincided with the modal frequencies of the stator.
The motion differential equation determined by the state vector method as follows [24]:
where M is the mass matrices, C is the damping matrices, K is the stiffness matrices, x is the displacement vector, and F is the radial electromagnetic force vector.
Because of the shell, the armature windings and rigid support affected the natural frequency of the mode. Thus, only the stator was considered, and it was assumed that the stator was completely unconstrained. The models of the three-dimensional stators of PMSM in different cases are shown in Figure 6.
The motion differential equation determined by the state vector method as follows [24]: where M is the mass matrices, C is the damping matrices, K is the stiffness matrices, x is the displacement vector, and F is the radial electromagnetic force vector. Because of the shell, the armature windings and rigid support affected the natural frequency of the mode. Thus, only the stator was considered, and it was assumed that the stator was completely unconstrained. The models of the three-dimensional stators of PMSM in different cases are shown in Figure 6.
The material of the stator was DW310-35, its stacking coefficient was 0.95, and its length was 142 mm. The mode of the stator in 12s10p and its resonant frequencies of the 2nd and 4th orders were obtained by modal simulation, as shown in Figure 7. The material of the stator was DW310-35, its stacking coefficient was 0.95, and its length was 142 mm. The mode of the stator in 12s10p and its resonant frequencies of the 2nd and 4th orders were obtained by modal simulation, as shown in Figure 7.

Modal Analysis of Stator
It was important that the modal analysis of the stator structure is completed. This determined whether the resonance occurred on the stator when the frequencies of the radial harmonic components of electromagnetic force coincided with the modal frequencies of the stator.
The motion differential equation determined by the state vector method as follows [24]: where M is the mass matrices, C is the damping matrices, K is the stiffness matrices, x is the displacement vector, and F is the radial electromagnetic force vector. Because of the shell, the armature windings and rigid support affected the natural frequency of the mode. Thus, only the stator was considered, and it was assumed that the stator was completely unconstrained. The models of the three-dimensional stators of PMSM in different cases are shown in Figure 6.
The material of the stator was DW310-35, its stacking coefficient was 0.95, and its length was 142 mm. The mode of the stator in 12s10p and its resonant frequencies of the 2nd and 4th orders were obtained by modal simulation, as shown in Figure 7. Compared with the resonant frequencies of the different small teeth width in Figure 7, it was found that the resonant frequencies of the same mode order decreased weakly in Cases 1-4, and the small teeth had a weak effect on the resonant frequencies in the same mode shapes. Compared with the resonant frequencies of the different small teeth width in Figure 7, it was found that the resonant frequencies of the same mode order decreased weakly in Cases 1-4, and the small teeth had a weak effect on the resonant frequencies in the same mode shapes.

Vibration Analysis of Stator
The radial electromagnetic force on the stator teeth was coupled to the model of the stator by the finite element simulation of the magnetic field, and the vibration simulation of the stator was analyzed when the stator was excited by a radial electromagnetic force at two times the rated frequency (2f = 600 Hz). By observing the average value of vibration velocity on A, B, C, and D in Figure 1 under different teeth widths and load times, the curve can be obtained, as shown in Figure 8.
Four curves with the different colors are used to represent four cases of small teeth width. Solid lines are used to represent the radial electromagnetic force on big teeth alone and the discontinuous curves are used to represent the radial electromagnetic force on big and small teeth together. Comparing the curves of Cases 2-4 in Figure 8, the velocity difference curve of stator vibration with load times can be obtained, as shown in Figure 9.
Comparing the curves of the solid and discontinuous lines in Figures 8 and 9, it can be found that with the larger width of the small teeth, the amplitude of vibration velocity caused by the radial electromagnetic force acting on big teeth alone increases and correspond to the increasing of |a l2 | (l = 2, 10) when the radial electromagnetic force on the big teeth work only.  With increasing load times, the amplitude of the 10th-order radial electromagnetic force on the teeth is almost constant, and the amplitude of the 2nd-order radial electromagnetic force gradually increases.
When the load time is nearly 1, the amplitude of the 2nd-order radial electromagnetic force of  Y ( 2) l p l is bigger, so the amplitude of the vibration velocity produced by the radial electromagnetic force on the big and small teeth together is bigger than that on the big teeth only. The radial electromagnetic force of  Y ( 2) l p l on the small teeth increases the amplitude of  2 a l , and the larger the width of the teeth tip, the larger the increase of velocity.  With increasing load times, the amplitude of the 10th-order radial electromagnetic force on the teeth is almost constant, and the amplitude of the 2nd-order radial electromagnetic force gradually increases.
When the load time is nearly 1, the amplitude of the 2nd-order radial electromagnetic force of  Y ( 2) l p l is bigger, so the amplitude of the vibration velocity produced by the radial electromagnetic force on the big and small teeth together is bigger than that on the big teeth only. The radial electromagnetic force of  Y ( 2)  When the load times are small, the amplitude of vibration velocity produced by the radial electromagnetic force acting on big and small teeth together is smaller than that on big teeth alone, and the larger the width of the small teeth, the larger the velocity difference.
When the load time is equal to 1, the amplitude of vibration velocity caused by the radial electromagnetic force acting on big and small teeth is bigger than that on big teeth alone, and the larger the width of the small teeth, the larger the velocity difference.
Under the light load, the amplitude of the 10th-order radial electromagnetic force of p Yl (l = 10) is much larger than the amplitude of the 2nd-order force of p Yl (l = 2), so the radial electromagnetic force on the teeth is mainly of the 10th order. The radial electromagnetic force of p Yl (l = 10) on the small teeth reduces the amplitude of a l2 (l = 10), and the larger the width of the teeth, the larger the reduction of vibration velocity.
With increasing load times, the amplitude of the 10th-order radial electromagnetic force on the teeth is almost constant, and the amplitude of the 2nd-order radial electromagnetic force gradually increases.
When the load time is nearly 1, the amplitude of the 2nd-order radial electromagnetic force of p Yl (l = 2) is bigger, so the amplitude of the vibration velocity produced by the radial electromagnetic force on the big and small teeth together is bigger than that on the big teeth only. The radial electromagnetic force of p Yl (l = 2) on the small teeth increases the amplitude ofa l2 (l = 2), and the larger the width of the teeth tip, the larger the increase of velocity.
As shown in Figure 8, the amplitude of vibration is largest at the rated load when the load times is from 0 to 1. Hence, in order to ensure the low vibration of the motor overall, it is necessary to guarantee that at a rated load.

Conclusions
Through theoretical derivation and simulation analysis, the main conclusions regarding the stator with different widths of small teeth are as follows: (1) The coefficient amplitudes a lk Y obtained by Fourier expansion of the radial electromagnetic force p Yl are different. The radial electromagnetic force with a lower order and bigger amplitude will produce a higher amplitude of vibration. (2) For p Yl , when l = 2, the radial electromagnetic force on small teeth increases the amplitude of |a l2 |, and the larger the small teeth width, the larger the amplitude of |a l2 |. When l = 10, the radial electromagnetic force on small teeth decreases the amplitude of |a l2 |, and the larger the small teeth width, the smaller the amplitude of |a l2 |. (3) The small teeth have a weak effect on the resonant frequencies in the same mode shapes. (4) The amplitude of vibration velocity caused by the radial electromagnetic force acting on big and small teeth together is smaller than that on big teeth alone when the load times are small, and the larger the width of the small teeth, the larger the velocity difference. When the load times is equal to 1, the amplitude of vibration velocity caused by the radial electromagnetic force acting on big and small teeth together is bigger than that on big teeth alone, and the larger the width of the small teeth, the larger the velocity difference.
When vibration and noise requirements are strict, for example, for underwater military vehicles, the drive of that needs to be as low as possible when running. Based on the above conclusions about DRPMSM, the minimum width of the small teeth should be chosen to obtain the minimum vibration velocity at a rated load to ensure its mechanical strength and electromagnetic function.
Author Contributions: All the authors have contributed significantly. Y.C. and Y.Y. put forward the idea and theoretical verification. Y.S. and Y.Y. built mathematical models and performed the simulation. Y.C. and Y.Y. analyzed the simulation results. Y.Y. made all graphics and wrote the paper.
Funding: This research was funded by the National Natural Science Foundation of China (No. 51377114).