Abstract

Hub direct-drive electric vehicles are driving with hub motors, which cause negative vibration directly on wheels and suspensions. To better analyze the problems of dynamic deterioration resulted by coupling excitations from random road and unbalanced magnetic pull (UMP), a cooperative hub direct drive-air suspension (HDD-AS) system is developed, and evaluation indices based on vertical and longitudinal dynamic coupling mechanism are proposed. The vertical and longitudinal characteristics in the time domain under different damping coefficients and motor speeds are demonstrated, and comparison between traditional electric wheel and suspension-integrating models and HDD-AS system model is conducted. Then, the responses of random road excitation and UMP are analyzed, and explanations of vertical and longitudinal characteristics in the frequency domain under different damping coefficients are given. These results provide important insights into the influence of coupling excitations and damping coefficient on the HDD-AS system model. Finally, the effectiveness of the proposed model is verified by the experiment on the test bench.

1. Introduction

With the increasing serious problems of energy shortage and environmental pollution, electric vehicles (EVs) have become the focus of the development on automobile research. EVs with the hub direct-drive system, namely, the distributed driving system, have the advantages of flexible control, efficient transmission, and compact structure, which have been recognized as the ideal configuration of EVs in the future [1].

As the core component of the distributed driving vehicle, the hub direct-drive system is a complex electromechanical-integrating system composed of wheel, in-wheel motor, and suspension. However, the integration of the in-wheel motor into the wheel will increase the unsprung mass. The torque ripple and UMP of the in-wheel motor will lead to the deterioration on vehicle ride comfort, controllability and stability [2], and motor fault [3]. Besides, different from road excitation, high-frequency excitations transmitting from the in-wheel motor not only have influence on the vertical vibration [4] but also on the longitudinal vibration of the electromechanical-integrating system [5].

To suppress the negative vertical vibration generated by the in-wheel motor, numerous research studies have been conducted, and the vertical negative effects of increased unsprung mass and main solutions were summarized. Based on the simulations for influences of the motor system parameters in the frequency domain, the vertical performances of the vehicle with the hub direct-drive system were analyzed [6, 7]. Novel configuration of the multifunctional-integrated electric wheel with the in-wheel suspension was proposed to avoid the influence of the increased unsprung mass [8]. Apart from the structure optimization, multiobjective optimization control for the active suspension system was presented, and the negative vertical vibration of the hub direct-drive system was suppressed [9, 10]. In these studies, structural optimization and multiobjective optimization control for the suspension system were utilized for solving the negative vertical vibration of the hub direct-drive system.

Recently, researchers have shown an increased interest in the longitudinal vibration of the electric wheel system. A longitudinal and torsional coupling dynamic model of the electric wheel was established, and the influence of the rotor position error caused by motor torque fluctuation on the longitudinal vibration of the electric wheel system was analyzed [11]. And, the spatial and temporal distribution characteristics of in-wheel motor radial force which induced the electromagnetic vibration and noise was derived [12]. Besides, an 11-DOF dynamics’ model was developed to study the magnetic force influence on vehicle vertical and lateral coupling dynamics, and the results show that magnetic force makes all dynamic response variables deteriorate in different degrees [13]. To suppress the magnetic force oscillation of the in-wheel motor and negative vibration of the sprung mass, the structural parameters of the in-wheel motor were optimized by the BP neural network [14], and a reliable robust controller for the active suspension system and a hybrid controller for the semiactive suspension system were designed [15, 16]. In these studies, based on the coupling dynamic model of the electric wheel, the reasons of magnetic force oscillation caused by the in-wheel motor were analyzed, and effective suppression methods were proposed.

Although extensive research studies have been carried out on the vertical and longitudinal characteristics’ analysis of the hub direct-drive system, the research on the effect of coupling excitations and hub direct-drive system-controllable parameters are still limited. Therefore, this paper introduced the analysis of the characteristics of the HDD-AS system on damping coefficient and UMP from vertical and longitudinal directions. The paper is organized as follows: literatures are introduced in Section 1. After that, the electromechanical-integrating HDD-AS model is developed in Section 2. Then, the verification of the HDD-AS system model is presented in Section 3. In Section 4, the vertical and longitudinal evaluation indices of the HDD-AS system are proposed, and comparison between the HDD-AS system model and two traditional quarter car models is conducted. Subsequently, the characteristics of the HDD-AS system model under different damping coefficients in the frequency domain are analyzed. Finally, the vertical and longitudinal characteristics of the HDD-AS system on damping coefficient and UMP are concluded.

2. Electromechanical-Integrating Model

The hub direct-drive EV is a complex electromechanical-integrating system, consisting of the in-wheel motor, suspension, and wheel. Due to the advantages of the brushless direct current (BLDC) motor, this paper develops a quarter-integrating model for the hub direct-drive vehicle. The quarter-integrating model includes two subsystems: one is the PMBLDC motor model, providing vertical UMP; the other is the three-DOF quarter car model considering real-time UMP.

2.1. PMBLDC Motor Model

As the key component of the electric drive system, the permanent magnet BLDC (PMBLDC) motor has advantages of high efficiency and low noise [17]. These features make the PMBLDC motor become the optimal choice for the hub direct-drive electric vehicle driving system. In this paper, the outer-rotor PMBLDC motor is adopted to drive the electric vehicle, as shown in Figure 1.

The simplified circuit equation of the outer-rotor PMBLDC motor is expressed as follows:where ua, ub, and uc are the external voltages applied to the three motor electrical connections,R is the equivalent resistance of each stator winding, L is the self-inductance of the stator windings, M is the mutual inductance of the stator windings, ia, ib, and ic are the currents flowing in the stator windings, and ea, eb, and ec are back electromotive forces (EMF).

The equations of EMF are expressed as follows:where,, and are the coefficients of EMF and is the rotor mechanical rotational speed.

According to equations (1) and (2), the electromagnetic torque equation and motor mechanical-motion equation can be expressed as follows:whereis the motor electromagnetic torque, is the motor load torque, is the inertia of the motor, and is the friction coefficient.

The UMP of the outer-rotor PMBLDC motor produces electromagnetic vibration and noise. In order to obtain the UMP, the instantaneous magnetic field distribution in the air gap region of the outer rotor PMBLDC is developed, which is based on the superposition of component fields due to the permanent magnet and armature reaction [18]. After that, the UMP is calculated.

The magnetic field distribution in the air gap region of the outer rotor PMBLDC motor can be expressed as follows:whereis the air gap radius, is the stator angle, is the time, and are radial and tangential air gap magnetic fields in the eccentric state, respectively, and are radial and tangential magnetic fields of the permanent magnet, respectively, and are radial and tangential magnetic fields of the armature reaction, respectively, and is the correction coefficient of magnetic conductivity in the eccentric state.

In the outer-rotor PMBLDC motor, the permanent magnet is attached to the inner surface of the rotor, and its structure is shown in Figure 2.

Assuming that the magnetization in the section of the permanent magnet is uniform, the permeability of the iron core is infinite, and the magnetization direction of the permanent magnet is radial. The magnetic field of the permanent magnet can be obtained through the equivalent remanence method. The equation is as follows:where is the permanent magnet magnetic field, is the permanent magnet remanence, is the vacuum permeability, and is the polar arc coefficient.

When , the radial and tangential magnetic field of the permanent magnet can be expressed as follows [19]:whereis the relative permeability, is the number of pole pairs, is the radius of the permanent magnet position, is the inner radius of the rotor, is the outer radius of the stator, and is the initial angle of the stator.

Radial and tangential magnetic fields of the armature reaction can be expressed as follows:whereandare the amplitude of the radial and tangential magnetic flux of order k, k is the winding distribution order, is the number of single-phase series windings, is the number of slot, ,, and are the vectors of the three-phase winding armature-reaction magnetic-potential vector direction, and ,, and are the angles of the magnetic potential vector.

The magnetic potential-vector star diagram of the three-phase armature winding is shown in Figure 3. The vectors of the three-phase winding armature-reaction magnetic-potential vector direction and the angles of the magnetic potential vector are shown in Table 1.

The amplitude of the radial and tangential magnetic flux can be described as follows:whereis the number of winding turns, is the number of winding branches, is the slot width, is the winding pitch of the motor, and is the slot angle.

The eccentric direction of the stator is assumed to be positive on the y-axis, as shown in Figure 4(b). The correction coefficient of magnetic conductivity in the eccentric state is as follows:whereis the effective eccentricity, and it is defined as follows:whereis the length of the actual air gap, is the thickness of the permanent magnet, and is the eccentricity.

Based on equations of magnetic field distribution in the air gap region, the UMP can be obtained:where is the axial length of the motor. The parameters of the outer-rotor PMBLDC motor are shown in Table 2.

2.2. The Three-DOF Quarter Car Model considering Real-Time UMP

When the motor is running, the motor load torque generated by the HDD-AS system can be calculated in the following equations:where is the rolling resistance coefficient, is the tire-rolling radius, is the tire vertical force, and is the longitudinal driving force.

The physical structure and equivalent model of the HDD-AS system are shown in Figure 5. The mathematical three-DOF quarter car model considering real-time UMP can be expressed as follows:where is the sprung mass, is the stator mass, is the rotor and wheel mass, is the sprung-mass vertical acceleration, is the stator vertical acceleration, is the rotor and wheel vertical acceleration, is the sprung-mass vertical velocity, is the stator vertical velocity, is the sprung-mass vertical displacement, is the stator vertical displacement, is the rotor and wheel vertical displacement, is the air spring force, is the acceleration of gravity, is the damping coefficient, is the bearing stiffness, is the tire stiffness, is the road excitation, P is absolute air-spring pressure, is the atmospheric pressure, is the initial sectional area of the air spring, is the active area conversion rate, is the initial relative air-spring pressure, is the initial air-spring volume, is the active volume conversion rate, and κ is the air polytropic exponent.

Generally, the commonly used tire models include two categories, namely, the Magic Formula model [20] and UniTire model [21]. In these models, the fitting formulas of tire longitudinal force, tire lateral force, and other tire forces can be obtained through the tire test under various working conditions. Although these models provide varieties of tire forces, large amounts of testing data are needed [22, 23]. To avoid this, the rigid ring tire model is selected [24], which is suitable for the theory analysis of low-frequency vibration from 1 to 100 Hz. For the convenience of analyzing the characteristics of the HDD-AS system, the rigid ring tire model is simplified properly by transforming the tire into a mass block with stiffness. Tire vertical force and longitudinal driving force can be expressed as follows:whereis the road friction coefficient, is the tire slip ratio, is the tread-element horizontal stiffness, and is the half length of the tire.

The parameters of the three-DOF quarter car model considering real-time UMP are shown in Table 3.

2.3. The Electromechanical-Integrating Process of the HDD-AS System Model

The electromechanical-integrating model logical diagram based on the PMBLDC motor model and three-DOF quarter car model considering real-time UMP is shown in Figure 6. According to the vehicle reference speed Vref and motor load torque calculated by the motor load equations, the motor current and motor speed can be obtained by the circuit equations. The UMP is calculated by the motor current , motor speed , and eccentricity e through the UMP equations. Then, it is applied to the vertical and longitudinal dynamics’ equations to obtain tire vertical force , longitudinal driving force , and eccentricity e. Finally, the vehicle vibration response is observed through the electromechanical-integrating model.

3. Experiment Validation

In order to validate the HDD-AS system model, a test bench for a quarter vehicle equipped with the electric wheel has been established, as shown in Figure 7.

The electric wheel fixture is supported by the passive suspension. A PMBLDC motor is installed in the wheel, and the original motor controller is selected. The PMBLDC motor and motor controller are powered by the ZhiDou D1 electric vehicle [25]. The test bench is constructed on the basis of an INSTRON 8800 single-channel electrohydraulic servo test system. The vibration exciter of the INSTRON 8800 is an electrohydraulic excitation table, which can impose multiple excitations on wheels, including sine waves, triangular waves, and random road profiles’ excitations classified by ISO.

Three types of sensors are placed on this test bench: two single-axis acceleration sensors, two current sensors, and one pressure sensor. Two single-axis acceleration sensors are placed on the inner and outer side of the bottom electric wheel fixture, to measure the vertical acceleration of the unsprung mass. Two current sensors are placed on two power wires to measure the single-phase current of the PMBLDC motor. The pressure sensor is placed between the passive suspension and suspension guide rod to measure the suspension force. Measured signals are recorded by an LMS SCADAS Mobile SCR05 data acquisition instrument.

Because of the limitations of the experiment conditions, there are differences between the HDD-AS system model and test bench, and adjustments of the simulation model are needed. For the reason that the passive suspension rather than the air suspension is installed in the test bench, the parameters of sprung mass, damping coefficient, and suspension stiffness in the HDD-AS system model should be adjusted.

Considering the sprung mass is replaced by top electric wheel fixture in the test bench, the equivalent sprung mass of the test bench can be obtained through initial suspension force measured by the pressure sensor. The initial suspension force is 2534.5 N, and the equivalent sprung mass is 258.6 kg.

The damping coefficient and stiffness of the passive suspension have been measured before the establishment of the HDD-AS system model. The damping coefficient used in the air suspension is equal to that in the passive suspension. Generally, the air spring (diaphragm type) can be established as the linear model for the requirements of the controller design [26, 27]. In the study by Zhu et al. [28], the linear stiffness of the air spring is defined as the sum of and , and the equations can be expressed as follows:

To match the HDD-AS system model with the test bench, the initial stiffness of the air spring in simulation is assumed to be equal to that of the spiral spring. By substituting the measured stiffness of the spiral spring (18000 N/m) into equation (15), the adjusted can be obtained with unchanged. Comparison of suspension force produced by the spiral spring and adjusted air spring around the balanced position is shown in Figure 8.

To analyze the vibration and power characteristics of the electromechanical-integrating model, excitation, which simulates a driving condition on a random road profile classified by ISO as “A” at a speed of 20 km/h, is imposed by the electrohydraulic excitation table. The comparison of vertical acceleration of the unsprung mass and single-phase current is shown in Figures 9 and 10.

According to Figure 9, there are mainly three peaks in the frequency domain from 1 Hz to 70 Hz. Further studies show that the peak around 5 Hz is caused by the resonance of the vertical vibrations of the sprung mass. The peak around 8 Hz is caused by the resonance of the unsprung mass, which is the combination of the motor and wheel. The peak around 40 Hz is caused by the resonance of the rotor and wheel. The frequency and amplitude of the simulation data roughly coincide with that of the experiment data. However, there is a great difference between the amplitude of experiment and simulation data around 10 Hz. A possible explanation for this may be due to the mass of the bottom electric wheel fixture. According to Figure 10, the simulated single-phase current of the PMBLDC motor is relatively consistent compared with the experiment data as well. Overall, the simulation results are similar to the test results, and so, the reliability of the simulation model is proved to a large extent. As there is structural difference of spiral spring and air spring between the simulation model and test bench, acceptable error exists between the simulation results and test results.

4. HDD-AS System Characteristics on Damping Coefficient and UMP

4.1. Evaluation Indices under Different Damping Coefficients and Motor Speeds

Suspension deflection, tire dynamic load, and sprung mass vertical acceleration are significant indices to evaluate vehicle ride comfort and driving safety. The integration of the in-wheel motor into the wheel not only deteriorates vehicle ride comfort and driving safety but also induces the UMP, which ultimately affects the dynamical performances by transmitting to the wheels. To study the vertical and longitudinal characteristics of the HDD-AS system caused by the coupling excitations from random road surface and UMP, six vertical and longitudinal evaluation indices are proposed.

As common evaluation indices are used to evaluate vehicle performance in suspension research, suspension deflection fd, tire dynamic load Fd, and sprung mass vertical acceleration are selected as the vertical evaluation indices of the HDD-AS system. The equations of vertical evaluation indices can be expressed as follows:

As the controllable parameters of the HDD-AS system, damping coefficient and motor speed are selected to analyze the vertical and longitudinal characteristics of the HDD-AS system, and the RMS of three vertical evaluation indices is shown in Figure 11. Due to the limitation of PMBLDC motor power, the range of motor speed is set as 100–700 rpm. The range of the damping coefficient is set as 1200–3600 N·s/m according to the limitation of the damping ratio (0.25–0.75). The road excitation imposes a random road excitation classified by ISO as “B.”

Figure 11 reveals that the RMS of three vertical evaluation indices rises with the continual increase of motor speed. Higher motor speed produces greater coupling excitations from random road surface and UMP, which result in the increasement of the RMS of vertical evaluation indices. With the increase of the damping coefficient, there have been steady fall on the RMS of suspension deflection and tire dynamic load, while the RMS of sprung mass acceleration rises steadily. The results above have been well known on the characteristics of the two-DOF quarter car model, which also shows the rationality of the HDD-AS system model to some extent.

To study the dynamical performance of the HDD-AS system affected by coupling excitations under the stable driving condition, three longitudinal evaluation indices are proposed. As introduced in equations (9)–(10), the eccentricity e is a significant factor to increase amplitude of UMP. Therefore, the eccentricity is selected to represent the variation characteristics of UMP. Meanwhile, the longitudinal driving force fluctuation σtx and motor load fluctuation σtl are selected to assess the stability of vehicle dynamical performance. The longitudinal evaluation indices are expressed as follows:

The RMS and fluctuation values of three longitudinal evaluation indices are shown in Figure 12.

What is interesting in these figures, with the increase of motor speed, all of three longitudinal indices have two peak values when the motor speed is at 100 rpm and 700 rpm, and the lowest point of RMS and fluctuation values appears when the motor speed is at 200 rpm. This discrepancy could be attributed to the rotation frequency of the motor resonance with the natural frequency of the electric wheel system when the motor speed is at 100 rpm, and the effect of resonance decreases with the motor speed rising up. And, the coupling excitations lead to the peak value of three longitudinal evaluation indices when the motor speed is at 700 rpm.

With the increase of the damping coefficient, the RMS of eccentricity decreases to its lowest point before increasing, and longitudinal driving force fluctuation shows a decrease trend. Ignoring the exception data, motor load fluctuation shows an increase trend with the increase of the damping coefficient. Compared to longitudinal driving force fluctuation, the change of RMS of eccentricity and motor load fluctuation is not obvious in numerical terms, which means that the RMS of eccentricity and motor load fluctuation is not sensitive with damping coefficient and motor speed changed.

4.2. Comparison with Traditional Electric Wheel and Suspension-Integrating Models

In order to simplify the electric wheel model, considerable research studies have been conducted. Some researchers simplified the electric wheel system as a fixed hub-drive system [29]. The motor mass (mmo) and wheel mass (mti) were combined into the unsprung mass. The road excitation and UMP were taken as the external coupling excitation of the system. Further research considered the stiffness of the motor bearing. The motor was divided into two components, namely, rotor and stator, and attached to the wheel and axle, respectively [30]. However, the UMP was still taken as external excitation. These traditional electric wheel and suspension-integrating models are shown in Figure 13.

To analyze the difference between HDD-AS system model and traditional models, the RMS and fluctuation values of evaluation indices on the two-DOF quarter car model and three-DOF quarter car model are divided by that of the HDD-AS system model under the same condition. The road excitation imposes a random road excitation classified by ISO as “B.” Comparison results are shown in Figure 14.

What can be seen from Figure 14 is that there are no obvious differences on RMS ratios of vertical evaluation indices between traditional quarter car models and HDD-AS system model. However, ratios of longitudinal driving force fluctuation and motor fluctuation between the traditional two-DOF quarter car model and HDD-AS system model are larger than those between the traditional three-DOF quarter car model and HDD-AS system model. A possible explanation for these results may be that the traditional two-DOF quarter car model ignored motor bearing stiffness. It means that there is no vibration isolator inside the motor and the vibration of the unsprung mass will be deteriorated. Though it is not obvious, the RMS ratio of tire dynamic load between the traditional two-DOF quarter car model and HDD-AS system model has shown the worst trend of the unsprung mass vibration in Figure 14(b). From the perspective of ratios of longitudinal driving force fluctuation and motor fluctuation, the difference between the traditional two-DOF quarter car model and HDD-AS system model is obvious. Additionally, the abnormal tendency in the RMS ratio of eccentricity and ratio of longitudinal driving force fluctuation is caused by the resonance between the motor and electric wheel system as well. Hence, these findings suggest that the introduction of real-time UMP have an effect on the HDD-AS system model compared to traditional quarter car models to some extent.

4.3. Transmissibility Ratios on UMP and Road Excitation

To analyze the effect of UMP and road excitation contributing to vertical and longitudinal evaluation indices in the frequency domain, the HDD-AS system transmissibility ratios on UMP and road excitation are calculated. The simulated motor speed and damping coefficient is set as 500 rpm and 2300 N·s/m. The results are shown in Figure 15.

What can be clearly seen in this figure is that the influence of road excitation is more than that of UMP. Road excitations have an effect on all evaluation indices in frequency around 1 Hz, 8 Hz, and 40 Hz, which are caused by the resonance of the sprung mass, combination of the motor and wheel, and combination of the rotor and wheel, respectively. The UMP contributes more to evaluation indices in the frequency domain around 40 Hz. According to Figure 15(e), the magnitude of the transmissibility ratio for longitudinal driving force in the frequency domain around 40 Hz is 25 dB, which indicates that UMP have great influence on vehicle dynamical performance rather than vehicle ride comfort.

4.4. Transmissibility Ratios with Different Damping Coefficients

According to Figures 11 and 12, the effect of motor speed on evaluation indices is clear. Apart from the abnormal peak on longitudinal evaluation indices when the motor speed is at 100 rpm, the RMS and fluctuation values of all evaluation indices rise up with the increase of the motor speed. However, there is a difference between the trend of RMS and fluctuation values of evaluation indices with the increase of the damping coefficient. To analyze the effect of different damping coefficients contributing to vertical and longitudinal evaluation indices, the HDD-AS system transmissibility ratios on coupling excitations with different damping coefficients are calculated. The results are shown in Figure 16.

The vibration responses of vertical evaluation indices mainly present in the low-frequency range from 0.1 Hz to 20 Hz. With the increase of the damping coefficient, the magnitude of suspension deflection shows a decrease trend in the frequency domain from 0.1 Hz to 20 Hz. In other words, suspension deflection will be suppressed when the damping coefficient is larger. Because the resonance frequency range of the unsprung mass is around 8 Hz, the magnitude of tire dynamic load shows a decrease trend with the increase of the damping coefficient in this frequency domain. Therefore, tire dynamic load will be suppressed when the damping coefficient rises up as well. The magnitude of sprung mass vertical acceleration is larger when the damping coefficient is smaller in the resonance frequency range of the sprung mass. However, according to Figure 11(c), the RMS of sprung mass vertical acceleration increases with the increase of the damping coefficient. This rather contradictory result may be attributed to the reason that the sprung mass vertical acceleration is influenced not only by the resonance vibration but also the vibration in the frequency domain from 2 Hz to 8 Hz. Therefore, the magnitude of sprung mass vertical acceleration shows a decrease trend with the decrease of the damping coefficient in the frequency domain from 2 Hz to 8 Hz.

The vibration responses of longitudinal evaluation indices are mainly present in the frequency domain around 40 Hz, which is the resonance frequency range for the combination of the rotor and wheel. In this frequency range, the magnitude of longitudinal driving force is smaller when the damping coefficient is larger, which means that increasing the damping coefficient can suppress the longitudinal driving force fluctuation. Due to the small magnitude in the frequency domain around 40 Hz, the effect of changing the damping coefficient on eccentricity and motor load is limited. Overall, these findings are basically consistent with results shown in Figures 11 and 12.

5. Conclusions

In order to better analyze the negative vibration problems of the electric wheel and suspension-integrating system, the HDD-AS system model considering coupling excitations and real-time UMP is established based on vehicle dynamics and electromagnetic mechanics. The effectiveness of the HDD-AS system model has been demonstrated experimentally. Analysis results of the HDD-AS system from vertical and longitudinal directions can be explained as follows:(1)According to the analysis requirements of vehicle ride comfort, driving safety, and dynamical performances, six vertical and longitudinal evaluation indices are proposed. The vertical and longitudinal characteristics of the HDD-AS system under different damping coefficients and motor speeds are demonstrated. And, the difference of the HDD-AS system model considering real-time UMP is presented by comparing with the traditional electric wheel and suspension-integrating models.(2)In the frequency domain, the effect of road excitation and UMP contributing to vertical and longitudinal evaluation indices is analyzed. The responses of vertical and longitudinal evaluation indices are influenced more from road excitation than UMP, and the effect of UMP contributes more to evaluation indices in the frequency domain around 40 Hz. Meanwhile, the reasons of changing trends for vertical and longitudinal evaluation indices under different damping coefficients are explained through the analysis in the frequency domain.

In the future, the HDD-AS system model will be modified by adding the complete rigid ring tire model. Based on the modified model, longitudinal evaluation indices such as sprung mass and unsprung mass longitudinal acceleration will be considered.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was financially supported by the National Natural Science Foundation of China (nos. 51975254 and 51775245).