#### Abstract

Current researches show that the constant speed mode adopted by the existing commercial blood pump may cause damage to the body. The way to solve this problem is to produce pulsating flow by changing the speed of the blood pump’s impeller. But at present, the flow field of the blood pump is not clear, when it changes speed, and the coupling between blood pump and body has not been considered in the simulation of the flow field. A multiscale coupling model combining hemodynamics (0D) and Computational Fluid Dynamics (3D) was established in this paper to solve the problem, and a speed change curve consistent with the ventricular motion was selected. The hemodynamics, shear stress, and hemolysis changes of 6000 rpm at different amplitude (2000, 3000, and 4000 rpm) were simulated, analyzed, and compared with the constant speed (7000 rpm). The results show that the pressure difference obtained by simulation is consistent with the experimental results, and the flow generated by the natural heart still flows through the blood pump, thus changing the working point of the blood pump. When the blood pump works at the changing speed, it could produce more pulsation, and the shear stress and hemolysis in the blood pump increase with the rising of speed and flow. But according to the hemolysis score of a single cardiac cycle, the hemolysis value of the changing speed model at an amplitude of 4000 rpm is only 11.71% higher than that of constant speed at 7000 rpm.

#### 1. Introduction

At present, there are more than 20 million heart failure patients in the world [1]. In the case of poor drug treatment and lack of donors, implanting a blood pump to assist the heart has become an effective treatment method. In 2013, it was listed as the treatment of heart failure (type: IIA) by the ACCF/AHA Guideline for the Management of Heart Failure [2]. In the past ten years, nearly 20000 patients have been implanted the blood pump and achieved a certain survival rate [3]. However, the current commercial blood pump generally adopts the constant speed working mode. The existing researches believe that this mode will reduce the pulsatility of blood flow and then cause damage to the body, such as thrombosis or aortic insufficiency [4], renal failure [5], and vascular dysfunction [6]. It has been proposed to generate pulsating blood flow by adjusting the speed of the blood pump’s impeller to solve this problem [7, 8].

After the blood pump is implanted in the body, it is coupled with the cardiovascular system and then transports the blood by the rotation of the impeller. During the process, there will be a shear flow field inside the blood pump, and the red blood cells may rupture under the effect of strong shear flow, thus inducing hemolysis [9, 10]. To verify the rationality of blood pump structure and speed design, Computational Fluid Dynamics (CFD) is widely used in the flow field simulation of blood pump [11–13]. At present, there are many kinds of research focused on the flow field and hemolysis analysis of the blood pump at the constant rotation speed. However, the study of the flow field under the variable speed of the blood pump is less, and the characteristics of the flow field under the coupling of the blood pump and cardiovascular system are not clear. Drešar et al. [14] conducted a turbulence model simulation of the HeartAssist5 blood pump based on CFD and analyzed the changes in the flow field and hydraulic characteristics at the constant speed. Lopes and Guilherme et al. [15] took an axial-flow blood pump as the research object and analyzed the hemolysis and distribution of the shear stress, while the researches of Drešar and Guilherme did not consider the condition of variable speed. Nammakie et al. [16] analyzed the flow and pressure changes of the blood pump under variable speed based on CFD and proved that the required pulsation could be generated by controlling the rotating speed. Chen et al. [11] analyzed the transient flow field of the axial-flow and centrifugal pump under the pulsating flow condition based on CFD and pointed out the high-risk areas of hemolysis inside. Sonntag et al. [17] carried out CFD simulation on a pulsatile total artificial heart, and the results show that the left ventricular assist pump has good pulsation and washout performance; Yang et al. [18] used a high-speed camera to photograph the flow line of pulsating blood pump. The results show that pulsation will cause periodic confusion in the flow field. Wang et al. [19] analyzed the hemolysis and hemodynamics of the blood pump in the pulsation process by CFD and hemodynamic models, respectively. The results showed that the speed modulation of the impeller and the optimization of the blood pump structure are important. Although the pulsing was analyzed in the above-mentioned studies, the coupling of the blood pump and cardiovascular system was not considered. However, when the blood pump works in the human body, it is dynamically coupled with the cardiovascular system. The blood flow generated by the natural heart will flow through the blood pump, which will cause the working point of the blood pump to deviate from the design conditions. Therefore, the coupling effect of the blood pump and cardiovascular system should be considered to obtain more realistic flow field information.

In this paper, a multiscale coupling simulation model of hemodynamics (0D model) and CFD (3D model) was established in this study to solve the above-mentioned problems. Different from previous studies, this paper analyzed the flow field with the coupling of the pump and cardiovascular system. A speed regulating curve consistent with the ventricular action was selected. The hemodynamic, shear stress, and hemolysis changes were analyzed through the coupling model, and the results were compared with the constant speed. The method of this study can provide theoretical support for obtaining the flow field under the coupling of the blood pump and body. Also, it has reference significance for the formulation of the blood pump’s variable speed strategy.

#### 2. Materials and Methods

##### 2.1. Blood Pump Model

The blood pump used in this paper is developed by the research group of the author, and the pump is an axial-flow blood pump. The structure is shown in Figure 1, which is composed of the front vane, impeller, and back vane. When a specific alternating magnetic field is outside the impeller, the permanent magnet will drive the impeller to rotate and then transport the blood. The front and back vanes are used to correct the flow direction. The diameter of the blood pump is 16.5 mm, and the length is 73 mm.

##### 2.2. Establishment of Coupling Model of Hemodynamics and CFD

The coupling model established in this paper is shown in Figure 2. The hemodynamic model was based on the research of Korakianitis and Shi et al. [20]. Compared with some simplified models [21, 22], the model shows more details for the simulation of the human blood circulation system, which includes systemic circulation, pulmonary circulation, the left and right ventricles, and atria. The systemic circulation and pulmonary circulation circuits include arteries, veins, arterioles, capillaries, and aortic sinus and pulmonary sinus, respectively. The following equivalences and assumptions were made in the model: the blood flow resistance is equal to resistance, blood vessel volume is equivalent to capacitance, blood vessels is equal to compliance inductance, pressure is equal to voltage, and flow rate is equivalent to current. The model simplifies the complex human vascular system. Thus, the model may not be able to simulate the specific lesions produced by some tissues. However, the purpose of introducing this model is to obtain the coupling effect between the blood pump and the cardiovascular system. From the previous studies on the influence of the blood pump on hemodynamics [23, 24], this simplified model can meet the needs, so this simplification can be considered reasonable. The model of the systemic circulation is shown in (1), where elv(*t*) and ela(*t*) represent ventricular and atrial activation functions. Models of pulmonary circulation are similar to the systemic circulation, and all details can be found in Shi’s paper [20]. One has

In the hemodynamic model, the blood pump was equivalent to the current source. We improved the blood pump model by introducing the derivative terms of pressure and speed so that the blood pump model can reflect the characteristics of variable speed more accurately. The original model [25] and the improved model are shown in (2) and (3), where is the flow through the blood pump, *K*_{0}∼*K*_{8} is the coefficient of the equation, is the pressure difference between the inlet and outlet of the blood pump, and is the speed of the blood pump. Besides, the hemodynamic model in this paper is established in MATLAB 2016b, CFD model is built in FLUENT 17.0, and the data transmission within the software is realized by m-function and UDF. The realization process of the coupling model is shown in Figure 3. First, the hemodynamics and blood pump’s flow can be obtained in the hemodynamic model. After that, MATLAB transmits the blood pump’s flow and speed information to the CFD model. Finally, this information is used as boundary conditions for CFD simulation, and the flow field information of the blood pump can be obtained.

To verify the correlation between the input and output of the model, we used the SIMLAB software to analyze the global sensitivity of each input of the equation. The Spearman coefficient (SPEA) was used to analyze the sensitivity of each input [26]. If the absolute value of the SPEA coefficient is larger, the relationship between input and output is stronger:

##### 2.3. Construction of Test Bench

As shown in Figure 4, an experimental test system was built. There are two primary purposes to build the experimental system: one is to obtain the hydraulic characteristics of the blood pump under different conditions and then to fit the parameters of the blood pump model; the other is to add additional flow to adjust the operating point of the blood pump, so as to verify the results of the blood pump pressure difference obtained by the coupling model. The auxiliary pump was used to provide extra flow to regulate the operating point of the blood pump, and the gear flowmeter was installed to measure its flow. The damping valves were used to adjust the resistance of the pipeline. The pressure sensors were installed at the inlet and outlet of the blood pump. The ultrasonic flowmeter was installed in the circuit to measure the total flow. The thermostat water bath cauldron was used to ensure that the temperature of the fluid is stable at a specific value. The characteristics of the equipment can be found in Table 1. In the experiment, a 33.3% aqueous glycerol solution was used, and its density and viscosity are 1.08 kg/m^{3} and 3.5 cps at 25 degrees, which was close to the whole blood [27].

In the experiment of fitting the parameters of the blood pump model, the auxiliary pump does not work, and the damping valve ten is completely closed. Besides, to make the sample more comprehensive, the sine wave (see (4)), square wave (see (5)), triangle wave (see (6)), and ventricular activation function (see (7)) were chosen as the speed regulating waveform of the blood pump, where *T* is the cycle, *N*_{b} is the basic speed (6000/7000/8000 rpm), *N*_{a} is the speed amplitude (500–4000 rpm), and the speed range is 4000–12000 rpm through combination. After obtaining the sample, the least square method (see (8) and (9)) was used to fit the model parameters, where *x*_{i} represents the *i*th term of the equation, *m* represents the sample size, and *n* represents the number of the ith term of the equation. As shown in (10) and (11), root mean square error (RMSE) and the mean absolute percentage error (MAPE) are used to evaluate the fitting effect, where *Q*_{TEST, m} represents the predict flow rate, *Q*_{REAL, m} represents the actual flow rate, and *m* represents the number of dates. One has

##### 2.4. CFD Simulation Settings

The CFD simulation was based on FLUENT 17.0. The SST k- model was selected as the turbulence model, for the SST k- model contains the modified turbulent viscosity formulas and considers the effect of turbulent shear stress; it can simulate the near-wall flow field more accurately and proved to be suitable for blood pump simulation [15].

The expression of the *k*- SST model [28] is shown in (12) and (13), where *k* is the turbulent kinetic energy, is the fluid density, is the dynamic viscosity, is the specific dissipation rate, is turbulent eddy viscosity, and *F*_{1} is the blending function, which is used to blend the k- model (*F*_{1} = 0) and k- model (*F*_{1} = 1), and the model coefficients are as follows: = 0.075, = 0.09, *σ*_{k} = 2, = 2, and = 1/0.856. One has

The boundary condition of the inlet was set as flow, the outlet boundary condition was set as pressure, and the interface method between impeller and guide vanes was set as the frozen rotor. A second-order upwind advection scheme was set in all cases, with double precision. For time discretization, a second-order backward Euler transient scheme was used. The residual target was set to 10^{−5}.

##### 2.5. Mesh Generation

In this study, ICEM CFD was used to generate the mesh of the blood pump. The unstructured tetrahedral mesh was used for the grid generation, which has good adaptability to the complex shape. The requirement of Y + for the SST k- turbulence model is less than or equal to 1. Then the boundary layer grid was divided according to the requirement, the growth factor was 1.05, and the total number of layers was 20. After verifying the independence of the number of grids, it was determined that the total number of mesh was 4.5 million. The mesh is shown in Figure 5.

##### 2.6. Selection of Speed Curve of the Blood Pump

In the study, the speed of the blood pump was divided into two cases: constant speed and variable speed. 7000 rpm was chosen as the constant speed. According to Pirbodaghi et al. [29], when the amplitude, phase, and frequency are fixed, different waveforms have little effect on the performance of the blood pump. As the ventricular activation’s function can reflect the change of ventricular action, so it was selected to regulate the speed. The speed waveform is shown in Figure 6. The basic speed was set to 6000 rpm to prevent the blood from flowing back. The amplitude of speed change was set as 2000 rpm, 3000 rpm, and 4000 rpm. The cycle *T* of speed change was equal to the cardiac cycle. That is, the change of blood pump speed is consistent with the ventricle action: the speed increases when the ventricle contracts and decreases when the ventricle dilates.

##### 2.7. Calculation of the Scalar Shear Stress and Hemolysis (*HI*)

As shown in (14), Giersiepen et al. [30] established a hemolytic model through a large number of experimental studies. The formula gives the functional relationship between hemolysis, shear stress scalar, and exposure time, where represents shear stress scalar and *t* is the exposure time under the shear stress. The scalar shear stress can be calculated by (15) and (16) [31], where is the fluid viscosity; *u*, , and are the velocity vector in *x*, *y*, and *z* directions; is the shear stress; and *s*_{ij} is the Reynolds shear stress. The above-mentioned values can be extracted from the CFD simulation results. In the hemolysis calculation, 500 streamlines were selected at each time point, and the hemolysis of each streamline was calculated in MATLAB. The values of these 500 streamlines were averaged to obtain the hemolysis value at the time point:

#### 3. Results

##### 3.1. Fitting Results of Blood Pump Model

In the parameter fitting experiment, 119000 groups of data were collected, 100000 groups of data were selected as training samples, and the remaining data samples were used to test the accuracy of the model. The model parameters obtained are shown in Table 2.

The sensitivity coefficient of each input of the model is shown in Figure 7. In the figure, IN_{i} (*i* = 1–8) corresponds to the second to ninth input of (3), respectively. It can be seen from the results that IN_{2} has the most significant correlation with the output (correlation coefficient is 0.6563). The SPEA of the remaining inputs distributes between 13% and 20%, while the sensitivity coefficient of IN_{6} is the smallest, which is 13.90%. The RMS error and the mean absolute percentage error of the improved model are 6.3148 and 0.075, respectively, which are better than those of Shi’s model (9.679 and 0.128).

##### 3.2. Changes in Blood Pump Flow and Physiological Pressure under Coupling Condition

The left ventricular pressure (P_{lv}) and aortic pressure () obtained by the model are shown in Figure 8(a). During the simulation, the natural heart was operated in the state of heart failure by setting the model parameter. Before the assistance of the blood pump, is 49.83 mmHg-70.51 mmHg, and Plv is 19.56 mmHg–71.13 mmHg. After the support of the blood pump, recovered to 102.50 mmHg–105.09 mmHg at the constant speed of 7000 rpm and to 89.72 mmHg–96.46 mmHg, 91.67 mmHg–100.04 mmHg, and 93.59 mmHg–103.48 mmHg at the variable speed (corresponding to the amplitude of 2000 rpm, 3000 rpm, and 4000 rpm, resp.). Besides, Plv decreased at all rotational speeds, indicating that the ventricle was unloaded, which is helpful for the recovery of ventricular function.

**(a)**

**(b)**

The change of blood pump flow is shown in Figure 8(b). At the constant speed of 7000 rpm, the changing trend of blood pump flow is consistent with the change of ventricular action. At 0–0.3 T, the ventricle is in systole; the flow increases with time and reaches the maximum value of 6.42 L/min at 0.3 T; from 0.3 T to 0.45 T, the ventricle is in the diastolic phase, and the flow rate decreases with time; during the 0.45T-T, the flow of the blood pump changes little because the ventricle enters isovolumic diastolic phase. In variable speed, the changing trend of blood pump flow is the same as that of constant speed. The peak flow corresponding to the amplitude of 2000 rpm, 3000 rpm, and 4000 rpm is 8.71 L/min, 9.87 L/min, and 11.04 L/min, respectively. Also, the flow rate during the constant volume contraction period is lower than that at the constant speed (7000 rpm), which is due to the lower speed (6000 rpm). In the above-mentioned results, the flow of blood pump showed the same change trend with the ventricular action; the reason for this phenomenon is that the flow generated by the natural heart flows through the blood pump.

##### 3.3. Verification of the Blood Pump’s Pressure Difference

The blood pump’s hydraulic characteristics in the final CFD simulation were selected to verify the model. In the verification process, several operating points were chosen in each speed range, and the cardiac flow at this point was extracted through the hemodynamic model. In the experiment, the auxiliary pump was adjusted to make the output flow equal to the cardiac flow obtained from the hemodynamic model. Then the pressure difference of the blood pump was recorded. The simulation and experimental results are shown in Figure 9, the maximum error is about 15%, and most of the errors are within 10%. It can be seen from Figure 9(d) that when the pump worked at the constant speed, the pressure difference of blood pump decreases with the increase of blood pump flow, which is related to the working characteristics of the blood pump. The pressure difference generated by the pump will decrease with the rise of flow at a certain speed. Although the power of the pump will increase when the speed rises, the pressure difference in Figure 9(a)-9(c) still represents the same change trend with the constant speed, which shows that the power increased is limited compared with the natural heart flow through the blood pump.

**(a)**

**(b)**

**(c)**

**(d)**

##### 3.4. Shear Stress Distribution on the Impeller

The distribution of the shear stress on the impeller obtained by the model is shown in Figure 10. For comparison, four points of 0.15 T, 0.3 T, 0.45 T, and 0.8 T were selected, which correspond to the four stages of flow change (i.e., rise, peak, fall, and basically unchanged). In general, the maximum shear stress is concentrated at the beginning of the spiral blade on the impeller, which is due to the thin thickness here (about 0.1 mm). When the impeller is rotating, there will be a large velocity gradient and shear force in this area. From the changing trend of shear stress, the shear stress at all speeds reaches the maximum at 0.3 T, decreases at 0.15 T and 0.45 T, and reaches the minimum at 0.8 T. Also, under the condition of variable speed, it can be observed that the shear stress and rotation speed show a positive correlation.

**(a)**

**(b)**

**(c)**

**(d)**

##### 3.5. Shear Stress Distribution inside the Blood Pump

To further analyze the shear stress distribution in the blood pump, the proportions of fluid in different shear stress range were extracted, and the selection of time points was the same as that in Sections 3 and 4. According to Chen et al. [32], the scalar shear stress was normal in the range of 0 to 25 Pa, slightly nonphysiological in the range of more than 25 Pa and less than 125 Pa, and very nonphysiological in the range of more than 125 Pa. The simulation results are shown in Table 3. At 0.15 T, the shear stress at all speeds is basically distributed between 0 and 25 Pa, only about 3%–5% of the stress is in the range of 25Pa-50 Pa, and few are in the field of 50Pa–100 Pa. At 0.3 T, the proportion in 0–25 Pa decreases, and the shear stress at constant speed is smaller than that at variable speed; under the condition of variable speed, the percentage of this range (0–25 Pa) will decrease with the increase of amplitude. Correspondingly, the shear stress level will increase with the ascend of amplitude in each range higher than 25 Pa. At 0.45 T, the shear stress at all rotational speeds falls back to the range of 0–25 Pa generally. At 0.8 T, the shear stress of all rotating speeds is still mostly in the range of 0–25 Pa, but in this range, the proportion of variable speed is higher than that of constant speed. The reason is that the working condition of variable speed is 6000 rpm, which is lower than that of constant speed (7000 rpm).

##### 3.6. Hemolysis Change in a Single Cardiac Cycle

The hemolysis index (*HI*) change of each speed mode in a cardiac cycle is shown in Figure 11(a). Generally speaking, the changing trend of the curve is basically consistent with that of rotation speed and flow. For the constant speed condition, the *HI* rises with the increase of flow rate at 0–0.3 T and reaches the maximum of 0.0097 at 0.3 T; within the time of 0.3 T to 0.6 T, the *HI* decreases with the decrease of flow rate; after 0.6 T, the *HI* maintains near 0.0016.

**(a)**

**(b)**

Compared with the constant speed condition, the *HI* of the variable speed condition is lower at the early acceleration stage. With the rising of speed, the *HI* increases, and it begins to be larger than the constant speed at about 0.2 T and reached the maximum value at 0.3 T; the maximum hemolytic values of each amplitude are 0.0119, 0.0145, and 0.017 (corresponding to 2000 rpm, 3000 rpm, and 4000 rpm, resp.). Subsequently, the *HI* decreased with the reduction of rotating speed in the range of 0.3 t-0.6 t. After 0.6 T, the *HI* of each amplitude remains at about 0.0004 due to the stable speed and flow. In a single cardiac cycle, the cumulative hemolysis value (*HI*_{acc}) can be calculated after integrating the *HI* and time. The obtained results are shown in Figure 11(b), the *HI*_{acc} of constant speed is 0.00333, and the value of variable speed is 0.00295, 0.00323, and 0.00372 under the amplitude of 2000 rpm, 3000 rpm, and 4000 rpm, respectively.

#### 4. Discussion

The selection of blood pump speed is based on the return of aortic pressure to a normal level without any backflow. It can be seen from Figure 8 that, with the assistance of the blood pump, the aortic pressure in the state of heart failure can be restored to normal. It shows that the selection of rotation speed is within a reasonable range. At constant and variable speed, the flow’s changing trend is consistent with the ventricular’s action. The flow fluctuation of this part comes from the natural heart. For the variable speed condition, the pulsation of blood pump flow and aortic pressure is higher than that of the constant speed condition, and its value will increase with the rising of amplitude. The simulation research of Gao et al. [33] obtained similar conclusions, and the experiment results of Bozkurt [34] also proved that.

It can be seen from Figure 9 that the pressure difference of the blood pump decreases with the increase of flow, which is easy to understand at the constant speed. But we also observed this rule under the condition of variable speed. Different from the constant speed, if the resistance is specific and there is no external flow, the pressure difference and flow will increase with speed rising (determined by the characteristics of the blood pump; Shi et al. simplified the vascular resistance in the model and regarded it as a constant value). Still, the pressure difference decreases, which is obviously due to the flow generated by the natural heart flowing through the blood pump. Therefore, it should be noted that when the blood pump is coupled with the cardiovascular system, the flow generated by the natural heart will flow through the blood pump. This will cause the flow of the blood pump to be higher than the design condition in vitro. When the blood pump in this paper is tested in a separate circulating pipe (similar to Figure 4, but without the auxiliary pump), the flow rate is about 4.9 L/min at the rated speed of 8000 rpm, and the pressure difference is about 13.2 kPa. However, the simulation results show that once the pump is coupled with the cardiovascular system, the flow rate can reach 6.42 L/min at 7000 rpm; this means that the actual working condition deviates from the designed working point. Therefore, the confluence effect of the natural heart should be considered in the design of the blood pump, so as to achieve better flow field and structure optimization.

The distribution of shear stress obtained in Table 3 is consistent with the changing trend of rotation speed and flow rate. With the increase of rotation speed, the proportion of high shear stress will increase, which is consistent with reality, because with the growth of rotation speed, the work of impeller on the fluid will increase, resulting in a higher velocity gradient. Also, the results show that most of the shear stress is distributed in the range of 0–50 Pa under different rotating speeds, and the shear stress in this range is considered to have little damage to cells in physiology.

In a single cardiac cycle, the hemolysis results in Figure 11 show that there is a significant fluctuation in the hemolysis value of variable and constant speed conditions, and the maximum value is about three to five times the minimum value. It can be seen that both the speed of the blood pump and the flow of the natural heart have a significant influence on hemolysis. In terms of cumulative hemolysis in a single cardiac cycle, we can see that the variable speed condition (at the amplitude of 4000 rpm) is only 11.71% (0.00039) higher than the constant speed (7000 rpm). Still, the variable speed condition can achieve better pulsation, so it can be considered that the aortic pressure and flow pulsation can be achieved without a significant increase of hemolysis by setting a reasonable variable speed range.

#### 5. Conclusions

In this paper, a multiscale coupling model combined with hemodynamics and CFD model is established. In the design of the structure and speed curve of the variable speed blood pump, the flow field and hemolysis caused by the coupling between the blood pump and the heart should be considered, and the impact of this coupling makes design more difficult. The model can obtain the flow field in the pump under the coupling of the pump and the cardiovascular system, which is helpful in analyzing the above-mentioned process.

Through model simulation, we can see that the blood flow generated by the natural heart will flow through the blood pump, regardless of whether the pump is working at a constant speed or a variable speed. This will change the actual working point of the blood pump, and this should be considered during the design and optimization of the blood pump. Through the comparison between the speed change curve and the constant speed, the cumulative hemolysis in a single cardiac cycle does not increase significantly. Therefore, by choosing a reasonable speed control mode, the pressure fluctuation can be achieved by changing the speed of the blood pump, and the hemolysis level will not be significantly increased. This study provides a simulation method for the analysis of the flow field, when the blood pump is coupled with the body and provides support for the research of the speed modulation of the blood pump.

The following limitations also exist in this study: Due to the limitations of model development, this study only analyzed the axial blood pump. However, because the working principles of the centrifugal pump and the axial flow pump are similar, the obtained model can still be applied to the centrifugal blood pump. The rotation speed selected in this paper restores the maximum value of aortic pressure to the normal level of about 110 mmHg. If a higher aortic pressure level is required, it can be achieved by increasing the rotation speed. This study mainly discusses the analysis and establishment of coupling models.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This research was supported by the National Nature Foundation of China (Grants nos. 51475477 and 31670999) and the Postgraduate Independent Exploration Project of Central South University (Grant no. 2018zzts144); all support is gratefully acknowledged.