Electromagnetic Hysteresis Based Dynamics Model of an Electromagnetically Controlled Torque Coupling

This paper proposes a novel computationally efficient, easy-to-implement electromagnetic hysteresis based dynamics model of a kind of intelligent electromagnetic torque controlled coupling (EMTC), which, with drag torque under consideration, first models the electromagnetic hysteresis existing in the primary clutch with the classical Preisach model, and then models the transferred torques in the three friction elements of the center coupling in the slipping and locked modes, respectively. The performance of the model is verified by simulation and experiment jointly, which lays the basis for the development of advanced control algorithm.


Introduction
Although electric all-wheel drive vehicles have been a hot topic recently [1][2][3][4][5], the majority of the market share still vests in traditional mechanical AWD vehicles. The last few years have witnessed the evolution of the three generations of EMTC, the high fuel economy and excellent AWD performance make it widely used in many AWD cars, see Figure 1. Controlling of the AWD vehicles has attracted much attention of researchers [6][7][8][9][10], however, a solid and effective model of the AWD effecter remains critical. To address that issue, a lot of literatures have been focused on modeling EMTC along the years.

Introduction
Although electric all-wheel drive vehicles have been a hot topic recently [1][2][3][4][5], the majority of the market share still vests in traditional mechanical AWD vehicles. The last few years have witnessed the evolution of the three generations of EMTC, the high fuel economy and excellent AWD performance make it widely used in many AWD cars, see Figure 1. Controlling of the AWD vehicles has attracted much attention of researchers [6][7][8][9][10], however, a solid and effective model of the AWD effecter remains critical. To address that issue, a lot of literatures have been focused on modeling EMTC along the years. Whole vehicle layout and the EMTC; 1-final reducing drive, 2-assembly of front axle differential, 3-power transfer unit, 4-driving shaft, 5-center coupling, 6-rear axle reducer, 7-assembly of rear axle differential. Torii [11] simply proposed the modeling method of an electronically controlled torque-splitting multi-plate clutch, which is the prototype of today's EMTC. Mohan [12] made an attempt to incorporate the four-wheel drive and AWD systems and their driveline architectures, and to explain the control tactics and philosophies behind the controlling of these traction systems. Cheli [13] introduced a type of driveline torque distribution management device with the corresponding control strategy, which was verified by using a precise vehicle model developed in Simulink. Lee [14] proposed a control system design method for a 4WD vehicle that distributed the front and rear driving torques through a clutchless center limited slip differential, the control system consisted of a control map with three different control modes and vehicle steering estimation method. Piyabongkarn [15] proposed a modeling method for the center coupling and the electronically controlled limited slip differential (ELSD), the controlling of which was designed in hierarchical architecture. Ando [16] came up with two 3D maps to describe the relations between vehicle speed, clutch opening, front-rear wheel speed difference and the driving torque output by EMTC which is quite implementable in the controller. In [17] the advantage of high efficiency of BorgWarner NexTrac TM was introduced, with new friction plate materials and innovatively compact layout, the newly released NexTracTM was capable of transferring larger driving torque with smaller configuration. Vantsevich [18] presented a hypothetical driveline power-dividing unit with the modeling method, the control algorithm of which allowed for flexible controlling of the power flow to the drive wheels.
The EMTC modeled in this paper consists of three friction elements that are conditioned and lubricated by viscous fluid or oil, namely, the EMTC is a wet clutch internally. To model the wet clutch, various researches and modeling methods have been proposed. Ando [19] described the third generation of electronically controlled AWD coupling that contributed higher fuel efficiency by reducing the drag torque due to the viscosity of the lubricant while the clutches were not engaged. Iqbal et al. [20] presented a mathematical model for the estimation of drag torque in multi-disk wet clutches considering the laminar flow in the clearance between the friction disks, and the simulation results validated the model's capability of predicting the drag torque in higher clutch speed range. Völkel [21] investigated the influences of such factors as oil flow rate and design of grooves on the cooling properties and thermal behavior of the wet multi-disk clutches, the flow capacity determinations were also presented, and the scoop grooves conveying effect was finally explained.
Despite the various researches on the modeling of the electronically or electromagnetically controlled center coupling, some key factors that can influence the driving torque governed by EMTC greatly, such as the hysteresis in the electromagnet and the torque allocation in the three friction elements, haven not been taken into consideration or discussed in detail [22], and a more explicit way of dynamically modeling the three friction elements is still missing. Thus, to accurately control the driving torque transferred by EMTC, these issues will be addressed in this paper.
This paper is organized in three parts: First, the electromagnetic force generated by the primary clutch apply system (coil and armature) is modeled utilizing the Preisach model in the environment of Matlab/Simulink; then, the central clutch status will be decided between 'slipping' and 'locked' modes based on the vehicle states detected by sensors, and the driving torques transferred in condition of the two modes will be calculated; and finally, the model built in this paper will be verified through a real-car experiment regarding the slipping and locked modes.

Hysteresis Modeling
The effect of magnetic hysteresis in the electromagnet materials is crucial in predicting the magnetic compression force generated. To model the magnetic hysteresis phenomenon presented in ferromagnetic materials, various modeling methods have been developed. For the electromagnet in EMTC, which is made of soft magnetic materials, the Preisach model is suitable to be implemented using the numerical implementation method mentioned in [23,24]. The electromagnetic force depends Processes 2019, 7, 557 3 of 17 on the input current and the airgap between the electromagnet and the armature, whereas the hysteresis effect of the airgap is negligible, then for the general current input trajectory, the output values are as follows: In Equation (1),γ αβ u(t) = 0, ∀(α, β) ∈ S o (t),γ αβ u(t) = 1, ∀(α, β) ∈ S + (t), the areas denoted by S o (t) and S + (t) can be found in Figure 2, α, β represent the ascending and descending value of the input current in each instant of time, respectively. Equations (2) and (3) calculate the output values in condition of decreasing and increasing input respectively. In the equations, f u(t) denotes the output value on the limiting ascending branch corresponding to the input value u = u(t), in which u(t) is the current time input value; M i , m j indicate the global maximum and minimum inputs of each instant of time, respectively; f M i m j corresponds to the output value at β = m j that lies on the first-order transition curve attached to the limiting ascending branch at the point α = M i . ferromagnetic materials, various modeling methods have been developed. For the electromagnet in EMTC, which is made of soft magnetic materials, the Preisach model is suitable to be implemented using the numerical implementation method mentioned in [23,24]. The electromagnetic force depends on the input current and the airgap between the electromagnet and the armature, whereas the hysteresis effect of the airgap is negligible, then for the general current input trajectory, the output values are as follows: In Equation (1),  Using the method mentioned in [25], as the limit working currents of the electromagnet are Using the method mentioned in [25], as the limit working currents of the electromagnet are α 0 = β 0 = 4.5A, the triangle T containing a grid of points is devised, the intervals between two adjacent points are expressed as ∆α = ∆β = 0.5A. The experimental settings are shown in Figure 2, the forces are measured with an FCS system (manufactured by Beijing FCS Corp.) containing a force sensor, which is mounted on the ferromagnetic armature. The airgap between the armature and the electromagnet can be varied by controlling the up and down movement of the shaft. Using that experimental setting, the forces at each point of inside the triangle T, see Figure 3, can be determined.  For every instant of time, the extrema of the input will be updated, and the output will be calculated. If the input current doesn't lie at the grid point of the triangle, the following interpolation equations can be applied [25]: if the input current extrema ( n n M m ) belongs to a square cell. And if the extrema lie in a triangular cell that is adjacent to the line α β = , where i C can be acquired by matching the values of f αβ in the cell vertices.
As for the effect of airgap, which changes from 0.2 mm to 1.4 mm, on the electromagnetic force, seven calibration points with the step of 0.2 mm are selected, and the electromagnetic forces along the general current trajectory at each point are measured as in the triangle gridding mentioned above. Then the force output at a certain airgap that doesn't lie at the gridding point, can be calculated with the method of second-order polynomial interpolation.

Verification and Inverse Compensation of Hysteresis Model
To verify the effectiveness and validity of the hysteresis model, the electromagnet set in the primary clutch of EMTC is measured for output electromagnetic force with input current. The current in the coil was energized by linear power supply in a range of 0-4.5 A, and the magnetic forces were measured utilizing a torque sensor mounted on the armature. With the experimentally acquired data, as shown in Figure 4, the hysteresis model can be established. The comparison results are shown as follows. For every instant of time, the extrema of the input will be updated, and the output will be calculated. If the input current doesn't lie at the grid point of the triangle, the following interpolation equations can be applied [25]: if the input current extrema (M n m n ) belongs to a square cell. And if the extrema lie in a triangular cell that is adjacent to the line α = β, where C i can be acquired by matching the values of f αβ in the cell vertices.
As for the effect of airgap, which changes from 0.2 mm to 1.4 mm, on the electromagnetic force, seven calibration points with the step of 0.2 mm are selected, and the electromagnetic forces along the general current trajectory at each point are measured as in the triangle gridding mentioned above. Then the force output at a certain airgap that doesn't lie at the gridding point, can be calculated with the method of second-order polynomial interpolation.

Verification and Inverse Compensation of Hysteresis Model
To verify the effectiveness and validity of the hysteresis model, the electromagnet set in the primary clutch of EMTC is measured for output electromagnetic force with input current. The current in the coil was energized by linear power supply in a range of 0-4.5 A, and the magnetic forces were measured utilizing a torque sensor mounted on the armature. With the experimentally acquired data, as shown in Figure 4, the hysteresis model can be established. The comparison results are shown as follows.  In Figure 5, the performance of the proposed model in normal scenario, where the amplitude of the current input first decreases then increases, proves the effectiveness and accuracy in predicting the corresponding force output with respect to the current input. While in Figure 6, the result shows the effective prediction performance when the airgap doesn't lie at the gridding point. Now it's possible to design an inverse controller to offset the hysteresis effects and increase the accuracy of the electromagnetic force commanded.  In Figure 5, the performance of the proposed model in normal scenario, where the amplitude of the current input first decreases then increases, proves the effectiveness and accuracy in predicting the corresponding force output with respect to the current input. While in Figure 6, the result shows the effective prediction performance when the airgap doesn't lie at the gridding point. Now it's possible to design an inverse controller to offset the hysteresis effects and increase the accuracy of the electromagnetic force commanded.  In Figure 5, the performance of the proposed model in normal scenario, where the amplitude of the current input first decreases then increases, proves the effectiveness and accuracy in predicting the corresponding force output with respect to the current input. While in Figure 6, the result shows the effective prediction performance when the airgap doesn't lie at the gridding point. Now it's possible to design an inverse controller to offset the hysteresis effects and increase the accuracy of the electromagnetic force commanded.   In Figure 5, the performance of the proposed model in normal scenario, where the amplitude of the current input first decreases then increases, proves the effectiveness and accuracy in predicting the corresponding force output with respect to the current input. While in Figure 6, the result shows the effective prediction performance when the airgap doesn't lie at the gridding point. Now it's possible to design an inverse controller to offset the hysteresis effects and increase the accuracy of the electromagnetic force commanded.  Based on the theories and formulas in [25], the corresponding current in the case of decreasing target force output, m l < β = i(t + ∆t) < m l+1 , can be acquired, and for increasing force output, M l+1 < α = i(t + ∆t) < M l , the corresponding input current can be obtained by the follows equations: Upon obtaining the accurate compensation current to be input in the next instant of time, an open-loop compensation of electromagnetic hysteresis can be devised to achieve the desired output profile, as in Figure 7. Based on the theories and formulas in [25], the corresponding current in the case of decreasing target force output, be acquired, and for increasing force output, , the corresponding input current can be obtained by the follows equations: Upon obtaining the accurate compensation current to be input in the next instant of time, an open-loop compensation of electromagnetic hysteresis can be devised to achieve the desired output profile, as in Figure 7. The verification of the performance of the inverse model can be conducted by comparing the output of the electromagnet and the input calculated force that can also be seen as the output of some anhysteretic material. The comparison of the measured data and the input to the model can be seen in Figure 8. The verification of the performance of the inverse model can be conducted by comparing the output of the electromagnet and the input calculated force that can also be seen as the output of some anhysteretic material. The comparison of the measured data and the input to the model can be seen in Figure 8. Based on the theories and formulas in [25], the corresponding current in the case of decreasing target force output, be acquired, and for increasing force output, , the corresponding input current can be obtained by the follows equations: Upon obtaining the accurate compensation current to be input in the next instant of time, an open-loop compensation of electromagnetic hysteresis can be devised to achieve the desired output profile, as in Figure 7. The verification of the performance of the inverse model can be conducted by comparing the output of the electromagnet and the input calculated force that can also be seen as the output of some anhysteretic material. The comparison of the measured data and the input to the model can be seen in Figure 8. As is shown in the figure, there are some minor errors detected between the reference values and the measured values. The main possible error sources can be measurement error, interpolation error of the model and computational delay, which tend to be easily eradicated by increasing the number of reference points in the Preisach triangle. For this reason, it is safe to say that the output corresponding to the current input converted follows excellently the target force generated by the control algorithm.

The Electromagnetic Actuator Model
The schematic diagram of the electromagnetic effector is illustrated as follows in Figure 9  As is shown in the figure, there are some minor errors detected between the reference values and the measured values. The main possible error sources can be measurement error, interpolation error of the model and computational delay, which tend to be easily eradicated by increasing the number of reference points in the Preisach triangle. For this reason, it is safe to say that the output corresponding to the current input converted follows excellently the target force generated by the control algorithm.

The Electromagnetic Actuator Model
The schematic diagram of the electromagnetic effector is illustrated as follows in Figure 9 [26]: Namely, where, and the airgap can be expressed as:

The Torque Transfer Model
The two clutches and the third friction element can work in two different modes apart from the 'separated' mode, like other kinds of clutches, slipping and locked modes. The fact that the torques Set K spring as the stiffness of the spring, F k0 the original spring force, c spring the spring damping and m 0 the mass of the armature. Then the dynamics equation of the electromagnetic effector can be expressed as follows: Namely, where, and the airgap can be expressed as:

The Torque Transfer Model
The two clutches and the third friction element can work in two different modes apart from the 'separated' mode, like other kinds of clutches, slipping and locked modes. The fact that the torques transferred in these two modes are different even with other conditions remaining the same, makes the precise modeling of the torque transmitting process of great significance [26].
The driving plate of the primary clutch is connected to the housing of the input shaft (interference fit with the input shaft), the engaged plate to the ball-cam mechanism by spline connection, the driving cam of which is then connected to the output shaft by spline, thus, under the static and quasi-static circumstances, there is no rotational speed difference between the driving and engaged plates of the ball-cam mechanism, which in turn leads to the equality of the speeds of the engaged plate of the primary clutch and the output shaft. The driving and engaged plates of the secondary clutch are inner and outer spline-connected to the housing of the input shaft and the output shaft, respectively. And the driving plate of the third friction element is fixed on the input shaft, the engaged plate to the driving plate of the ball-cam mechanism. In a nutshell, the three friction elements are zero-speed-differently connected to the input and output shafts at the both ends of them, the whole system can be seen as a parallel structure, while the primary clutch and the third element comply to series connection of the output ends at the ball-cam mechanism.
The layout of the three torque transfer elements is shown in Figure 10.
Processes 2019, 7, x FOR PEER REVIEW 8 of 17 transferred in these two modes are different even with other conditions remaining the same, makes the precise modeling of the torque transmitting process of great significance [26]. The driving plate of the primary clutch is connected to the housing of the input shaft (interference fit with the input shaft), the engaged plate to the ball-cam mechanism by spline connection, the driving cam of which is then connected to the output shaft by spline, thus, under the static and quasi-static circumstances, there is no rotational speed difference between the driving and engaged plates of the ball-cam mechanism, which in turn leads to the equality of the speeds of the engaged plate of the primary clutch and the output shaft. The driving and engaged plates of the secondary clutch are inner and outer spline-connected to the housing of the input shaft and the output shaft, respectively. And the driving plate of the third friction element is fixed on the input shaft, the engaged plate to the driving plate of the ball-cam mechanism. In a nutshell, the three friction elements are zero-speed-differently connected to the input and output shafts at the both ends of them, the whole system can be seen as a parallel structure, while the primary clutch and the third element comply to series connection of the output ends at the ball-cam mechanism.
The layout of the three torque transfer elements is shown in Figure 10. The modes transition strategy is determined in Figure 11 [26]. As regards the slipping mode, with the increase of the compression force of the electromagnet, the locking mode steps in if the torque capacity of the clutch cap T exceeds the torque required when being locked lock T , and it only happens when the speed difference fades gradually to zero. One thing worth mentioning is that, generally during the mode transition, shock occurs because of the sharp change of the torque transferred, as the static friction coefficient is slightly larger than that of the slipping friction, which is also an essential factor that should be taken into consideration in the control algorithm. The modes transition strategy is determined in Figure 11 [26]. As regards the slipping mode, with the increase of the compression force of the electromagnet, the locking mode steps in if the torque capacity of the clutch T cap exceeds the torque required when being locked T lock , and it only happens when the speed difference fades gradually to zero. One thing worth mentioning is that, generally during the mode transition, shock occurs because of the sharp change of the torque transferred, as the static friction coefficient is slightly larger than that of the slipping friction, which is also an essential factor that should be taken into consideration in the control algorithm. is the torque capacity in lock mode.

Torque Transfer Model for Slipping Mode
In condition of slipping mode, the torques where, ki μ are the kinetic friction coefficients of the primary, secondary clutch and the third friction element, respectively, 1 2 , n n F F are the compression force applied on the primary and secondary clutches, 1 1 , z z are the number of friction pairs of the two clutches, and i R are the effective radiuses of the three elements. Neglect the moment of inertia of the driving plate in the ball-cam mechanism, the compression force applied on the secondary clutch and the third friction element is, In the equation above, in T is the torque transferred by the ball-cam mechanism, R is the working radius of the mechanism, and λ is the amplifier ratio of the mechanism that can be measured experimentally and seen as a constant. Substituting those formulae, the torques transferred and torque capacities of each friction element can be obtained: Based on the empirical values, it is also plausible to set 1.02 · T cap < T lock as the threshold of the clutch's transition into slipping mode, as is illustrated in Figure 10, where, according to Newton's Law, F n2 = F n3 , ∆ω = ω f − ω r , ∆ω pre is the speed difference of the last step, T cap_S is the torque capacity of the clutch in slipping mode, and T cap_L is the torque capacity in lock mode.

Torque Transfer Model for Slipping Mode
In condition of slipping mode, the torques T ci (i = 1, 2, 3) transferred by the three friction elements are kinetic friction torques T ki that can be calculated as follows: where, µ ki are the kinetic friction coefficients of the primary, secondary clutch and the third friction element, respectively, F n1 , F n2 are the compression force applied on the primary and secondary clutches, z 1 , z 1 are the number of friction pairs of the two clutches, and R i are the effective radiuses of the three elements. Neglect the moment of inertia of the driving plate in the ball-cam mechanism, the compression force applied on the secondary clutch and the third friction element is, In the equation above, T in is the torque transferred by the ball-cam mechanism, R is the working radius of the mechanism, and λ is the amplifier ratio of the mechanism that can be measured experimentally and seen as a constant. Substituting those formulae, the torques transferred and torque capacities of each friction element can be obtained: where µ si denotes the kinetic friction coefficients of the primary, secondary clutch and the third friction element, respectively.
Thus, the total torque transferred and the torque capacity are,

Torque Transfer Model for Locked Mode
When the clutches are locked in the sense of speed, the static friction torques transferred T li can vary and shift among the three friction elements with the total torque changing, considering the parallel structure of the clutches shown in previous sections. The torque distribution of the three friction elements can be divided into three conditions, as illustrated in Figure 12.
where si μ denotes the kinetic friction coefficients of the primary, secondary clutch and the third friction element, respectively. Thus, the total torque transferred and the torque capacity are,

Torque Transfer Model for Locked Mode
When the clutches are locked in the sense of speed, the static friction torques transferred li T can vary and shift among the three friction elements with the total torque changing, considering the parallel structure of the clutches shown in previous sections. The torque distribution of the three friction elements can be divided into three conditions, as illustrated in Figure 12.
2. The capacity of the third element, In Figure 11, the distribution is shown as an example under the circumstances of T cap1 > T cap3 , while the contrary condition also stands a chance in the real car. To explain the distributions in detail, we set T li , T cap_li as the static friction torque and the torque capacity in locked mode for the primary clutch, the secondary clutch and the third friction element, respectively.

1.
The torque transferred is small enough for the three elements to divide equally: 2.
The capacity of the third element, T cap_l3 , passes through as it doesn't exceed T lock /3, and the rest will be taken over by the other two clutches: 3.
The primary clutch and the third friction element both transmit the maximum of them, and the secondary clutch takes care of the rest: In condition of T cap1 > T cap3 , a switch of T l1 and T l3 will meet the requirement and the equations are quite similar to those listed above. The total torque transferred and the total torque capacity are: Based on the previous discussions, the torque distribution logic can be devised.

Friction and Drag Torque Model
The friction coefficient between two friction discs can vary due to the variation of the factors such as the temperature inside the clutches. Modeling the friction situation is a challenging task. Despite the difference of the friction types, the Lugre model has made it possible to unify the representation of these types in one formula [27], in which the friction coefficient is the function of the bristle deflection z and the speed difference between the friction surfaces v. .
where, σ 0 and σ 1 are the tangential direction stiffness and damping coefficient, σ 2 the viscous friction coefficient, µ c the Coulomb friction coefficient, µ p the peak friction coefficient, and v s the Stribeck velocity. The six parameters are typically bristle-related, and therefore can be experimentally determined. Then the static and kinetic friction coefficients can both be identified with varied speed difference between the friction surfaces. It was discovered in experiments that there is torque transferred in the clutches even in disengaged phases. This uncontrollable and undesired torque is defined as drag torque, which is induced by the resistance of the lubricant fluid filled in the clutches to the relative rotation of the friction discs. The factors that influence the drag torque of the multi-plate clutch include the number of discs, the clearance between two consecutive discs, lubricant flow rate, internal temperature of the clutch and the radius of the friction discs.
The calculation of drag torque can be expressed as follows, where, N denotes the number of friction pairs, µ the dynamic viscosity of the clutch fluid, ω the relative rotational speed of the clutch discs, h the axial clearance between two consecutive discs, R o the equivalent radius and R 1 inner radius of the discs. The detailed parameter selection and calculation were given in [28]. And the drag torque variation with respect to the relative rotation speed is shown as s in Figure 13.
where, N denotes the number of friction pairs, μ the dynamic viscosity of the clutch fluid, ω the relative rotational speed of the clutch discs, h the axial clearance between two consecutive discs, o R the equivalent radius and 1 R inner radius of the discs. The detailed parameter selection and calculation were given in [28]. And the drag torque variation with respect to the relative rotation speed is shown as s in Figure 13. Then the overall driving torque the EMTC transfers is expressed as, where, 0 η = denotes the clutch is in disengaged mode, and 1 η = means the clutch is engaged.

The Experiment Settings
When the EMTC model is connected with a whole-car model, the transferred lock torque lock T will be decided by the driving torque input from the engine and the road condition jointly. Then according to Figure 1, where, , Then the overall driving torque the EMTC transfers is expressed as, where, η = 0 denotes the clutch is in disengaged mode, and η = 1 means the clutch is engaged.

The Experiment Settings
When the EMTC model is connected with a whole-car model, the transferred lock torque T lock will be decided by the driving torque input from the engine and the road condition jointly. Then according to Figure 1, ω cr (27) where, J w f , J wr are the net integrated moments of inertia of the front and rear drive axle and the wheels, respectively; J s f , J sr are the moment of inertia of the driving and engaged shaft of the EMTC; T trans is the output torque of the transmission and i 0 the gear ratio for the front retarder; T c f , T cr and ω c f , ω cr denote the rotational speed and driving torques of the input and output shaft of the EMTC, respectively; T f f , T f r indicate the rolling resistance torques originated from the revolving drums and exerted on the front and rear wheels.
To assure the performance of the locked mode, the gear ratios of the front and rear axles against the center coupling should be equal, Substituting (28) into (27), we get where, And to improve the robustness and the practicability of the model, a Torsen reversal damping model is added to the output lock torque, where, K s , D s denote the equivalent stiffness and damping coefficient of the driving system. To give a reference for the input current of the EMTC, a two-input control algorithm of estimated output torque is designed, where T cal (t) is the calculated reference value of the transferred torque by EMTC, thr indicates the throttle opening, and v x denotes the vehicle longitudinal velocity. The calculated value and the measured value in the experiment will be compared to verify the effectiveness of the dynamic model.

Real-Car Experiment
The car used in this experiment is an FWD (front-wheel driving) based real-time AWD sport utility vehicle, the system configuration is as shown in Figure 1, and the main parameters are listed in Table 1. The real-car experiment is conducted on a Mustang revolving drum test table, see Figure 14, which consists of eight pairs of independent revolving drums and one control system, it is able to adjust both the rotating speeds and resistant torques of the drums while measuring the driving torques working on them. In order to test the output torque of the EMTC in slipping and locking modes, the speed differences of the two axles can be set to a considerable positive constant and some value close to zero, respectively. torques working on them. In order to test the output torque of the EMTC in slipping and locking modes, the speed differences of the two axles can be set to a considerable positive constant and some value close to zero, respectively. To eliminate the spikes and some unknown perturbations [9], two first-order responses ( 1 2 , τ τ ) are implemented for the control of the current [10], as shown in Figure 15. Then the test results for slipping and locked modes are illustrated in Figures 16 and 17, respectively. To eliminate the spikes and some unknown perturbations [9], two first-order responses (τ 1 , τ 2 ) are implemented for the control of the current [10], as shown in Figure 15. torques working on them. In order to test the output torque of the EMTC in slipping and locking modes, the speed differences of the two axles can be set to a considerable positive constant and some value close to zero, respectively. To eliminate the spikes and some unknown perturbations [9], two first-order responses ( 1 2 , τ τ ) are implemented for the control of the current [10], as shown in Figure 15. Then the test results for slipping and locked modes are illustrated in Figures 16 and 17, respectively. Then the test results for slipping and locked modes are illustrated in Figures 16 and 17, respectively.   Figure 16 shows the states of the vehicle and the torques transferred in condition of slipping mode. It can be seen that the difference between the front and rear wheel speeds almost remained constant throughout the experiment time; The transferred torques in each of the three elements are proportional to the overall transferred torque, and the calculated torque from the model matches excellently the measured torque to a large extent; Figure 17 denotes those states and torques in condition of locked mode, in which the torques transferred in the three friction elements display the three situations depicted in the chapter above, and the calculated torque follows the measured values precisely. The effectiveness of the dynamic model can be verified through all those results together.

Conclusions and Future Work
An EMTC dynamics model has been developed considering the electromagnetic hysteresis effect in the primary clutch and the torque distribution when the clutches are in disengaged, locked and slipping modes, respectively, which lays the groundwork for advanced control algorithm for EMTC. The simulation and experiment results verified the validity and effectiveness of the model.
The contributions of this paper can be summarized into two points: First, the model proposed in this paper takes the electromagnetic hysteresis into consideration, which can have a significant effect on the accuracy of the calculation of the output driving torque; and second, three working modes of EMTC are classified and the corresponding driving torques calculation is given mathematically.
As an important part of the fail-safe system, the thermal factor of the clutches is essential to the control effectiveness of the EMTC, especially in the trapped situation. Thus, the thermal model of EMTC is also worth investigating. In addition, advanced control strategies based on mechanical 4WD system aiming at improving the mobility and stability of the vehicle are essential for the development of a complete 4WD system, all of which will be discussed in later literature.
Author Contributions: J.F. and Z.Q. conceived and designed the experiments; J.Z. performed the experiments; J.F. and Z.Q. analyzed the data; S.C. supervised the whole project; Z.L. contributed analysis tools; J.F. wrote the paper.
Funding: This research was funded by National Key R&D Program of China: 2017YFB0102600. And the APC was funded by Zhiquan Qi.

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