Abstract
This paper presents a novel design for a bearingless axial flow blood pump based on the magnetic-hydrodynamic double levitated concept. In the axial direction, the magnetic levitation system consisted of two pairs of permanent magnet rings offsets the force of fluid. The hydrodynamic shell mounted on the impeller rotor is designed for generating dynamic pressure, which can balance the radial force like gravity when the blood pump is working. Because of the unsteady force and torque acting on the rotor and the passive suspension, the position of the rotor is not steady. The suspension force, stiffness, and torque of the rotor are calculated by the theoretical method and finite element method. Then, the dynamics of the rotor are analyzed. Arrangements of Hall-effect sensors with the corresponding data acquisition system which can measure the axial displacement of the rotor are explained. The sensorless drive control system for the blood pump is described too. With a prototype pump, an external circulation experiment system is built and then the axial and radial displacements of the rotor are measured by using Hall-effect sensors and the laser vibrometer under different working conditions.
1. Introduction
The high speed of the mechanical bearing in the blood pump could cause hemolysis, thrombus, abrasion, and other sequelae [1, 2]. These disadvantages of mechanical bearings limit long time use of a blood pump. In the third generation of the blood pump, noncontact bearings (bearingless motor) are utilized for suspension of the rotor [3]. For compact structure, no mechanical contact, no mechanical wear, no friction heat, and other superiority, noncontact bearings are applied in many fields such as machine tool spindles, molecular pumps [4], high-speed motors [5], and blood pumps [6]. The magnetic force, electromagnetic force, and hydrodynamic (HD) force are often applied in the suspension structure of a blood pump. In the Terumo DuraHeart [7] and INCOR pump [8], electromagnetic bearings are employed to axial balance and magnetic bearings are applied for rotors’ radial suspension. Hydrodynamic bearings are intended in the CorAide pump [9] and LVAS [10] for the suspension of impeller rotors. The hydrodynamic bearing together with the magnetic bearing solves the five degrees of freedom suspension in HVAD [11]. The hydrodynamic bearing combines with the electromagnetic bearing, and MVAD realized the unsupported rotation [12].
Generally, magnetic and hydrodynamic suspension is passive and uncontrollable but the electromagnetic suspension is active and controllable. For the additional control system, the volume of electromagnetic suspension is much larger than that of others. Compared with traditional bearings, there is no mechanical fixing structure such as mechanical bearing and driveshaft in noncontact bearings, and the space is saved, but the position of the suspension rotor may be more changeable. For assessing the performance of a blood pump, the movement state of the blood pump’s rotor should be analyzed and measured. This paper proposes the design of an axial blood pump with a magnetic-hydrodynamic double levitated (MHDL) impeller rotor. Bearings in the pump are all passive control which makes the structure of the pump more simplified. The structure of noncontact bearings is modeled, and then the bearing force, stiffness, and torque of the rotor are calculated by the theory method and analyzed by the three-dimensional finite element method (FEM). With a prototype pump and sensorless drive control system, the external circulation system of the blood pump is built. Special arrangements of Hall-effect sensors with the data acquisition system (DAS) are assembled to measure the axial and radial displacements of the rotor. The result indicates the rotor can be suspended well at rated speed, and the bearing system is stable enough for working and shock.
2. Mechanical Design
The schematic and protype of the MHDL axial blood pump are shown in Figure 1 [13]. First, the blood flow is straightened by the front guide vane and then pressured by a rotating impeller, at last outflow along the rear guide vane. The motor consists of a conductor, stator core, coil winding, and permanent magnet (PM) of the rotor. PM ring I-2, PM ring II-1, magnetic isolation ring I, magnetic isolation ring II, end cap, PM of the rotor, and HD shell are all fixed on the impeller, which make up the rotor of the motor and can be driven to rotate together. Four PM rings are magnetized axial N50-grade neodymium-iron-boron (NdFeB), and the PM of the rotor is magnetized radial NdFeB. The magnetic pole direction of four PM rings I-1, I-2, II-1, and II-2 in Figure 1 are S-N, N-S, S-N, and S-N, respectively.
(a)
(b)
Magnetic forces between PM rings I-1 and I-2 and II-1 and II-2 make up the axial bearing force. Electromagnetic force of the stator also can play a little part for the axial bearing when the rotor leaves from a special position. As shown in Figure 2, four spiral grooves separate the HD shell to four sheets whose diameter is ranged from small to large in order to form wedge gaps with the inner surface of the pump housing. The average width of the wedge gaps is about 0.4 mm. The blood in wedge gaps generates HD pressure which can suspend the rotor in the radial direction. The injury of blood in the domain HD pressure can be lightened for the exchange of blood through spiral grooves. The shaft, whose diameter is 0.36 mm less than the diameter of center holes of the rotor, can limit the initial position and posture of the rotor. Geometric parameter values of noncontact bearings are presented in Table 1. There is not any control system for bearings for both the magnetic and HD bearings are used passively, so the structure of the blood pump is simplified and compact.
(a)
(b)
(c)
(d)
3. Torque and Force of Bearings
The position of the rotor can affect the bearing force and torque which are generated from PM, electromagnet, and HD. Depending on the geometric parameters of parts of the blood pump, the axial space, radial space, and deflection space are 0–1.2 mm, −0.18–0.18 mm, and −0.0136–0.136, respectively.
3.1. Magnetic Bearing
Two pairs of PM rings are utilized for axial stabilization. The relative position of two PM rings can be modeled as presented in Figure 2. Establish coordinate systems o-xyz and o’-x’y’z’ with the center and the center axis of each PM ring as the origin and z-axis, respectively. Positional relationship between two origins of coordinate can be expressed by oo’ = (ex, ey, ez). The angles between axes of two coordinate systems are listed in Table 2. The axial thicknesses, inner diameters, outer diameters, and magnetization of two PM rings are h1 and h2, Rd1 and Rd2, RD1 and RD2, M and M’, respectively. Two end faces of one PM ring are denoted as numbers 1 and 2, and those of the other PM rings are denoted as numbers 3 and 4. Because of the shape of PM rings, it is convenient to apply cylindrical coordinates for position representation of the element on the PM rings.
Depending on the theory of the equivalent magnetic charge, the PM is regarded as the superposition of many discrete magnetic units, and the interaction between the two PM rings is also regarded as the superposition of the interaction of the two respective magnetic units. [14] The magnetic Coulomb’s theorem is employed to calculate the magnetic force between units. The interaction of the two magnetic rings is regarded as the sum of the actions of the end faces 1 and 3, 1 and 4, 2 and 3, and 2 and 4, for only the torus of the end face has magnetic charges. Taking the interaction between the end faces 1 and 3 as an example: P is an element on end face 1, the location in coordinate systems o-xyz is (r1cosα, r1sinα, z1), and the magnetic charge q1 can be calculated by the following:where Br1 is the area density of magnetic charge and dAp1 is the area of P. The magnetic charge of an element Q on the end face 3 q3 can be obtained for the same principle as q1:
According to the magnetic Coulomb’s theorem and coordinate transformation formula, the magnetic force between faces 1 and 3 F13 can be expressed by the following:
F13, x, F13, y, and F13, z are the components of F13 in the direction of axis. µ0 is the vacuum permeability. A13, B13, and C13 are
The magnetic force between two PM rings FPM is the sum of forces between face 1 and face 3 F13, face 1 and face 4 F14, face 2 and face 3 F23, and face 2 and face 4 F24:
The force between elements P and Q is d FPQ, the distance vector is PQ, the torque of element Q under the action of element P is d MPQ = d FPQ × PQ, and using the same integration area and superposition as magnet force calculation, the torque between two PM rings can be computed.
3.2. Electromagnetic Driven System
There are six windings, which are three-phase six-state in one turn in the stator of the blood pump’s driven system. The electromagnetic force and torque of the stator can drive the PM assembled in the impeller. The force and torque of PM in the electromagnetic field generated by winding can be calculated by the following [15]:where Fpx, Fpy, and Fpz are the components of force act on PM; r1 and r2 are the inner and outer diameter of PM, respectively; l is the length of PM; Mi, Mj, and Mk are the components of magnetization; ΔBx = Bx∣o1 − Bx∣o2, ΔBy = By∣o1 − By∣o2, ΔBz = Bz∣o1 − Bz∣o2, o1, and o2 are the center point of the interface between two magnetic poles; Bx, Bx, and Bx are the components of magnetic intensity generated by windings which can be calculated according to the Biot–Savart law; dMp is the torque acting on an element of PM; and θ is the rotation angle of PM.
3.3. HD Bearing
The HD bearing is applied for radial support of the rotor in our design of the blood pump. The cross section of HD gap between the HD shell and inner of pump housing is shown in Figure 3. The theory of fluid lubrication is applied for the analysis of the HD bearing. Based on the two-dimensional Reynolds equation for HD lubrication and Ocvirk solution for infinitely short bearing assumption, the local oil film pressure p acting on the HD shell is expressed by the following [16]:where z and B are the axial position and axial length of HD gap, respectively; η is the dynamic viscosity of fluid which passes through the HD gap; is the tangential speed of the element of HD shell’s blades where the HD pressure acts on; Cr is the radial clearance, given by (D − Di)/2, where D is the bearing diameter and Di is the journal diameter; r is the bearing radius; and ε is the eccentricity ratio, given by e/Cr, where e is the absolute bearing eccentricity.
The bearing force dWi on the ith blade with axial length of dz can be obtained by the integration of equation (8):where dWix and dWiy are the components of dWi. In the interval of [φi1, φi2], the pressure of the oil film is continuously positive. cb is the bearing coefficient related to the structure of HD bearing. The total bearing force W is as follows:where Wx and Wy are calculated by the following:
3.4. FEM of the Bearing System
In order to verify the theoretical analysis of the bearing system, FEM is used for force and torque calculation of the magnetic bearing, electromagnetic driven system, and HD bearing. Maxwell 3-D analysis is employed in the simulation of the magnetic bearing and the electromagnetic driven system. Parametric modelling is applied for the position change of PM in a certain range, and torque and force are set to be output options. In the simulation of the coupling of PM ring II-1 with II-2, the magnetic flux density at radial 9 mm and different axial positions is researched for looking for a position where the gradient is largest so that the position change of the rotor can be sensed more sensitively by Hall-effect sensors. In the simulation, the Br of PM is set to be 1.43 T. The stator is made of pieces of silicon steel. The number of winding is 75 per stator pole with the copper wine of 0.32 mm diameter.
Computational fluid dynamics (CFD) are performed in the analysis of the HD bearing. The axial and radial forces of the rotor are computed in ANSYS Fluent with the rotation speed in the range of 5000–10000 rpm and rotor eccentricity in the range of 0–0.18 mm. Pressure inlet and pressure outlet are set up in boundary conditions for the simulation of demanded conditions.
4. Results of Force and Torque
The axial forces of fluid flow acting on the impeller are evaluated to be 1.90 N, 1.98 N, 2.06 N, and 2.14 N when the speeds are 6000 rpm, 7000 rpm, 8000 rpm, and 9000 rpm, respectively, according to the CFD and hydraulic experiment. For the external factors in actual situation and the speed of the rotor is detected and controlled all the time, there is often a little disturbance in the force of fluid flow.
The reluctance forces on the rotor response to the axial displacement are displayed in Figure 4, and the balance position of the rotor is 0.22 mm when the rotor is stopped. The magnetic bearings can supply axial −0.73–3.53 N bearing force when the axial position of the rotor varies from 0.1 to 0.9 mm. The axial balance positions of the rotor are 0.559 mm, 0.575 mm, and 0.590 mm when the speeds of the rotor are 6000 rpm, 7000 rpm, and 8000 rpm, respectively, based on the axial force of the fluid flow and magnetic bearings.
As illustrated in Figures 5–8, the radial displacement and deflection have little effect on the axial and radial reluctance forces and yawed torques especially when the rotor at an axial working balance position.
As shown in Figure 9, the increase in rotation speed and eccentricity of the rotor can result in a larger HD bearing force. The rotor begins to be suspended at a rotor speed of 5000 rpm (the HD bearing force is larger than the gravity of the rotor when the radial eccentricity reaches the maximum value). The eccentricity of the rotor is 0.05 mm when the rotor speed achieves rated speed. The HD bearing force can reach five times the gravity of the rotor when the eccentricity of the rotor is 0.12 mm at rated speed and so the HD bearing is stable enough for working. The eccentric direction of the rotor is near 225° for the x component and y component of the HD bearing forces which are approximately equal as shown in Figure 10.
The radial movement of the rotor is mainly determined by the balance of gravity and HD bearing force. For the relationship among the directions of the radial eccentricity of the rotor, HD bearing force, and gravity of the rotor, the direction of radial movement of the rotor, which is assumed to be the axis of the coordinate system for simplifying the dynamic analysis, is chiefly 45° to the direction of gravity. The electromagnetic force changes with rotation angle of the rotor at different speeds are shown in Figure 11. The increase in rotation speed leads to the decrease in the radial eccentricity of the rotor, and then the radial electromagnetic force of the stator acting on the rotor decreases.
The comparison of simulation results with theoretical results is shown in Figures 12 and 13; the magnetic force between PM ring II-1and II-2 and the HD bearing force of the rotor with the HD shell are selected as a representative. The difference in treatment of air gap permeability and meshing between the simulation and theoretical calculation leads to the different results in the calculation of magnetic force, but the deviation is smaller than 6%. In the process of HD bearing force calculation, the distinction of solving the Navier–Stokes equation causes the simulation results are about 5% larger than theoretical results. The comparison shows that the theoretical study on the magnetic and HD bearing force is reliable, and data obtained can be used for the following research.
4.1. Dynamic Model of the Rotor
The schematics of the different forces acting on the rotor of the pump and corresponding coordinate systems are shown in Figure 14. The magnetic forces exerted on the rotor by PM rings I-1 and II-2 are denoted by Fex1, Fey1, and Fez1 and Fex2, Fey2, and Fez2, respectively. The electromagnetic forces and torques generated from stator windings are represented by Fex3, Fey3, and Fez3 and Mx, My, and Mz, respectively. The distance between the centers of PM rings I-1 and PM of the rotor is l1. The distance between the centers of PM rings II-2 and PM of the rotor is l2. The HD force acting upon the HD shell is denoted by Fdx and Fdy. The fluid force exerted on the impeller is represented by Fldx, Fldy, and Ff. The center of mass and gravity of the rotor and the center of mass and gravity of the PM of the rotor are designed to be coincident at cm. The x-axis dynamics of the rotor can be expressed as follows:where m is the mass of the rotor and b is the viscous friction coefficient which is assumed to be constant in all directions. Ix and Iz are the mass moments of inertia about θx and θz axis, respectively. The Fdy on the axial element of the HD shell is dFdy, and lz is the axial distance between the element with cm. The angular velocity of the rotor is ωm. The blades of the impeller are axisymmetric about x-axis, and Fldx and Fldy are the drag forces on the blades of the impeller, so the Fldx and Fldy can only have z-axis torque effect which is offset by driven torque Mz [17].
(a)
(b)
(c)
The y-axis dynamics of the rotor can be written as follows:where is the acceleration of gravity, dFdx is the Fdx acting on axial element of the HD shell, and Iy is the mass moments of inertia about θy axis.
When the blood pump is working, the forces and torques of magnetic and HD bearings acted upon the rotor are static if the position of the rotor is stable, but the electromagnetic force and torque are changed periodically and so can be regarded as external excitations. As shown in Figures 8 and 10, the radial force of magnetic bearings and HD bearings can be treated as linearly proportional to the rotor’s radial displacement when the rotor is near the radial working balance position. The scale factor between the radial force of the HD bearing and the radial movement of the rotor is changed with the speed of the rotor as described in Figure 10. The natural frequency of the rotor in radial displacement is evaluated to be 411 Hz, 506 Hz, 559 Hz, and 588 Hz when the speed of the rotor is 6000 rpm, 7000 rpm, 8000 rpm, and 9000 rpm, respectively. The natural frequency of the rotor is far from the frequency of rotation speed.
All the bearings including the stator can provide torque to prevent rotor deflection around z-axis. There is not any external yaw moment excitation, so the rotor can be stable without any deflection.
The y-axis dynamics of the rotor are given as follows:
Ff, which is related to the rotating speed of the impeller, is axial force of fluid flow and can be treated as external disturbance. Fez1, Fez2, and Fez3 have a relationship with the axial position of the rotor, so the axial equilibrium position of the rotor is changed with Ff. The rotation of the rotor is controlled under the balance of Mz and drag torque caused by Fldx and Fldy.
As shown in Figure 5, the axial force of magnetic bearings is linearly proportional to the rotor’s axial displacement. The stiffness of axial force of magnetic bearings can be approximated from the slope of the line shown in Figure 5. The natural frequency of the rotor in radial displacement is calculated to be 569 Hz. The frequency of external disturbance Ff is lower than the speed of the rotor in actual situation, so the resonance will not happen in the axial of the rotor system.
4.2. Control System
The overall block diagram of the sensorless driven control system for the bearingless motor is shown in Figure 15. The six-arm full-bridge drive circuit is designed with STM32F103 as the main controller chip based DSP. There is no sensor in the control system for the position and the speed of the rotor is obtained by zero-crossing detection. The closed-loop control mode based on the PI controller is applied for the speed control of the rotor. The voltage, current, and flow of the pump are monitored during working. The drive control system responds fast and has a less than 1% steady error [18]. The pump with a driven control system can fulfil the requirements of flow and differential pressure of human blood circulation.
4.3. Rotor Position Sensing
It is hard to install a position sensor on the rotor directly for the rotor which is rotating and suspended in fluid when working. The noncontact sensor like a capacitive sensor, ultrasonic sensor, eddy current sensor, Hall-effect sensor, and laser sensor is needed for the position detection of the rotor. The capacitive sensor is not reliable here for the velocity of blood can affect the dielectric constant. The irregular structure of the pump can generate echo noise for the ultrasonic sensor. The eddy current sensor requires extra system which would increase the volume of the pump. In order to observe the rotor, pump housing is made of transparent material and so the radial displacement of the rotor can be sensed by a laser. But, axial end faces of the rotor are blocked by the front and rear guide vane. For the magnetic bearing in the pump, the Hall-effect sensor is convenient to sense the change of magnetic flux density caused by the axial movement of the rotor [19]. In our research, Allegro linear Hall-effect sensors A1324 with a sensitivity of 5 mV/Gauss and an internal bandwidth of 17 kHz are employed for radial displacement measurement of the rotor.
In order to measure the axial position of the rotor, four linear Hall-effect sensors are fixed evenly around the radial surface of the outside of pump housing as shown in Figure 16. The sensing faces of sensors are parallel to the axis. According to FEM, the gradient of magnetic flux density can reach biggest value when the axial distance between the sensors with PM ring II-1 is about 4.5 mm. The magnetic flux density can be regarded as linear when the movement of the PM ring is small. The linear Hall-effect sensor can provide linear output in response to a ramping applied magnetic field. The output voltage of sensors in response to the displacement of the PM ring can be modeled as follows [20]:where ΔUx1, ΔUx2, ΔUy1, andΔUy2 are the change of output voltage of the sensors x1, x2, y1, and y2, respectively; ΔBx1, ΔBx2, ΔBy1, andΔBy2 are the change of flux density sensed by four sensors; kE is the constant which is related to the sensor; c0 and c1 are constants; (x, y, z) is the location of the PM ring before movement; and (Δx, Δy, Δz) is the displacement of the PM ring. Equation (16) can be transformed to be the following:
(a)
(b)
A low-pass filter is designed for the sampling of signals of four sensors in order to reduce the high-frequency noise. The cutoff frequency of the filter is set to be 4 kHz according to the rotating speed of the rotor. The scheme of axial displacement measurement is low-cost, compact-in-size, and has reasonably good accuracy.
5. Experimental Results and Discussion
A prototype with an external circulation experiment system was built, and the performance and operation of systems mentioned above were investigated too. For the transparent material of pump housing, the Julight laser vibrometer VSM1000 based on the interferometric configuration can be employed to measure the radial displacement of the rotor. The experimental setup for measuring the axial displacement of the rotor is shown in Figure 17. The experimental setup for measuring the radial displacement of the rotor is shown in Figure 18. The prototype was fixed on the platform, and the Hall-effect sensors were stuck to the pump housing in order to reduce the error of measurement. The signals of Hall-effect sensors and the laser vibrometer were acquired and analyzed by the data acquisition (DAQ) system (the sample frequency is set to be 2 kHz–4 kHz in the experiment).
The axial oscillation of the rotor when the pump was running at 6000 and 8000 rpm is shown in Figure 19. The error of axial displacement is about 0.03 mm when the rotor at rate speed, which means the axial space of the gap between the rotor and rear guide vane is in the range of 0.58–0.605 mm. The increase in rotation speed of the rotor causes larger axial movement of the rotor, axial disturbance force, and frequency of axial disturbance force. Compared with electromagnetic actuators applied in other blood pumps, the accuracy and the axial gap between the rotor and guide vane of axial passive control of the PM bearing applied in this paper is a little worse, but there is no additional component and corresponding control system in the PM bearing. The passive control of the PM bearing has a compact structure and is good for long-term use.
(a)
(b)
The y-axis performance of the rotor at the speed of 4800 rpm and 5000 rpm is shown in Figure 20. The rotor begins to be suspended at 4800 rpm and is fully suspended at 5000 rpm. The speed of the rotor should reach 5000 rpm as soon as possible during the startup of the blood pump for the rotor has contact with shaft lower than this speed.
(a)
(b)
The x-axis and y-axis performances of the rotor at speed of 6000 rpm and 8000 rpm are shown in Figures 21 and 22, respectively. The x-axis and y-axis eccentricities of the rotor at the speed of 6000 rpm are all near −0.05 mm. Both the x-axis and y-axis eccentricities of the rotor at the speed of 8000 rpm are near −0.03 mm. The vibration amplitude and frequency of the rotor in the x-axis and y-axis are similar. For the precision of HD shell, impeller, pump housing, stator, and so on, there are some differences between the irregular waveforms of oscillation of the rotor in x-axis and y-axis. The arrangement of laser sensors can affect the waveforms too. The vibration amplitude of the rotor in the radial direction is about 0.015 mm at the speed of 6000 rpm and 0.011 mm at the speed of 8000 rpm.
(a)
(b)
As shown in Figures 23 and 24, fast Fourier transform (FFT) is employed in the analysis of the magnitude spectrum of the radial oscillations. The domain frequency is about 100 Hz at 6000 rpm and 140 Hz at 8000 rpm corresponding to the operation speed of the rotor. The eccentric distance of the rotor is much larger than the amplitude of displacement, so the FFT magnitude is large at zero frequency which represents the DC component of the signal. The vibration of radial displacement of the rotor is the response of the electromagnetic force acting on the rotor which is shown in Figure 11. The frequency of system input (electromagnetic force acting on the rotor) increased with increase in the rotating speed, and then the frequency of system response would increase too. So, the domain frequency of the radial displacement of the rotor increases when the rotating speed increases from 6000 rpm to 8000 rpm. In the radial direction, the natural frequency of the rotor is about four times as big as the domain frequency, so there is no resonance in the radial displacement of the rotor when working which ensures the stability of the radial bearing system. The HD bearing has a good performance in radial suspension for small eccentricity, tiny error, large suspension force, and high natural frequency. But, the high stress in the HD gap may injure the red blood cells, and spiral grooves on the HD shell are designed in this paper to reduce the injury of red blood cells.
6. Conclusions
A novel design of a bearingless rotor for the axial blood pump has been proposed in this paper. Two pairs of PM rings are utilized for the axial suspension of the rotor, and the HD shell is employed for the radial suspension of the rotor. Both the axial suspension and radial suspension are passive control. The dynamic model of the rotor and the analytical model for the axial and radial suspension forces were demonstrated. For the special environment of the rotor, four linear Hall-effect sensors were fixed up to measure the axial position of the rotor. A prototype with an external circulation system was assembled to guide the design and optimization of the blood pump. The experiment results indicated that the rotor can be successfully suspended when the rotation speed of the rotor is above 5000 rpm. There are oscillations in both the axial and radial displacements of the rotor because of the passive control of noncontact bearings and external disturbance. The errors of axial and radial placement are about 0.03 mm and 0.02 mm, respectively, when the rotor is at rated speed. The domain frequency of the rotor which corresponds to operation speed is much smaller than the natural frequency.
Data Availability
The displacement accuracy of rotors of different bearings in blood pumps is used to support the findings of this study and is available at 10.1109/TIE.2017.2762626 and 10.1109/TIE.2009.2017095. These prior studies (and datasets) are cited at relevant places within the text as references [17, 20]. The displacement data of the rotor of the blood pump researched in this paper 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 with respect to the research, authorship, and/or publication of this article.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (grant no. 31670999), Natural Science Foundation of Hunan Province, China (grant no. 2019JJ40364), and Zhejiang Provincial Top Key Discipline of Mechanical Engineering (grant no. GK170201201003/002).