Abstract
A spatial kinetics model for supercavitating vehicles including the coupled motion of pitch, yaw, and rolling was presented. Mathematical simulations were performed to investigate dynamic performance of supercavitating vehicles. An intrinsic conic-like oscillation in the motion of supercavitating vehicles was discovered. A lakebed experiment was carried out to validate the oscillation. Good agreement is achieved between simulation and experiment in terms of the amplitude and frequency of attitude angles. The steady phase difference between the motions in different planes is also accurately captured which is about . The conic-like oscillation is attributed to the occurrence of the tail-slap motion in the vertical and horizontal planes and the steady phase difference.
1. Introduction
With most surfaces encompassed by supercavities, the drag force of supercavitating vehicles is greatly reduced compared to conventional fully wetted vehicles. This results in tenfold increase in speed. Logvinovich [1, 2] made the pioneering work on the theory of supercavity flow. Semenenko [3], Savchenko [4], and Vasin and Paryshev [5] further expanded this theory and carried out a set of experimental investigations on supercavity shape, ventilated cavitation, hydrodynamic force, and control strategy.
However, challenges still exist in kinetics model and motion control of supercavitating vehicles due to their unique characteristics, for example, time-delay effect of supercavity, and nonlinearity of planing force [6, 7]. Therefore, studies on hydrodynamic characteristics, kinetics model, and control strategy for supercavitating vehicles are increasing. Nguyen et al. [8–10] established a 6-Degrees-of-Freedom (6-DOF) kinetics model of supercavitating vehicles considering the time-delay effect and the nonsteady planing. The application of inner-loop and outer-loop control schemes in the maneuvering of the vehicle is studied. They concluded that the hysteresis and nonlinearity in the system are mutually related and pose a significant influence on the motion of the vehicle. Nguyen et al. [11] also used a moving Chimera grid method to predict the real-time motion of water entry bodies by combining the 6-DOF rigid body motion model and the numerical calculation of the multiphase flow field, by which the coupled effect of supercavity and moving body can be obtained. S. Kim and N. Kim [12] performed integrated dynamics modeling of supercavitating vehicles and established the 6-DOF equations by defining the hydrodynamic forces and moments. The simulations also demonstrated the instability of supercavitating vehicles caused by the asymmetric immersion depth of coupled horizontal or vertical fins. Mirzaei et al. [7] held the idea that the existing planing force models weaken the nonlinear interaction among the solid, liquid, and gaseous phases and are often too simple and therefore inaccurate. They established the 6-DOF motion model of supercavitating vehicles by redefining the spatial planing force and developed a hybrid control scheme based on the online planing force identification which can increase the stabilization of supercavitating vehicles.
To validate the theoretical method, experimental investigations of supercavitating vehicles are also carried out. Mokhtarzadeh et al. [13] studied the hydrodynamic characteristic of disk and cone cavitator and provided the empirical correlation to calculate the lift and drag force acting on the cavitator. Yen et al. [14] conducted experiments and concluded the method to determine the planing force of cylinder on the free surface of water. Sanabria et al. [15] established an experimental system to validate the mathematical models and control systems of supercavitating vehicles by utilizing a scaled-down model, of which the cavitator is fixed in the water tunnel, and the after-body and fins are all rotatable. The experimental system is demonstrated to be effective in assessing the mathematical model of supercavitating vehicles including actuation and time delay and be suitable for evaluating the real-time control and guidance algorithms. Choi et al. [16] experimentally studied the hydrodynamic characteristics of the real size supercavitating vehicle model in a large cavitation tunnel. It is found that the drag coefficient of the vehicle is rarely changed with the Reynolds number and the drag force is largely determined by the level of the body covered by the supercavity.
With the achievements obtained in the kinetics modeling and experimental investigation, attentions are increasingly drawn by the control strategy of supercavitating vehicles. Kirschner et al. [17] analytically studied the hydrodynamic forces acting on the control plane of supercavitating vehicles. A control strategy considering the angle of rudder and the position of the supercavity is proposed. This method was employed by Goel [18] to investigate the control strategy of supercavitating vehicles. Euteneuer [19] presented a layout scheme of supercavitating with cruciform orientation rudders. A nonlinear model to design the control law was developed. It is concluded that the decrease of the planing force is helpful to both the sailing control and the drag-reduction of the vehicle. Hassouneh et al. [20] emphasized the significance of the time-delay effect of advection during the modeling of supercavitating. The delayed time can lead to destabilization of the motion of supercavitating vehicles and cause the periodic motion. It is asserted that the control law designed by neglecting the delayed time is useless. Dzielski and Kurdila [21] parametrized the vertical motion equations by the geometrical and sailing features of supercavitating vehicles and derived the benchmark problem for the control of the vehicles and presented several nonlinear methods to design the active control system. Kirschner et al. [22] employed a discontinuous slope to describe the hydrodynamic force acting at the hull of supercavitating vehicles. The time-delay effect is regarded as a disturbance when investigating the behavior of the system under various conditions. Vanek et al. [23, 24] developed a linear, parameter-varying controller and improved the performance by considering the coupling of governing equations. The inner-loop control technique is proposed to handle the switched, time-delay dependent behavior of the system.
The intrinsic unsteady motion characteristic of supercavitating vehicles is widely accepted and is called the tail-slap phenomenon. Lin et al. [25] studied the bifurcation behavior with respect to a quasi-static variation of the cavitation number. The tail-slap phenomenon is identified as the consequence of a Hopf bifurcation followed by a grazing event. Considering dominant nonlinearities associated with planing forces, Lin et al. [26] presented a dynamic model. This provided a theory to interpret the tail-slap phenomenon as a limit cycle motion.
However, only the motion of supercavitating vehicles in a single plane (vertical or horizontal) is modeled when investigating hydrodynamic forces, motion model, and sailing control. The strong coupled relationship between the vertical and horizontal motion, caused by the time-delay effect, is typically ignored. This paper presents a spatial kinetics model for supercavitating vehicles. The coupled motion (pitch, yaw, and rolling) and the planing forces are explicitly included in this model. A lakebed experiment is then performed to validate the developed spatial kinetics model and the simulation results.
2. Coordinate Systems
Supercavitating vehicles are composed of a cavitator, a conic section, a cylindrical section, vertical rudders at the tail, and a nozzle, as shown in Figure 1. During sailing, pitch motion is controlled by the disk cavitator; yaw motion is adjusted by the vertical rudders at the tail and rolling motion is regulated by the vertical rudders through operating differential angle.
To describe the 6-DOF motion of supercavitating vehicles, a set of coordinate systems are defined, and the positional relations are depicted in Figure 2.(i)Earth-fixed coordinate system is an inertial coordinate system and its origin is associated with the launching site. is defined as the horizontal plane; lies on the direction of initial motion.(ii)Body-fixed coordinate system is attached to the vehicle and its origin is the vehicle’s mass center; the axis is along the vertical axis of the vehicle.(iii)Translational coordinate system : the axes are all parallel to the corresponding axes of the earth-fixed coordinate system, and its origin locates at the mass center.(iv)Velocity coordinate system : the origin is the same as , and the axis points to the velocity direction of the mass center.
To express the attitude and velocity of the vehicles, Eulerian angles and hydrodynamic angles are defined as follows.(i)Angle of pitch is the intersection angle between the axis and the plane .(ii)Angle of yaw is the intersection angle between the axis and the projection of the axis onto the plane .(iii)Angle of rolling is the intersection angle of the axis and the vertical plane including .(iv)Trajectory inclination angle is the intersection angle between the axis and the plane .(v)Trajectory deflection angle is the intersection angle between the axis and the projection of the axis onto the plane .(vi)Angle of attack is the intersection angle between the axis and the projection of the axis onto the plane .(vii)Angle of sideslip is the intersection angle between the axis and the plane .
For these angles, , , describe the spatial attitude of the vehicle, and describe the velocity direction of the mass center, and and are closely associated with hydrodynamic forces of the vehicle.
3. Forces on Supercavitating Vehicles
For ventilated supercavitating vehicles, hydrodynamic forces only act on wetted parts including the cavitator, the partial hull surface (planing zones), and the vertical rudders. This section briefly describes the definition of each hydrodynamic force.
3.1. Forces on Cavitator
The cavitator is employed as a head rudder. The force acting on the cavitator is largely determined by the deviation angle of the cavitator, the angle of attack , the pitch rate in the vertical plane, and the sailing speed [13]. Under the conditions of nonyaw and nonrolling, the force on the cavitator is defined as along the axis if , , and . However, when , , and , the force is divided into drag and lift , as shown in Figure 3. In the body-fixed coordinate system, the magnitudes of the pitch moment produced by drag and lift are both and with the opposite direction. Therefore, the angle of attack provides no effect on the pitch moment, but lift and drag forces are generated.
When the pitch rate and , , an additional velocity and an additional angle of attack are generated at the vertical plane as shown in Figure 4. The drag along the direction of the composite velocity is , and the lift along the perpendicular direction of the composite velocity is . The projections of these two forces onto the axis are both except at the reversed directions. Thus, no additional lift or moment on the cavitator is produced by .
Accordingly, the lift force , pitch moment and lateral force on the cavitator can be described as
3.2. Forces on Vertical Rudders
The lateral force , yaw moment , and rolling moment of the upper-vertical rudder arewhere and represent the additional rudder angles generated by the rolling motion. The lateral force , yaw moment , and rolling moment from the lower-vertical rudder are
Equations (3) indicate the hydrodynamic forces on the vertical rudders that are determined by the parameters of , , , and . The offset between the vehicle and the supercavity changes the wetted length of the rudders during the maneuvering motion. This further influences the hydrodynamic forces on the rudders [15]. Therefore, , , , and can be obtained explicitly once the relative position between the rudders and the supercavity is confirmed.
For weak maneuvering supercavitating vehicles, the trajectory is approximately a straight line [27]. The offset between the center of the vehicle and the supercavity at the tail becomes negligible. This leads to a small change in the wetted length of vertical rudders. Due to the slight change of wetted zone of the rudders, and are regarded as constants. In addition, the hydrodynamic coefficients and are obtained through Computational Fluid Dynamics method in this paper.
3.3. Planing Forces
The planing forces are hydrodynamic forces acting on the wetted hull [5]. For the wetted zone of the vehicle, the lift force , pitch moment , lateral force , yaw moment , and rolling moment arewhere is positive when the acting point of the planing force is behind the mass center and vice versa. Typically, the acting location of the damping forces (caused by the planing zone) is at the surface of cylindrical section of the vehicle. Moreover, the arm of the damping force is represented by product of (a correction coefficient slightly less than 1 and manifested by Computational Fluid Dynamics method) and . and are the hydrodynamic position force and damping force acting at the planing zone, respectively, and are obtained by Computational Fluid Dynamics method.
4. Relative Position between Vehicle and Supercavity
Hydrodynamic forces acting on the rudders and hull are associated with the offset between the vehicle and the supercavity. The offset can be calculated from the principle of independent expansion of the supercavity sections [1, 2]. The time-delay effect is unavoidable and is the interval where the disturbance at the cavitator propagates to the supercavity at the tail. Additionally, the position and size of supercavity at an arbitrary section are only determined by the previous states of the cavitator. The transformation matrix from to at the time iswhere , , , , , , , , and .
The transformation matrix from to at the time is
The transformation matrix from to at the time is
The transformation matrix from to at the time iswhere , , , , , , , , and .
In , the center of the cavitator and tail are () and (), respectively. The tail center in at the time iswhere is the position of the mass center at the time .
The center of the cavitator in at the time isSimilarly, is the mass center at the time .
According to the principle of independence of the supercavity sections expansion, the lift and lateral force on the cavitator at the time leads to the offset of supercavity center along and axis at the time . This offset at the time can be represented by in and in or . Furthermore, (10) is modified by considering the supercavity axis deviation induced by gravity. The center position of the corresponding supercavity section at the time in isConsidering (8), the last item of (11) at the right-hand side is
When sailing at the fixed-depth and fixed-direction mode, the angle of attack and angle of sideslip can be regarded as small quantities. This is also confirmed by the numerical simulation results in Section 6.3 (simulation results show the magnitude of and are less than 1.5°). Therefore, the small quantity simplifications are made, that is, , , , , and . Equation (12) is further rewritten as
Accordingly, a vector pointing from the vehicle center to the supercavity center at the tail is represented by in and at the time . This vector in is
This vector describes the offset between the vehicle center and the supercavity center at the tail. The rolling control is necessary for an underwater vehicle, and its objective is to eliminate the angle of rolling. For supercavitating vehicles, the control of rolling is seldom influenced by the time-delay effect. The angle of rolling can be regarded as approximately zero [23, 25]. Therefore, , , and in (14) are further simplified as
The offset between the supercavity center and vehicle center at the tail is then obtained once the delayed time is confirmed. implies the cross section of supercavity propagates to the tail of vehicle exactly, and the corresponding delayed time is of interest in this paper. Because of the complexity of (15), simplifications are imperative to calculate . As mentioned before, when the vehicle sails at the fixed-depth and fixed-direction mode, , , , , and are small quantities. In addition, , , , , are all far smaller than . Considering the product of two or more small quantities is towards zero, (15) is simplified aswhere is the vehicle length excluding the engine nozzle. The delayed time is calculated as
Regarding (15) and (16), the vertical and horizontal components of the vector from the vehicle center to the supercavity center in , named as and , respectively, are calculated by considering (19). , , and are written as
The offsets , , and of the axis of supercavity induced by the lift and lateral force of cavitator and gravitational force are expressed as [3, 4, 28, 29]
5. Spatial Kinetics Equations
Now that the hydrodynamic forces acting on the supercavitating vehicle and the relative position are defined, the spatial kinetics equations are established in this section. The force balances of the vehicle at the axis in the vertical plane and horizontal plane are
The moment balances at the axes , , and arewhere
These relationships are obtained by matrix transformation at different coordinates [8]. The coordinates of the mass center in are
The kinetics model described by (24)–(29) is solved by inputting rudder angles , , and . The detailed process is shown in Figure 5.
6. Numerical Simulations
6.1. Initial Conditions
Mathematical simulations are performed based on the spatial kinetics model and the hydrodynamic forces. The detailed simulation conditions are listed in Table 1.
6.2. Control Strategy
For the control system of supercavitating vehicles, the closed-loop PD algorithm is utilized to regulate the sailing direction, sailing depth, and rolling motion. The pivot disk cavitator, upper-vertical rudder, and lower-vertical rudder are the actuators of the control system of supercavitating vehicles.
It should be emphasized that the motion of supercavitating vehicles is dynamically steady [8, 20]; therefore the fluctuated feedback signal (sailing depth and the angle of yaw at the mass center) results in the oscillated control forces. This effect increases the possibility of the amplified oscillation or even instability. For supercavitating vehicles, its instantaneous velocity center locates slightly behind the center of cavitator and the corresponding oscillation is approximately zero [30]. Therefore, the sailing parameters at the location of the instantaneous velocity center are taken as the characteristic parameters of the vehicle. The control effect is negligible when the vehicle flights are in a dynamically steady condition.
The control of the sailing direction can be realized by average rudder angle of vertical rudders:
The sailing depth is controlled by the inclination angle of cavitator:
The control of rolling motion is adjusted by operating the differential angle between the upper-vertical rudder and lower-vertical rudder:
Combining (30) and (32), the rudder angle of the upper-vertical rudder and the rudder angle of the lower-vertical rudder are
The motion of cavitator and vertical rudders are driven by direct-current motors. In practice, the response characteristic of rudder angle to the control command is inertial and the time constant is approximately 10 ms [31].
6.3. Numerical Results
The obtained results are shown in Figures 6–11. Here, the length scale is normalized by the total length of the vehicle and the time is normalized by the time delay of the supercavity at the tail . Figure 6 presents the speed of the vehicle in terms of the nondimensional time. The vehicle is launched at an initial speed of 70 m/s and is accelerated (boosted by a solid propellant rocket) to the speed of 100 m/s at approximately nondimensional time of 8. At this condition, the formed supercavity is large enough to envelop the whole vehicle. The upper-stream dynamic pressure is able to meet the operational requirement of the water-ramjet engine. The vehicle then proceeds to the cruise mode with an approximately constant speed of 100 m/s.
Figures 7 and 8 depict attitude (angle of pitch and yaw) and hydrodynamic (angle of attack and sideslip) angles during the sailing, respectively. It is demonstrated that both attitude and hydrodynamic angles oscillate and amplify gradually with the nondimensional time. The amplitude of the oscillation reaches roughly to the constant at the nondimensional time of 20. Due to the continuous effect of gravitational force, the oscillation in the vertical plane (angle of pitch and attack) occurs earlier and reaches to the steady state faster than in the horizontal plane (angle of yaw and sideslip). As shown in Figure 8, the oscillation periods of angles of attack and sideslip are approximately constant and are defined as and , respectively. Considering ≈ 26.9 ms for the vehicle, the oscillation periods are ~85 ms and the frequency is equal to about 11.8 Hz. Furthermore, a phase difference (between the two angles) of about is obtained at the ranges of nondimensional times 20–40 and 60–80. The same phenomenon is also found by comparing angle of pitch and angle of yaw as shown in Figure 7. For the nondimensional time 40–60, strong control is activated because of the significant change of attitude angles. This influences the oscillation characteristic of the vehicle correspondingly.
To show the existence of the relative movement between the supercavity and vehicle, the vector at the tail of vehicle from the supercavity center to the tail center is exhibited in Figure 9. This figure shows the relative position between the supercavity center and vehicle center at the tail. The horizontal and vertical axes represent the nondimensional relative positions and , respectively. A continuous revolution of the tail center of the vehicle with respect to the axis of the supercavity is clearly observed. Given that the axis of the vehicle intersects with the axis of the supercavity at the location near cavitator, this specific revolution leads to a conic-like movement and the semicone angle remains to be a constant.
The nondimensional distance between the centers of the vehicle and supercavity at the tail is the reflection of their intersection angle. Figure 10 shows this distance as a function of the nondimensional time. The magnitude changes between 0.01 and 0.015 during most of the sailing. The wetted length of the vehicle is associated with the hydrodynamic force and moment acting on the vehicle and is determined by the offset between the vehicle and supercavity. This wetted length is shown in Figure 11. The averaged wetted length is about half of the length of vehicle. The characteristic of planing force determines the motion mode of the vehicle and its dynamic process also shows the relative position between supercavity and itself. The net planing force is normalized by the half density of water, the maximum cross section area, and the square of real-time sailing speed. The dynamic process of coefficient of net planing force is depicted in Figure 12. As shown in Figures 10, 11, and 12, the zero value of the distance, the wetted length, and net planing force are seldom observed. This demonstrates that the oscillation of supercavitating vehicles does not exist in a single plane; otherwise the zero nondimensional distance and wetted length will occur periodically. This confirms that the “tail-slap” phenomenon exists at the vertical and horizontal planes simultaneously.
In summary, an intrinsic oscillation exists during the motion of supercavitating vehicles. The amplitude and frequency of attitude and hydrodynamic angles are approximately the same, but a phase difference of about is observed as shown in Figures 7 and 8. The axis of the vehicle is revolving around the supercavity axis. This results in a conic-like oscillation.
7. Experimental Investigation
To validate the oscillation characteristic as introduced in Section 6.3, a lakebed experiment is carried out. The experimental vehicle has the same geometric, weighting parameter as listed in Table 1 and the same control strategy described in Section 6.2. In the lakebed experiments, the supercavitating vehicle was launched from a stockade launch tube, which is mounted on a pontoon catamaran as shown in Figures 13 and 14 and exhibits the installation of the vehicle in the launch tube. An inertial measurement unit was used to measure the translational acceleration and angular velocity in , , and -axis. A PC104 was used as the central control unit to perform the sailing control and testing data storage.
A smoothing filter is adopted to remove the noise of the testing data. The dynamic process of the angle of yaw and angle of pitch is obtained by integrating the angular velocity in the -axis and -axis as shown in Figure 15. Again, the testing time is normalized by .
Figure 15 demonstrates that the simultaneous oscillation in horizontal and vertical plane exists in the motion of supercavitating vehicles. The angle of pitch and angle of yaw exhibit the approximately periodic oscillations. The period is about 76.9 ms if ≈ 26.9 ms. By comparing Figures 7 and 15, similar oscillation amplitude and frequency are found for the angle of pitch and yaw between mathematical simulations and the lakebed experiment. The oscillation frequency from the simulation (shown in Figures 6 and 7, 11.8 Hz) and the experiment (shown in Figure 14, 13 Hz) is close. Also, the same finding of an approximately steady phase difference (about a quarter of a period) between the angle of pitch and angle of yaw is found. This validates the suitability of the developed mathematical model in predicting the kinetics motion of supercavitating vehicles and concludes that a conic-like oscillation is an intrinsic characteristic.
The track of the tail center in the translational coordinate system is obtained by multiplying the tangential components of the time-varying angles by . The result at nondimensional times 28 to 38 is shown in Figure 16. The horizontal axis and vertical axis denote -position and -position of the tail center of vehicle (not including the nozzle), respectively. It is found that the tail-slap motion at the tail center of the vehicle occurs not only in the vertical plane but also in the horizontal plane simultaneously. This further validates the conic-like oscillation motion of supercavitating vehicles as presented in Section 6.2.
8. Conclusions
In this work, a spatial kinetics model including the coupled motion of pitch, yaw, and rolling for supercavitating vehicles was developed. The coupled relative position between the supercavity and the vehicle was described to compute the planing force. Key findings are as follows:(i)An intrinsic conic-like oscillation during the sailing is investigated with the developed kinetics motion model and validated by a lakebed experiment; and the oscillation frequencies from the mathematical simulations and the experiment are close.(ii)At the steady state sailing section, a phase difference of is observed between the motions in the vertical and horizontal planes.(iii)The tail-slap motion exists not only in the vertical plane, but also in the horizontal plane. This tail-slap motion in the two orthogonal planes and the phase difference result in a spatial conic-like oscillation.
Nomenclature
| : | Drag force of vehicle |
| : | Slope of lift force coef. of cavitator |
| : | Slope of lateral force coef. of cavitator |
| : | Slope of lateral force coef. of lower-vertical rudder |
| : | Slope of lateral force coef. of hull |
| : | Slope of lateral force coef. of upper-vertical rudder |
| : | Inclination angle of cavitator |
| : | Damping force on hull |
| : | Position force on hull |
| : | Forces on the cavitator along |
| : | Oscillation frequency |
| : | Gravitational force |
| : | Gravitational acceleration |
| : | Offset between mass center and buoyancy center |
| : | Rotational inertia about |
| : | Rotational inertia about |
| : | Rotational inertia about |
| : | Correction coefficient |
| : | Length of cavity |
| : | Distance from mass center to rudder |
| : | Distance from mass center to cavitator |
| : | Length of vehicle |
| : | Distance from mass center to tail |
| : | Arm of planing force about mass center |
| : | Pitch moment on cavitator |
| : | Rolling moment on upper-vertical rudder |
| : | Yaw moment on upper-vertical rudder |
| : | Pitch moment on hull |
| : | Rolling moment on hull |
| : | Yaw moment on hull |
| : | Rolling moment on lower-vertical rudder |
| : | Yaw moment on lower-vertical rudder |
| : | Mass of vehicle |
| : | Radius of cavitator |
| : | Radius of supercavity at tail |
| : | Maximum radius of vehicle |
| : | Force arm of upper rudder about |
| : | Force arm of lower rudder about |
| : | -position of supercavity center at tail in |
| : | Maximum cross-sectional area of vehicle |
| : | Propulsion |
| : | Time |
| : | Sailing speed |
| : | Component of velocity along |
| : | Lift force on cavitator |
| : | Lift force on cavitator at time |
| : | Lift force on hull |
| : | Desired sailing depth |
| : | Offset of supercavity axis along caused by cavitator |
| : | Offset of supercavity axis caused by gravity |
| : | -position of tail center of vehicle in |
| : | Lateral force on cavitator |
| : | Lateral force on cavitator at time |
| : | Offset of supercavity axis along caused by cavitator |
| : | Lateral force on upper-vertical rudder |
| : | Lateral force on hull |
| : | Lateral force on lower-vertical rudder |
| : | -position of tail center of vehicle in |
| : | Angle of attack |
| : | Angle of sideslip |
| : | Average angle of vertical rudders |
| : | Differential angle of vertical rudders |
| : | Angle of upper-vertical rudder |
| : | Angle of lower-vertical rudder |
| : | Trajectory inclination angle |
| : | Angle of pitch |
| : | Medium density |
| : | Time delay |
| : | Angle of rolling |
| : | Trajectory deflection angle |
| : | Angle of yaw |
| : | Desired angle of yaw |
| : | Rolling rate |
| : | Yaw rate. |
| : | Translational coordinate system |
| : | Body-fixed coordinate system |
| : | Earth-fixed coordinate system |
| : | Time . |
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The project is supported by the National Natural Science Foundation of China (NSFC) (Grant nos. 51579209, 51679202, and 51409215). The authors acknowledge them for their sponsorship.