2.2. Modeling of the Bearingless Motor
The magnetic poles of the motor rotor were fixed on the rotor backplate. Therefore, the motor was a non-salient pole structure, and the magnetic poles were magnetized according to the Halbach array arrangement. The specific composition of the magnetic suspension system and permanent magnet bearing are shown in
Figure 3a,b.
The two concentric, permanent magnet rings embedded in the upper-end face of the impeller and the two concentric, permanent magnet rings embedded in the upper casing cooperated with each other, and the permanent magnet bearing generated a magnetic pull force in the
z axial direction. Therefore, the component of the magnetic force generated by the motor stator and the permanent magnet rings in the
z axial direction jointly determined the axial position of the impeller. The
x,
y axial translation and
ϕx,
ϕy rotation of the impeller were passively constrained by the magnetic coupling between the permanent magnet bearing and the bearingless motor. Therefore, the translation and rotation of the impeller in
x and
y axis direction were limited. Let the upper and lower air gaps of the rotor be represented as
g1 and
g0, respectively. When the rotor is translated
z in the positive axial direction due to disturbance, the upper and lower air gaps becomes
g1 −
z and
g0 +
z. Therefore, in order to maintain the stable suspension and rotation of the impeller at the balance position, it is crucial to regulate the motor stator current. The specifications of the bearingless motor and permanent magnet bearing are shown in
Table 1.
The bearingless motor stator adopted the structure of the three-phase concentrated windings. The three-phase windings,
A,
B and
C, were arranged in a space at an interval of 120°, and the number of turns was
Ns. The rotor and stator structure of the bearingless motor, as shown in
Figure 4a, is used to analyze the rotation motion of the motor, expand the spatial structure of the motor along the circumferential direction and establish the coordinate system, as shown in
Figure 4b. The direct axis and the quadrature axis are synchronous, dynamic coordinate systems fixed on the motor rotor, which are used to describe the rotating motion of the motor. The magnetic field direction of the motor is the
d axis and the axis perpendicular to the magnetic field direction is the
q axis. Therefore, the spatial phase difference between the
d axis and the
q axis is
π/2. Permanent magnets with different magnetic pole directions are pasted on the lower surface of the rotor of the bearingless motor, and the number
p of rotor pole pairs was four. The three-phase coils of the stator were fed with a sinusoidal current to drive the motor. The magnetic flux distribution
ψ formed by the magnetic field interaction of the permanent magnet and the motor energizing coil is shown by the black dotted arrow in
Figure 4b. In order to study the rotor dynamic characteristics of blood pump, the magnetic field of the bearingless motor was theoretically analyzed using the inductance matrix method. The air gap magnetic energy was calculated by the magnetic flux distribution and inductance matrix, and the torque and axial electromagnetic force were obtained by deriving the rotation direction
θ and
z-axis direction. In order to calculate the inductance and deduce the electromagnetic force and torque model through the effective permeability, the magnetic field generated by the permanent magnet on the rotor was equivalent to a direct current excitation winding with the turns
Nf and current
if.
The electromagnetic force and the torque of the bearingless motor were theoretically modeled, and the assumptions were as follows:
- (a)
The spatial distribution of the magnetomotive force was sinusoidal.
- (b)
The stator surface was smooth and the slot harmonics could be ignored.
- (c)
Magnetic saturation was ignored.
- (d)
The permeability of the stator core was infinite, and the permeability of the permanent magnet and air was considered equal to that of vacuum.
- (e)
The iron loss was ignored.
For the non-salient rotor, the self-inductance of the stator coil was independent of the rotor position angle
θ. Therefore, the self-inductance of the stator coil was constant, and the
d axis inductance and
q axis inductance were equal. The self-inductances of the three-phase stator coils are given by Equation (1):
Therefore, the self-inductance of the stator was composed of
Laa0 and
Lal, where
Laa0 is the self-inductance component generated in the coil by the fundamental flux in the air gap; and
Lal is the additional component generated for the leakage flux of armature winding, which is also called the leakage inductance
Lal. The self-inductance of a winding is generally a function of the air gap
g [
20]. Therefore, the self-inductance of stator winding can be approximately expressed by Equation (2) with air gap
g:
where
L′aa0 is the stator magnetizing inductance multiplied by the air gap length, and the unit is Hm.
L′aa0 and
Lal are constants, which are obtained by the motor parameter identification and the phase inductance measurement. The mutual-inductance between the two windings is the product of the self-inductance and
cosα, whose phase angle difference is the electrical angle
α. Since the angle between the different phases of armature winding was 120°, the mutual inductances between the different phases of armature winding are equal, which is expressed by Equation (3):
Similar to the self-inductance of stator, the self-inductance
Lf of the rotor excitation winding is expressed in Equation (4):
where
Lff0 is the inductance generated by the fundamental component of magnetic flux in the air gap. The leakage inductance
Lfl of the rotor excitation winding is used to characterize the influence of permanent magnet leakage flux. The mutual inductance between the stator and rotor is related to the rotor position and changes periodically with
θ, so the mutual inductance between the excitation winding
f and phase winding
a is given by Eqaution (5):
where
L′m is the constant part of the mutual inductance of the motor stator and rotor. It was observed that the relationship between the self-inductance and the mutual inductance of the stator and rotor that the relationship
was satisfied. When the windings of the motor stator were three-phase symmetrical,
ia +
ib +
ic = 0, and thus
ib +
ic = −
ia. On substituting the self-inductance of the above phase
A and the mutual inductance of the excitation winding
f, the phases
B and
C into the flux linkage expression were obtained as shown in Equation (6):
If the synchronous inductance
Ls = 3/2
Laa0 + Lal is defined, the corresponding flux linkage equation of phase
A is:
Therefore, the phase voltage of the motor is composed of the voltage on the phase resistance and the back electromotive force, and its expression is:
The flux linkage was obtained by the product of the inductance matrix of the motor and the current, as shown in Equation (9):
where
L is the inductance matrix of the motor; the conversion matrix between the 3-phase current and
d axis and
q axis currents is given by Equation (10):
Then, the motor flux linkage equation was converted into:
where,
λm = Lmif is the flux linkage generated by the magnetic field of the rotor permanent magnet in the stator winding. The magnetic energy
W stored in the coil was calculated by Equation (12):
If the motor rotor deviated slightly from the balance position
g0 in the axial direction, the magnetic energy changed accordingly. Therefore, the axial magnetic force and torque were calculated by the partial derivative of the magnetic energy in the axial displacement
g and rotation direction
θ [
21]. The axial electromagnetic force between the motor stator and rotor is:
Additionally, the torque
Ts was expressed as:
2.3. Axial Force Modeling of the Permanent Magnet Bearing
The magnetic field of permanent magnetic bearing was solved by the static magnetic field. The permanent magnetic force was obtained by integrating the Maxwell stress tensor on the surface of a permanent magnetic bearing [
22]. The calculation formula used was:
where
n is the outward normal vector of the outer surface of the permanent magnet bearing and
T is the Maxwell stress tensor on the outer surface of the permanent magnet bearing. The axial component
fz of the permanent magnetic force
fp is the axial force of the permanent magnetic bearing on the rotor, and the direction is vertical upward along the
z axial direction. Its derivative relative to the axial displacement is the magnetic stiffness [
23], which was defined as:
The three-dimensional model of the permanent magnet bearing assembly was established, and the schematic diagram is shown in
Figure 5. The upper concentric magnetic rings were arranged on the blood pump casing, which were static magnetic rings. The lower concentric magnetic rings were arranged on the upper-end face of the impeller and were movable magnetic rings.
The material of the permanent magnetic bearing was the permanent, magnetic, rare earth material NdFeB35, which adopted an axial magnetization. The specific parameters of the permanent magnetic bearing are shown in
Table 2.
The magnetic field of the permanent magnet bearing was solved by the static field module in the COMSOL software, and the magnetic force was obtained by integrating the Maxwell stress surface. When simulating the static magnetic field of the permanent magnet bearing, the axial air gap scanning range was set to 0.1~3.1 mm. In the permanent magnet bearing, the movable magnetic ring was offset under the action of the external load and interference.
Figure 6a shows the magnetic density and magnetic induction line distribution when the air gap
g1 was 0.4 mm and 1.72 mm (Δ
z = 1.32 mm). The air gap was inversely proportional to the magnetic density and the strength of the magnetic field. The variation law of the axial restoring force
Fm1 is shown in
Figure 6b. There was a strong nonlinearity between the axial restoring force
Fm1 and the air gap
g1. The smaller the air gap between the movable magnetic ring and the static magnetic ring, the greater the axial magnetic force. The equilibrium position of the rotor was
g1 = 1.6 mm and the movable range was 1.5–1.7 mm. Within this range, the axial restoring force
Fm1 and axial offset were approximately linear. Therefore, the axial force was linearized near the equilibrium position to obtain the linearization model of the permanent magnet bearing.
The axial force at the different axial displacements was fitted, and the mathematical model between the axial force
Fm1 and the axial displacement
g1 was established, as shown in Equation (17). The root mean square error
RMSE of the fitting results was 0.042 and the correlation coefficient
R was 0.999:
2.4. Model Linearization of Magnetic Suspension System
Without considering the effects of fluid force, gravity and external disturbance, the total axial force on the rotor is given by Equation (18):
When the motor rotor is in a balanced and stable position, the air gap between the stator and the rotor surface is
g0. Assuming that the rotor has an axial offset
z and the rotor deviates from the blood pump inlet along the axial direction, the air gap becomes
g0 z. The axial force and electromagnetic torque applied to the rotor are linearized at the equilibrium point, and Equations (13) and (14) are expanded by the Lagrange linearization method, then the respective linearization models of axial force and torque are:
The air gap
g0 of the motor was 1.8 mm. Due to structural constraints, the maximum displacement on one side of the rotor was 0.1 mm and the axial displacement
z was relatively small during the actual axial position control of the rotor. Therefore,
z/
g0 was approximately 0. In addition,
kt = 3
pL′mif/2
g0 instead of the constant term in Equation (20), Equation (20) is simplified as:
According to Equation (19), the axial force
Fz was coupled by the direct axis current
id and the quadrature axis current
iq. If the coupling term is ignored, the axial force is only related to
id and can be controlled linearly. According to Equation (21), the electromagnetic torque is only related to the quadrature axis current
iq. After linearization, the axial force shown in Equation (19) can be expressed as in the form of
z and
id as:
where,
is the force gain;
is the axial displacement stiffness; other terms of the model can be regarded as disturbance terms, that is,
f . When the impeller was passively suspended in the equilibrium position,
f is 0 N. Without considering gravity and fluid force, according to Newton’s second law of motion, the axial motion of rotor can be expressed as:
According to Equation (23), it can be seen that the rotor of the bearingless motor exhibited a negative stiffness in the axial direction, which is unstable for the rotor suspension state and seriously affects the stable operation of the blood pump.
2.6. Verification of the Magnetic Suspension System Theoretical Model
The transient electromagnetic field of the bearingless motor was solved using the ANSYS Electronics software. The rotor backplate was made of titanium alloy. The permanent magnet material was NdFeB35, the residual magnetic flux density Br was set to 1.3 T, and the different magnetization directions were set according to the magnetic pole orientation. The material of the stator winding was copper and the number of turns was 115. The material of the iron core and yoke was silicon steel. The number of turns and sinusoidal current excitation were set for the coil. The AC current excitation amplitude im was 0.4 A and the spatial phase difference was 120 degrees. The electromagnetic field solution time was set to 10 ms, and the simulation solution time step was the time corresponding to each 1° rotation of the motor.
The axial force received by the rotor of the bearingless motor during rotation is shown in
Figure 8a. The axial force obtained by the FEM simulation fluctuated all the time during the rotation process, and the amplitude was essentially within 0.5 N. The axial force was negative, indicating that the motor rotor was always subject to the suction of the stator. The average axial force was calculated by FEM to be −5.2 N. Assuming that the rotor had no axial offset, the direct axis current
id was 0 A and the quadrature axis current
iq was 0.4 A; the axial force calculated by the theoretical model was −4.63 N. The two results were close, and the error of the theoretical model was 11%. The variation process of the electromagnetic torque during the motor rotation is shown in
Figure 8b. There is an obvious pulsation phenomenon of the rotor torque during rotation. The potential causes of torque ripple were the interactions of the stator and rotor magnetic field, current commutation, the cogging effect, and the armature reaction. The average value of the electromagnetic torque in the FEM simulation results was 6.6 N·mm. According to the theoretical Equation (21), the theoretical torque at 3000 rpm was 6.5 N·mm, and the error was 1.5%. The two values were very close, which verified the accuracy of the motor torque theoretical model.
When the motor rotor was offset under the action of external disturbance, the changes in the permanent magnetic force on the rotor are shown in
Figure 9a. Because the motor structure and the permanent magnet bearings were spatially symmetrical, the permanent magnetic force in the
x and
y direction was always 0 N, when the impeller was offset, and the axial force
Fz changed with the offset. When the impeller was offset to the inlet of the blood pump, the axial force on the impeller was upward, that is, the axial force direction was consistent with the offset direction, and its amplitude increased with an increase in the offset. Therefore, the impeller was unstable in the axial direction. The resultant axial force at the balance point was close to 0 N.
Figure 9b shows the theoretical value and the FEM simulation results of the axial force of the rotor within the axial movable range. It can be seen from
Figure 9b that the theoretical value of axial force was very close to the FEM simulation results, which verifies the accuracy of the theoretical model. Thus, the stiffness calculated by the theoretical model can be used as a reference index for the structural design of the hydrodynamic bearing. The position where the resultant force was 0 N was located below the theoretical equilibrium point. Furthermore, the linear relationship between the axial force and the axial displacement at the equilibrium point also verified the rationality of the linearization of the axial force theoretical model. According to the definition of stiffness, the axial stiffness was negative, which indicated that the impeller was unstable under the action of external disturbance.
2.7. Hydrodynamic Force and Stiffness Model of Hydrodynamic Bearing
The hydrodynamic bearing was located on the inner cavity surface of the upper and lower casings of the blood pump. The fluid force generated by the dynamic pressure effect acted on the rotor to realize the axial stable suspension. Due to the secondary flow field at the upper and lower clearance of the blood pump, the blood flow field was in a spiral form. Although the herringbone groove hydrodynamic bearing had a better stability, the groove type was inconsistent with the blood flow form. When the blood flowed through the hydrodynamic bearing, flow separation and turbulence occurred, which could easily damage the blood. Therefore, in order to obtain a better blood compatibility, the hydrodynamic bearing was designed in the form of a spiral groove. The specific structure is shown in
Figure 10a. The groove type of the hydrodynamic bearing is spiral. According to the rotation direction of the impeller and the flow direction of blood, the hydrodynamic bearing on the lower pump casing was designed as the “Pump-in” type, and the hydrodynamic bearing on the upper casing was designed as the “Pump-out” type. The structural parameters of the two hydrodynamic bearings were consistent. The inlet of the hydrodynamic bearing was located outside the spiral groove with radius
rs. When the impeller rotated in the direction shown in
Figure 10a, the high-pressure fluid emitted by the impeller flowed into the hydrodynamic bearing from the inlet of the spiral groove. As the groove width of the spiral groove decreased in the flow direction, the wedge-shaped spiral groove formed a dynamic pressure effect, resulting in the bearing capacity to realize the suspension of the impeller.
Figure 10b shows the cross-sectional view of the hydrodynamic bearing. The vertical distance between the lower surface of the impeller and the pump casing is the liquid film thickness
h1 of the hydrodynamic bearing, that is, the suspension clearance of the impeller. The groove depth of the spiral groove of the hydrodynamic bearing is
h2, and the sum of
h1 and
h2 is the distance
h between the bottom surface of the impeller and the groove depth. The liquid film thickness
h1 did not only affect the suspension performance of the hydrodynamic bearing, but also affected the hemolysis performance. Therefore, the liquid film thickness was designed after comprehensively considering the performance of the hemolysis and suspension. Here, this factor is not considered temporarily in the structural design of the hydrodynamic bearing, and only the influence of the groove depth
h2 on the suspension performance was analyzed. In addition to the groove depth, the groove number
N, the groove width ratio
δ = Bg/
B and spiral angle
β were the key parameters affecting the hydrodynamic bearing performance.
The spiral groove of the hydrodynamic bearing was in the form of a logarithmic spiral, and its structural shape was closely related to the spiral angle
β, which is expressed as:
where,
rg is the base circle radius of spiral line,
θ is the spiral line control angle, and
β is the spiral angle. If the pressure change in the clearance direction was ignored, the liquid film pressure distribution of the hydrodynamic bearing could be solved based on the Reynolds equation in lubrication theory. The steady-state Reynolds equation is shown in Equation (25):
where,
ρ is the fluid density,
h1 is the clearance,
h1 = hc + hb,
hc = −uc·nref,
hb = hb1 + ub·nref − ub·∇
thb1,
uc is the impeller displacement,
nref is the reference normal direction,
hb1 is the initial clearance,
vave is the average speed, and the calculation formula is:
where
μ is blood viscosity;
vc,t and
vb,t are the velocity vectors of the impeller surface and hydrodynamic bearing wall at the clearance, respectively. The specific calculation method is shown in Equations (27) and (28):
The absolute pressure
pA at any point on the bearing surface was calculated from the sum of the reference pressure
pref and the relative pressure
pf, which is
pA = pref + pf. The pressure distribution on the whole hydrodynamic bearing surface was obtained through an iterative calculation. The pressure of each node on the hydrodynamic bearing surface was double integrated along the radial direction and circumferential direction to obtain the bearing capacity
Ff of the liquid film:
where,
ri and
ro are the inner radius and outer radius of the hydrodynamic bearing, respectively. When the axial displacement of the impeller occurred, the clearance between the impeller and the hydrodynamic bearing changed correspondingly, as did the bearing capacity imposed by the hydrodynamic bearing on the impeller. The change rate is the stiffness of the hydrodynamic bearing, which represents the repulsion of the hydrodynamic suspension bearing to the disturbance, which is defined as: