**Numerical Study on the Influence of Step Casing on Cavitating Flows and Instabilities in Inducers with Equal and Varying Pitches**

#### **Lu Yu <sup>1</sup> , Haochen Zhang <sup>1</sup> , Hui Chen <sup>2</sup> , Zhigang Zuo 1,3,\* and Shuhong Liu 1,\***


Received: 17 August 2020; Accepted: 1 September 2020; Published: 4 September 2020

**Abstract:** It is known that cavitating flow characteristics and instabilities in inducers can greatly impact the safety and stability of a liquid rocket. In this paper, step casing optimization design (Model OE and Model AE) was carried out for two three-bladed inducers with an equal (Model O) and a varying pitch (Model A), respectively. The unsteady cavitation flow field and accompanied instabilities were studied via numerical simulations. Reductions of the cavity size and fluctuation were observed in cases with a step casing. A significant difference in cavity structures was seen as well. Referring to the pressure distributions on the blades and details of the flow field, the mechanism of cavitation suppression was revealed. This work provides a feasible and convenient method in engineering practice for optimizing the characteristic of the cavitating flow field and instabilities for the inducer.

**Keywords:** inducer; step casing; varying pitch; cavitating flow and instabilities

### **1. Introduction**

In hydraulic systems of liquid rocket engines, turbopumps are the main hydraulic components that convey fuel and oxidizers. Due to the demand for a maximum power/weight ratio of the main pump, cavitation may develop on the suction sides of the blades at the inlet of the main impeller [1,2], leading to performance degradation (i.e., sharp decrease in head and efficiency) and hydraulic instabilities (i.e., significant pressure fluctuations). Installing an inducer at the upstream of the main pump is a common solution to mitigate the aforementioned effects. Therefore, the inducer often operates under cavitation conditions accompanied by complex cavitation instabilities [3–7]. Among them, rotating cavitation (RC) is widely considered as a major cause for the premature cutoff of engines [8,9].

Great efforts including inducer impeller optimization and casing modifications were made to alleviate the influences of cavitation and its associated instabilities. Considering the manufacturing convenience and the degree of overall structural change, casing modification, especially step casing design, has been widely studied [10–14]. With a local clearance enlargement, it tends to be an effective and realizable way for performance improvement. Kamijo et al. [4] designed five casings with upstream/trailing edge enlargements and proposed a criterion for RC suppression in a LE-7 LOX (Liquid Oxygen) turbopump inducer. Hashimoto et al. [10] experimentally illustrated the influence of step casing with an upstream enlargement. It was observed that the onset and occurrence region of RC were modified effectively. Furthermore, Shimagaki et al. [15] described the mechanism for

RC suppression with step casing (i.e., a decrement of incidence angel by a thickened backflow vortex with a wider tip clearance) by combining particle imaging velocimetry (PIV) and high-speed photography. Moreover, Fujii et al. [16] carried out experiments of 8 step casings to study the effects of geometries (depth, location, etc.). Therefore, it can be concluded that step casing with an upstream enlargement may be a possible method for RC suppression. However, concrete analysis is still needed for specific inducers.

In this paper, step casing optimization design was carried out for two three-bladed inducers with equal and varying pitch, respectively. Emphasis was exerted on the unsteady cavitation flow field and the accompanied instabilities by numerical simulations. The characteristic of the cavitation area on the blades and three-dimensional cavity structures was analyzed. The reductions of the cavity size and RC fluctuation were observed in cases with step casing. Referring to the pressure distributions on the blades and the details of the flow field, the mechanism of cavitation suppression was revealed.

#### **2. Numerical Studies**

#### *2.1. Numerical Methods*

The computational domain is shown in Figure 1, which consisted of an inlet pipe, an annulus inlet casing, an inducer impeller, and an outlet pipe. An annulus casing with deflectors qas placed upstream of the inducer and aimed at forming a uniform inlet flow [17]. The inlet and outlet pipes were extended for a fully developed flow simulation (5 times the diameter of the annulus inlet casing and 7 times the diameter of the inducer blade tip, respectively).

**Figure 1.** Computational domain.

Reynolds averaged Navier–Stokes (RANS) equations were solved in ANSYS CFX 19.2 code. A SST (Shear Stress Transfer) *k* – ω turbulence model [18] and Zwart-Gerber-Belamri (ZGB) cavitation model [19] were utilized for turbulence and cavitation modeling with the boundary conditions of mass flow rate and pressure at the inlet and outlet, respectively.

The entire mesh system was developed in a commercial software package ICEM-CFD (Integrated Computer Engineering and Manufacturing code for Computational Fluid Dynamics). Unstructured hybrid grids were utilized for the annulus inlet casing and inducer. As for the inlet and outlet pipes, hexahedral grids were developed. A 10-layer boundary layer grid was utilized with a height of 0.01 mm at the first layer. More than 20 layers of mesh were set on the tip to accurately capture the clearance flow characteristic.

Unsteady simulations were applied to evaluate the characteristics of cavitation instabilities. The grid and time independence studies were carried out to find the best compromise between

accuracy and efficiency. The mesh system with 2.1 <sup>×</sup> <sup>10</sup><sup>7</sup> elements (Figure 2 shows the mesh of the inducer) and the time step corresponding to 2◦ of the inducer impeller rotations were chosen for the following simulations. The verification of the above numerical setup can be referred to in our previous research [20].

**Figure 2.** Mesh of the inducer.

#### *2.2. The Inducers and Step Casing Design*

Two inducers with an equal (Model O) and varying pitch (Model A) were selected to compare the influence of step casing design (Model O and Model OE, Model A and Model AE). Tables 1 and 2 summarize the geometries of the studied inducers. For the equal pitch inducer, the blade profile was a straight line (shown on the left of Figure 3) leading to an equal inlet and outlet blade angle. As for the varying pitch inducer, the inlet and outlet blade angle was reduced and increased, respectively. Therefore, the blade profile of the varying pitch inducer on the right side of Figure 3 was curved. It should be pointed out that both the inlet and outlet blade angles of the two inducers were different and thus results between the equal and varying pitch were not comparable. The study in this paper mainly focused on the influence of step casing on the equal and varying pitch inducers, respectively, and the cross-comparison was not involved.




**Table 2.** Geometrical parameters of the cases with the varying pitch inducer.

**Figure 3.** Schematic of the blade profile. (Left one for the equal pitch inducer and the right one for the varying pitch inducer).

Figure 4 shows the schematics of the two casings used in this study. Figure 4a terms a step casing, with enlargement near the leading edge of the inducer. While the other (shown in Figure 4b) is a straight one, owning a constant diameter from the inlet to the outlet.

**Figure 4.** Schematics of the studied casings: (**a**) Step casing; (**b**) straight casing.

#### **3. Results**

#### *3.1. Characteristics of Cavity Oscillation on the Blades*

The oscillation of the cavity on the blades was monitored to understand cavitating instabilities. Figure 5 shows the variation of a nondimensional cavity area *Scav '* over time on three blades of Model O. Tthe ratio of the cavity area on the blade to the area of the inlet, =*Scav*/*S<sup>f</sup>* ; from the perspective of cavitation suppression, was expected to be as small as possible. It can be seen that the cavity variation on each blade was similar. They changed in an approximate sinusoidal pattern with a frequency of 152 Hz (~*f* <sup>0</sup>) and a certain phase difference between neighboring blades. It can be inferred that there was rotating cavitation in the studied inducer. Such asymmetric distribution of cavities might influence the operation's safety and stability.

**Figure 5.** The variation of the nondimensional cavity area over time on three blades of Model O.

Taking the result of a certain blade as an example, the impact of the step casing design on the cavity oscillation is discussed below. Figures 6 and 7 show the comparison of cavity oscillation in Model O and OE (cases with the equal pitch inducer), as well as Model A and AE (cases with the varying pitch inducer). The influence of step casing tended to be similar both in equal and varying pitch inducers. Significant reduction of cavities was seen when step casing was applied, showing that the step casing design had a positive effect on cavitation suppression both in equal pitch and varying pitch inducers.

**Figure 6.** The variation of the nondimensional cavity area over time on Blade 1 of Model O and OE (cases with the equal pitch inducer).

**Figure 7.** The variation of the nondimensional cavity area over time on Blade 1 of Model A and AE (cases with the varying pitch inducer).

#### *3.2. Characteristics of Three-Dimensional Cavity Structures*

#### 3.2.1. The Results of Cases with the Equal Pitch Inducer

Three-dimensional cavity structures were analyzed with the help of the isosurface of the vapor volume fraction α*<sup>v</sup>* = 0.3. As shown in Figure 8a,b, the cavity mainly appeared near the leading edge in both Model O and OE. The cavity structures were similar. Cavities in both Model O and OE had two subregions with different characteristics: Type I resembled the attached cavitation region and Type II resembled cloud cavitation. With the step casing, substantial suppression of the cavity in both streamwise and axial directions was figured out. Next, we show details regarding the cavity on three typical sections along the streamwise direction within 1 revolution, thus demonstrating the varying features of cavity structure over time.

**Figure 8.** The 3D cavity structures of (**a**) Model O and(**b**) Model OE.

As shown in Figure 9a, the first section (S1) was in the middle of Type I. θ<sup>1</sup> = θ<sup>0</sup> + <sup>1</sup> 2 ∆θ<sup>I</sup> ,θ<sup>0</sup> indicates the inception position of Type I; ∆θ<sup>I</sup> indicates the region Type I occupies streamwise, indicating the characteristic of Type I cavitation. Figure 9b,c show the variation of cavities within 1 revolution (T) at S<sup>1</sup> for Model O and OE. Here, PS stands for the pressure side, whereas SS denotes the suction side. Furthermore, *r'* indicates the ratio of radial position and tip radius (=*r*/*r*T), while *Z'* indicates the ratio of axial position and the axial length of the blade (=*Z*/*L*B). With a weak radial fluctuation and almost no axial fluctuation, the variation pattern of Type I in each case tends to be similar. We found the cavity area as O ≈ OE. θ1*OE* > θ1*<sup>O</sup>* indicated that the application of step casing may result in downstream cavitation inception.

**Figure 9.** The variation of cavities over time at S1. (**a**) Schematic of S<sup>1</sup> , (**b**) Model O, and (**c**) Model OE.

The section at the end of the sweep was chosen as S2, thus finding that θ2*<sup>O</sup>* = θ2*OE* = 110◦ . For Model O and OE, S<sup>2</sup> located the region of Type II with significant radial and axial variations

(Figure 10a). Obvious morphological differences of cavities were observed. Cavities in Model O (Figure 10b) developed in a much wider range along the radial and axial direction and were limited by the casing. However, in Model OE (Figure 10c), with a local increment of clearance, the cavities were far away from the casing and presented a semi-elliptic shape.

**Figure 10.** The variation of cavities over time at S<sup>2</sup> . (**a**) Schematic of S<sup>2</sup> , (**b**) Model O, and (**c**) Model OE.

The section of average cavitation oscillation along streamwise was taken as S<sup>3</sup> (Figure 11a), which led to the lack of results. In Model O and OE (Figure 11b,c), the radial and axial fluctuations of the cavities were captured, and there was little difference in morphology as the blade was away from the step in Model OE.

**Figure 11.** The variation of cavities over time at S<sup>3</sup> . (**a**) Schematic of S<sup>3</sup> , (**b**) Model O, and (**c**) Model OE.

Table 3 summarizes the detailed characteristics of the cavity of cases with an equal pitch inducer along the streamwise, radial, and axial directions. It can be seen that the cavities in Model O and OE consisted of two regimes—Type I and II—and there were slight differences in the inception location and occupied range along the streamwise direction. However, there was an increase in the maximum range along the radial direction while the maximum range along the axial direction decreased.


**Table 3.** Geometrical parameters of the studied inducers.

3.2.2. The Results of Cases with the Varying Pitch Inducer

Three-dimensional cavity structures for the varying pitch inducer are shown in Figure 12a,b. Unlike the equal pitch inducer, a significant difference in cavity structures were seen in the varying pitch inducer when the step casing was applied. Cavities in Model A owned two subregions. In Model AE, only clearance cavitation was observed and nearly no leakage cavitation occurred in other cases. Thus, there was only the Type I cavity.

**Figure 12.** The 3D cavity structures of (**a**) Model A and (**b**) Model AE.

Cavity details varying pitch cases were analyzed. Figure 13b,c shows the characteristics at S<sup>1</sup> (Figure 13a). For Model A, the midsection of Type I (θ<sup>1</sup> = θ<sup>0</sup> + <sup>1</sup> 2 ∆θ*<sup>I</sup>* ) was taken. As for AE, θ1*AE* = θ0*AE* + <sup>1</sup> 2 ∆θ*IA*. A weak radial fluctuation and almost no axial fluctuation of cavity variation was observed in both cases. In terms of the cavity area, we found that A ≈ AE. θ1*AE* > θ1*<sup>A</sup>* indicates that the application of step casing could result in downstream cavitation inception.

**Figure 13.** The variation of cavities over time at S1. (**a**) Schematic of S<sup>1</sup> , (**b**) Model A, and (**c**) Model AE.

The results at S<sup>2</sup> (Figure 14a, θ2*<sup>A</sup>* = θ2*AE* = 110◦ ) are shown in Figure 14b,c. A striking influence of step casing was observed. In Model A, the cavity variations were similar to Model O and OE, indicating a Type II cavitation. While in Model AE, the fluctuation in the axial direction greatly reduced and tended toward Type I.

**Figure 14.** The variation of cavities over time at S2. (**a**) Schematic of S<sup>2</sup> , (**b**) Model A, and (**c**) Model AE.

At S<sup>3</sup> (Figure 15a), significant differences of morphology and fluctuation characteristics were observed (Figure 15b,c). Apparent radial and axial fluctuations of the cavities were captured in Model A, indicating a Type II cavitation. In Model AE, the cavities tended to be Type I

**Figure 15.** The variation of cavities over time at S<sup>3</sup> . (**a**) Schematic of S<sup>3</sup> , (**b**) Model A, and (**c**) Model AE.

Table 4 summarizes the detailed characteristics of the cavity along the streamwise, radial, and axial directions. For varying inducer cases (Model A and AE), close inception location and total occupied range along the streamwise were observed, whereas in Model AE, the maximum range along the radial and axial direction decreased significantly. As mentioned above, cavities in Model AE mainly manifested as Type. Model A, however, owned Type I and II.


**Table 4.** Geometrical parameters of the studied inducers.

#### *3.3. Characteristics of Blade Loading and Flow Field*

To further analyze the influence of the step casing on the flow field, Figures 16 and 17 show the blade loading Cp (equals to the difference between the static pressure and the inlet pressure at each position divided by the dynamic pressure) along the streamwise direction at a 95% span on certain blades. The distribution on the other two blades was similar. In both the equal and varying pitch inducers, with a local increment in the clearance in the vicinity of the leading edge (LE), the impact of the step casing was shown near the leading edge (0~0.25). The reduction of the maximum pressure and pressure difference between SS and PS was observed. In the subsequent area (>0.25, till the trailing edge (TE)), the step casing design showed little influence on blade loading.

**Figure 16.** Blade loading along the streamwise direction at a 95% span for a certain blade in Models O and OE.

**Figure 17.** Blade loading along the streamwise direction at a 95% span for a certain blade in Models A and AE.

Tables 5 and 6 demonstrate the characteristics of the clearance flow at S2. With the reduction of ∆*p* (pressure difference between SS and PS) and the increment of local tip clearance when step casing was applied, the average velocity at the clearance significantly reduced, while the leakage flow rate increased (145% in Model OE and 56% in Model AE).

**Table 5.** Characteristic parameters of clearance flow in the cases with the equal pitch inducer.


**Table 6.** Characteristic parameters of clearance flow in the cases with the varying pitch inducer.


Furthermore, the distribution of the vapor volume fraction α*<sup>v</sup>* and velocity vectors are shown in Figures 18 and 19 (Figure 18a for Model O, Figure 18b for Model OE. And Figure 19a for Model A, Figure 19b for Model AE). Clearance flow in cases without step casing (Model O and Model A) tended to penetrate upstream and interact with the main flow, leading to a more apparent axial extension of cavities. As for Model OE and AE, with lower velocity and higher leakage flow rates, the clearance flow tended to interact with the main flow in a closer region.

**Figure 19.** Flow fields at S2. (**a**) Model A and (**b**) Model AE.

#### **4. Conclusions**

To evaluate the influences of step casing design on cavitating flows and instabilities in inducers with equal and varying pitches, a numerical simulation was employed in four cases. We analyzed cavitation area characteristics on blades and three-dimensional cavity structures. The step casing design showed a positive effect on cavitation suppression in both equal and varying inducers with a significant reduction of cavity size and fluctuation. Its impact on cavity structures, however, was quite different. In Models O and OE (cases with the equal pitch inducer), there were two regimes: Type I (resembled attached cavitation) and Type II (resembled cloud cavitation). While in cases with the varying pitch inducer, a significant difference was observed. Model A owned two regimes (Type I and II), while in Model AE there was only the Type I cavity. Referring to the pressure distributions on the blades and the details of the flow field, it seemed that with the alteration of clearance flow characteristics, the cavity feature varied. Clearance flow in the straight casing tended to penetrate far more upstream, resulting in a more apparent axial extension of cavities. As for step casing, the clearance flow tended to interact with the main flow in a closer region with lower clearance velocity and higher leakage flow rate.

**Author Contributions:** Investigation, L.Y.; data curation, L.Y.; methodology, L.Y. and H.Z.; validation, L.Y. and H.Z.; writing—original draft preparation, L.Y.; writing—review and editing, H.C., Z.Z., and S.L.; supervision, H.C., Z.Z., and S.L.; funding acquisition, H.C., Z.Z., and S.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key RD Program of China, grant No. 2018YFB0606103 and the National Basic Research Program of China ("973" Project), grant No. 613321.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Research on the Dynamic Characteristics of Mechanical Seal under Di**ff**erent Extrusion Fault Degrees**

#### **Yin Luo \*, Yakun Fan, Yuejiang Han, Weqi Zhang and Emmanuel Acheaw**

Research Center of Fluid Machinery Engineering and Technology, Jiangsu University, Zhenjiang 212013, China; 2221811033@stmail.ujs.edu.cn (Y.F.); 2221711011@stmail.ujs.edu.cn (Y.H.); 2221811018@stmail.ujs.edu.cn (W.Z.); 5103190336@stmail.ujs.edu.cn (E.A.)

**\*** Correspondence: luoyin@ujs.edu.cn; Tel.: +86-132-2262-6939

Received: 24 July 2020; Accepted: 27 August 2020; Published: 30 August 2020

**Abstract:** In order to explore the dynamic characteristics of the mechanical seal under different fault degrees, this paper selected the upstream pumping mechanical seal as the object of study. The research established the rotating ring-fluid film-stationary ring 3D model, which was built to analyze the fault mechanism. To study extrusion fault mechanism and characteristics, different dynamic parameters were used in the analysis process. Theoretical analysis, numerical simulation, and comparison were conducted to study the relationship between the fault degree and dynamic characteristics. It is the first time to research the dynamic characteristics of mechanical seals in the specific extrusion fault. This paper proved feasibility and effectiveness of the new analysis method. The fluid film thickness and dynamic characteristics could reflect the degree of the extrusion fault. Results show that the fluid film pressure fluctuation tends to be more intensive under the serious extrusion fault condition. The extrusion fault is more likely to occur when the fluid film thickness is too large or too small. Results illustrate the opening force is affected with the fluid film lubrication status and seal extrusion fault degrees. The fluid film stiffness would not always increase with the rotating speed growth. The seal fault would occur with the increasing of rotating speeds, and the leakage growth fluctuations could reflect the fault degree.

**Keywords:** mechanical seal; dynamic characteristics; extrusion fault; numerical simulation; sealing performance; fluent

#### **1. Introduction**

As a widely used fluid transmission device, centrifugal pump plays a crucial role in the national economy. The proportion of pumps equipped with mechanical seals in industry has risen to 80%. The ratio is even greater among petrochemical industry, affecting up to 90%. Besides, statistics reveal that the fault caused by mechanical seals accounts for more than 40% in all machine faults. According to the statistics of centrifugal pump faults from the German Engineering Association [1], the sealing fault ratio is 31%, the rolling bearing fault ratio is 20%, the leakage fault ratio is 10%, the motor fault ratio is 10%, the rotor fault ratio is 9%, and the sliding bearing fault ratio is 8%. The sealing fault can cause an unpredictable waste of resources [2–4]. Mechanical seal faults can affect the internal flow of the centrifugal pump [5–7]. A serious situation would cause casualties and property losses. Thus, the research on mechanical seal fault diagnosis is necessary to ensure reliability and safety [8–10]. Moreover, it is necessary to research the dynamic characteristics detection of each stage when the mechanical seal fault occurs.

Great efforts have been made to do research on the dynamic characteristics of mechanical seals [11–14]. He et al. [15] studied that the viscous shear heat and frictional heat due to asperity contact decrease with an increase in the thickness of the tapered film. As the shaft decelerates, the wear distance rate increases with an increase in the axial stiffness. The axial damping only affects the duration of the oscillations. Zhang et al. [16] researched that lubrication reduces friction and wear and generates heat, but leakage has to be considered. The effects of the sealing surface characteristics on the leakage, and the effects of the external factors of the sealing device on the leak rate. Towsyfyan et al. [17] mainly introduces the fault detection of mechanical seal friction by acoustic emission technology. Zhang et al. [18] investigated the fluidic leak rate through metal sealing surfaces by developing fractal models for the contact process and leakage process. Gropper et al. [19] provides a comparative summary of different modeling techniques for fluid flow, cavitation, and micro-hydrodynamic effects. Migout et al. [20] studied the relationship between the temperature change of the sealing medium and the vaporization of the medium under the seal rings face deformation condition, and pointed out that the stability of the fluid film would be seriously affected when the temperature change gradient of the medium is large. Varney et al. [21] researched the influence of the installation misalignment of the seal rings, established a three-degree-freedom dynamic model of the stationary ring, analyzed the response characteristics formed by the force excitation in all directions, and pointed out that the increase of the excitation intensity would lead to the occurrence of the collision phenomenon of the seal rings. Zhu et al. [22] studied that the sealing performance is enhanced by increasing the spacing of adjacent sealing sheets. The sealing sheets with positive bending angle have less air resistance in the flow path, which would lead to larger leakage. The increase in the number of sealing sheets gives rise to an increase in the generation probability of recirculation zone and vortex, which aggravates the mainstream energy loss. Mosavat et al. [23] researched that the influences of the thermal radiation on the temperature distribution of the mechanical face seal are investigated. Also, the effect of the stretching and shrinking on the thermal performance of the fin with different profiles are comprehensively studied.

Three-dimensional models of single-cone and double-cone were established [24]. A comparison made between the flow velocity, a shear rate and shear stress in single-cone and double-cone zones. This paper revealed the CFD analysis of the flow of a polymeric material inside the double-cone plasticization-homogenization zone of the screw-disc extruder.

Plenty of investigation and discussion have been undertaken on the development of sealing performance, temperature change, and distribution of mechanical seal. At present, there is a lack of monitoring research on mechanical seal failure state. The difference of this paper is to judge the fault degree of mechanical seal according to the dynamic characteristics.

Current research about the dynamic characteristics of mechanical seals mainly focuses on the seal rings modality variation and internal flow field analysis of mechanical seals. However, there is deficiency existing in mechanical seal fault though pressure, opening force, and the leakage. Therefore, the dynamic characteristics and fault mechanism of mechanical seal in different stages need to be studied urgently.

#### **2. Analysis of Sealing Fault Mechanism**

#### *2.1. Normal State*

As shown in Figure 1, mechanical seal is an important part to prevent leakage in centrifugal pumps. Mechanical seal model in normal state is shown in Figure 1. The rotating and stationary rings stick to each other. The function principle of mechanical seal is that a thin fluid film is formed between the rotating ring and the stationary ring. Hydrodynamic effects would be formed because of the thin fluid film. The fluid film with proper thickness can improve the lubrication performance and seal performance. A certain distance between the rotating ring and the stationary ring would be formed because of the hydrodynamic effects. The normal state is that the thickness between rotating ring and stationary ring is proper. At this time, the fluid film could provide enough open force to prevent the direct contact between the rotating ring and stationary ring. The fluid film would form a certain resistance to prevent the medium from leaking out, and make the seal surface lubricated, so as to achieve the better sealing effect.

**Figure 1.** Schemes of mechanical seal model.

#### *2.2. Fault State*

The mechanical seal would often be in fault due to that the friction and wear behavior are universal phenomena in the rotating parts. Seal fault could be divided into many kinds, including seal surface extrusion fault, face temperature resistance fault, lateral load increasing. Many faults are caused by the squeezing of rotating and stationary rings. Figure 2 depicts the fault schematic diagrams, which show the extrusion fault in different degrees. The cause of the extrusion fault on the mechanical seal face might result from the relative deviation of the two seal ring surfaces during operation, the misaligned installation of the mechanical seal, or the mechanical seal external pressure.

Excessively small gap Normal state Excessive clearance

**Figure 2.** Physical pictures of different fault degrees.

#### *2.3. Extrusion Faults*

This paper mainly studied the fault caused by the rotating and stationary ring extrusion of the mechanical seal. The mechanical seal would form a stable fluid film during normal operation. However, the stable fluid film would be destroyed, when the mechanical seal is in fault state. From normal to fault operation state, the internal hydrodynamic pressure and the force of seal rings would change. When the relative movement of the rotating and stationary rings results in extrusion fault, the slight faults would affect the sealing performance and the serious faults would result in extrusion deformation, wear damage, and fault of the rotating and stationary rings.

Figure 3 showed that the characteristics and change of the rotating ring surface under different fault degrees. Figure 3a indicated that the seal surface would be smooth and intact in normal state. Abrasion and damage would be found in the rotating ring surface when the seal is in slight faults, which is showed in Figure 3b. The crack caused by light extrusion could be seen in the red circle of Figure 3c. Significant damage and obvious behaviors would occur when the seal was in severe faults, which are illustrated in Figure 3c. The end face damage caused by severe extrusion fault could be seen in the red circle in Figure 3c.

(**a**) Normal state (**b**) Slight extrusion fault (**c**) Severe extrusion fault

**Figure 3.** Rotating ring in different extrusion fault.

Change in mechanical seal rings was too complicated to be expressed by mathematical equations. Fault physical models were necessary to be carried out to describe the corresponding relationship between the extrusion fault and mechanical seal characteristics.

### **3. Establishment of Calculation Model**

#### *3.1. Fault Physical Models*

In order to study the causes of mechanical seal fault and the internal sealing mechanism, it was necessary to establish a reasonable method and a proper fault physical model for analysis. The actual mechanical seal fault problem was affected by various factors. Therefore, the physical model of mechanical seal fault should be simplified. The fault of mechanical seal was mostly caused by extrusion wear of rotating and stationary rings, which resulted from the long-term uneven stress and the change of relative position of rotating and stationary rings during the long-term operation. In this paper, a simplified physical model of mechanical seal fault was established to analyze the extrusion fault of rotating and stationary rings. To establish the fault physical model, two aspects need to be considered. One is to assume that the relative movement of the rotating and stationary rings only occurs along the axial direction, and the movement along the radial direction was assumed to be zero. The second is to assume that the material texture of the rotating and stationary ring is uniform. In this paper, the thickness of the fluid film between the rotating and stationary rings was used to represent the distance between the rotating and stationary rings when different extrusion faults occur in the axial direction of the mechanical seal. The dynamic characteristics of fluid film with different thickness could reflect the different fault degrees of mechanical seal when extrusion fault occurs. Finally, the fault model of mechanical seal under this fault was formed.

In this paper, the upstream pumping mechanical seal was selected as the research object. The physical model was composed of three main parts: rotating ring, fluid film, and stationary ring. Figure 4 showed the axial section diagram of the rotating ring modeling and the geometric parameters of the section. The main geometrical parameters of the fluid film were selected: the inner radius of the rotating rings *r*<sup>i</sup> = 25 mm, the outer radius of the rotating ring *r*<sup>o</sup> = 32 mm. Figure 5 indicated the axial section diagram of the stationary ring modeling and the geometric parameters of the section. The inner radius of the stationary ring *r*<sup>i</sup> = 24.5 mm, the outer radius of the rotating ring *r*<sup>o</sup> = 31 mm. Figure 6 showed the three-dimensional diagram of fluid film modeling. Besides, the dimensions of the inner and outer rings of the fluid film were shown in detail. Because the order of magnitude in the direction of film thickness was micrometer, Figure 6 showed the enlarged fluid film thickness. The thickness of the fluid film was set to 1, 2, 3, 4, 5, 6, 7, 8, 9 µm, respectively. Figure 7 illustrated the computational domains of the mechanical seals. μ μ μ

Axial section diagram Geometric parameters

**Figure 5.** Stationary ring modeling diagram.

**Figure 6.** 3D diagrams of fluid film modeling.

**Figure 7.** Computational domains of the mechanical seal (1-stationary ring, 2-fluid film, 3-rotating ring).

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation as well as the experimental conclusions that can be drawn.

### *3.2. Dynamic Calculation Model*

The fault of the mechanical seal was mostly manifested in the squeeze and wear fault of the rotating and stationary rings. The dynamic calculation model was selected according to the simplified fault physical model given above when the rotating and stationary rings fail due to extrusion fault. With the deepening of theoretical research on mechanical seals, the mechanical seal faults include various sciences such as mechanics, power, fluids, materials, chemistry, heat transfer, and tribology. In order to study the law between the mechanism and structural characteristics of mechanical seal when extrusion fault occurs, a dynamic calculation model suitable for fault physical model was adopted in this paper.

In order to make the simulation simple and clear, the following assumptions were made in this paper, considering the existence of fluid pressure, elastic force, solid deformation and the interaction force and heat transfer between fluid and solid in the mechanical seal.


#### 3.2.1. Control Equations

The flow and diffusion of the liquid film fluid inside the mechanical seal satisfies the momentum equation, energy equation, and continuity equation [25].

1. Fluid Domain Equation

The mass conservation equation

$$\frac{\partial \rho\_{\rm f}}{\partial t} + \frac{\partial}{\partial x\_j} (\rho\_{\rm f} \nu\_i) = 0 \tag{1}$$

In the above formula: ρ<sup>f</sup> is the fluid density; ν is the fluid velocity, subscript *i*, *j* = 1, 2, 3, representing the components in three directions, *t* is the time.

The momentum equation

$$\frac{\partial(\rho\_l \nu\_i)}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j}(\rho\_l \nu\_i \nu\_j) = -\frac{\partial p}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j}(\mu \bullet \frac{\partial \nu\_i}{\partial \mathbf{x}\_j}) \tag{2}$$

In the above formula: *p* is the pressure of fluid film, µ is the dynamic viscosity of the fluid. The energy equation

$$\frac{\partial \rho \mathbf{E}}{\partial t} + \nabla \bullet [\nu (\rho \mathbf{E} + p)] = \nabla \bullet \left[ k\_{eff} \nabla T - \sum\_{j} h\_{j} I\_{j} + (\Gamma\_{eff} \bullet \nu) \right] + S\_{h} \tag{3}$$

In the above formula: E is the total energy of the fluid micelle, ρ is the density of the fluid micelle, and *p* is the pressure of fluid film. Γ*e*ff is the effective stress of fluid domain, *h<sup>j</sup>* is the enthalpy of the component, *Ke*ff is the effective thermal conductivity of fluid film, *J<sup>j</sup>* is the diffusion flux of the component *j*, *S<sup>h</sup>* is the source term including other volumetric heat.

2. Solid Domain Equation

$$\stackrel{\bullet}{M}\_{\rm s} \stackrel{\bullet}{d} + \stackrel{\bullet}{C}\_{\rm s} \stackrel{\bullet}{d} + \stackrel{}{K}\_{\rm s} \stackrel{d}{d} + \stackrel{\bullet}{\tau\_{\rm s}} = 0$$

In the equation, *M*<sup>s</sup> represents for the mass matrix of the solid element, *C<sup>s</sup>* represents for the damping matrix of the solid element, *K<sup>s</sup>* represents for the rigidity matrix of the solid element, *d* represents for the displacement vector of the solid element, τ*<sup>s</sup>* represents for the stress on the rotating ring and stationary ring.

At the same time, the thermal deformation term caused by the temperature difference in the solid area as followed.

$$f\_T = \alpha\_T \bullet \mathsf{VT} \tag{5}$$

In the formula: α*<sup>T</sup>* is the coefficient of thermal expansion related to temperature. *T* is the temperature of the seal rings.

#### 3. Dynamic Model of Axial Movement

Based on the analysis of the kinematic relationship, the dynamical model of mechanical seal system was derived by using the D'Alembert principle. The force and axial movement of the mechanical seal are shown in Figure 8.

**Figure 8.** Axial dynamic model of mechanical seal.

From the dynamic point of view, the dynamic equation of mechanical seal axial movement was established as follows [26]. .. .

$$
\mathbf{\dot{m}}\ddot{\mathbf{Z}}\_{\mathbf{s}} + \mathbf{c}\dot{\mathbf{Z}}\_{\mathbf{s}} + \mathbf{F}\_{\mathbf{s}} = \mathbf{F}\_{\mathbf{l}} + \mathbf{F}\_{\mathbf{c}} \tag{6}
$$

 In the formula: *m* is the equivalent mass of the rotating ring and the spring. *Z<sup>s</sup>* is the axial relative displacement between the rotating ring and the stationary ring. *c* is the axial damping coefficient of spring and auxiliary seal rings. *F<sup>s</sup>* is the force of the spring on the rotating ring. *F<sup>l</sup>* is the axial force of the fluid film between the end faces of the rotating ring and stationary ring, and *F<sup>c</sup>* is the axial contact force between the rotating ring and stationary ring.

#### 3.2.2. Sealing Performance Parameters

#### 1. Seal Opening Force

The opening force of mechanical seal surface is the sum of the pressure exerted on the seal surface by the liquid film fluid. The opening force could be obtained by integrating the pressure field of the liquid film on the seal surface.

$$F\_0^{\Omega} = \int \int \mathbf{p} dA = \int \int \mathbf{p} r dr d\theta \tag{7}$$

<sup>0</sup> =p p In the formula: *p* is the pressure of fluid film, *r* is the radial coordinate, and Ω is the whole calculation area.

#### 2. Leakage of Mechanical Seal

Leakage is an important indicator for measuring the performance of mechanical seal. Leakage *Q* could be synthesized by this formula [27].

$$Q = \frac{\pi d\_m h\_o^3 \Delta p}{12\mu b} \tag{8}$$

Ω

12 *Q*-leakage, *dm*-average diameter of the sealing surface, *ho*-the thickness of fluid film, ∆*p*-pressure difference, µ-dynamic viscosity of the medium, *b*-effective width of the seal.

#### *Δ μ* **4. Dynamic Simulations and Analysis**

In this paper, the computational domains were meshed using the structured blocking hexahedral method. Due to the large difference between the radial and axial length of the fluid film parameters, the diameter of the fluid film is millimeter, but the thickness of the fluid film is micron, thus it is difficult to directly divide the three-dimensional mesh of the liquid film. Therefore, the fluid film corresponding to the center angle of 8.5 degrees was first meshed in this paper, and then the entire fluid film is formed based on this array. Figure 9 revealed coupling relationship among multiple fields in the dynamic simulations. To ensure the computational accuracy, mesh independence analysis was

conducted. Calculation of fluid film parameters with different number of grids was shown in the Table 1. The velocity was not growing when the number of grids increased to 3.6 million. As shown in Figure 10, the pressure became relatively stable after increasing the grid number to 3.6 million elements. After the grid independence test, the final grid number unit was 362,943. The fluid membrane grids were shown in Figure 11. In order to remove the influence of the order of magnitude of radial and axial fluid film grids on the simulation results, a grid-independent assessment was conducted. The numerical results were obtained by dividing the fluid film into different mesh numbers and repeating the simulation. The error of the final results is less than 1%, which shows that the mesh division has no effect on the simulation results.

**Figure 9.** Multi-field coupling.



**Figure 10.** Mesh independence.

**Figure 11.** Mesh of fluid film and mesh quality.

Mechanical seal material and boundary condition were set as followed:

The rotating ring material used in this paper is silicon carbide, and the stationary ring material is carbon graphite. The main properties of the materials used in the rotating and stationary rings are shown in Table 2.


The details about boundary, mesh, and calculation are listed in Table 3.


**Table 3.** Details about boundary, mesh and calculation.

When the residual values of all variables are reduced to 10−<sup>3</sup> , the calculation converged.

Entrance boundary conditions in numerical simulation: the pressure inlet boundary condition was adopted, and the inlet boundary position was set outside the fluid film.

Exit boundary conditions in numerical simulation: the pressure outlet boundary condition was adopted, and the outlet boundary position was set inside the fluid film.

Wall boundary conditions in numerical simulation: Standard wall functions were used.

It had been a hot issue that the fluid film flow is laminar flow or turbulent flow in the mechanical seal field. The fluid coefficient α method was adopted to determine whether it is laminar or turbulent [28] in this paper.

$$\alpha = \sqrt{\frac{\text{Re}\_{\text{c}}}{1600}^2 + \left(\frac{\text{Re}\_p}{1600}\right)^2} \tag{9}$$

α < 0.58, is laminar flow. α > 1, is turbulent flow.

$$\text{Re}\_{\text{c}} = \frac{2\pi\rho vrh}{60\mu} \tag{10}$$

$$\text{Re}\_p = \frac{\rho Q}{2\pi\mu r} \tag{11}$$

*<sup>Q</sup>* <sup>=</sup> 1.12 <sup>×</sup> <sup>10</sup>−<sup>5</sup> kg/s, *<sup>n</sup>* <sup>=</sup> 1500~6000 rpm, *<sup>r</sup>* <sup>=</sup> 0.031 m, *<sup>h</sup>* <sup>=</sup> 1~9 <sup>µ</sup>m, <sup>ρ</sup> <sup>=</sup> 998 kg/m<sup>3</sup> .

Through calculation could get, α = 0.29 < 0.58, the fluid film is laminar flow. Therefore, the laminar model was adopted for the relevant flow dynamic calculation.

The simulation part was set according to the theoretical physical fault model and the dynamic model set up in the Section 3.2.1., and the flow chart of simulation and the analysis details are shown in Figure 12. Simple C algorithm and steady-state solver were used during the calculation of the fluid domain. Firstly, the fluid domain results were calculated by CFD. Then, the fluid domain results were loaded to the solid domain through the Workbench platform. The dynamic equation of mechanical seal axial movement was used during the calculation of the deformation and forces in solid domain. Thus, the solid domain dynamic results, deformation, and forces were obtained. Finally, the fluid film dynamic characteristics, the force, and deformation of the seal rings could be analyzed from the results.

**Figure 12.** Flow chart of calculation.

The first group of simulation process is for different film thickness. The force and deformation simulation data of rotating ring and stationary ring were obtained, and the pressure of fluid film could be got at the same time.

The second group in the simulation process was carried out under different speeds. The hydrodynamic and structural dynamic characteristics were obtained.

μ

Repeat each simulation process. Find the characteristics under different conditions through data analysis. μ

μ

μ μ μ

μ

#### **5. Results and Discussion**

#### *5.1. Fluid Film Pressure Analysis*

The pressure distributions under different extrusion degrees are shown in Figure 13. The fluid film pressure would grow as the fluid film thickness increases when the rotating speed *n* = 2950 r/min, pressure inlet was 0.2 MPa. The fluid film thickness represents the extrusion fault degree. The dynamic characteristic parameters could be reflected by three typical cases. Two-µm fluid film could reflect the serious extrusion fault, 3-µm fluid film was a sign of slight extrusion fault, 4-µm-fluid film was on behalf of the normal state. Such large pressure change between 2-µm fluid film and 3-µm fluid film means that the fluid film thickness has strong recognition capability for the changes of the mechanical seal operation. When the extrusion fault of mechanical seal becomes severer, the fluid film pressure fluctuation would tend to be more intensive. Such tendency indicated that the thickness and the pressure could be treated as the indicator of extrusion fault. μ μ μ μ

**Figure 13.** The pressure of 2, 3, 4 and 5 µμm fluid film.

No significant fluctuation could be found in the 4-µm fluid film pressure contours. This was mainly due to that the 4-µm fluid film was the relatively proper thickness for this mechanical seal. From Figure 13, it could be seen that the pressure was relatively evenly distributed at this time. The thicker fluid film leaded to the bigger viscous shear flow. In addition, pressure bearing capacity could also become stronger with the rapidly increasing thickness. The results indicated that the thicker fluid film could efficiently be used to maintain the hydrodynamic effects and improve the lubrication performance. However, there was deficiency existing in this parameter, which was not sensitive to the light extrusion fault.

The pressure fluctuation under different fluid film is shown in Figure 14. The pressure mainly fluctuates at 0.189 MPa when the thickness of the fluid film is 2 µm. The pressure mainly fluctuates at 0.192 MPa when the thickness of the fluid film is 3 µm. The pressure mainly fluctuates at 0.195 MPa when the thickness of the fluid film is 4 µm. The pressure mainly fluctuates at 0.198 MPa when the thickness of the fluid film is 5 µm. The pressure fluctuation is largely affected by the operation of the rotating ring and stationary ring. When there is extrusion fault, the operation of the rotating ring and stationary ring will change. The pressure distribution and pressure fluctuation of the fluid film will change with the fault. To a certain extent, the pressure fluctuation could reflect the occurrence of extrusion fault. Besides, the fault degree could be preliminarily judged according to the pressure distribution and pressure fluctuation within a certain range.

**Figure 14.** Pressure fluctuation diagram under different fluid films.

#### *5.2. Seal Performance Analysis*

μ μ μ Compared with pressure, leakage *Q* showed a stronger ability in reflecting the slight extrusion fault. That was mainly due to the leakage as the main indicator that could measure seal fault. Figure 15 showed that the leakage increases with fluctuations when the fluid film becomes thicker. Besides, the minimum value of the leakage, shown in Figure 15, under different fluid film thickness. As shown in the red circle in Figure 15a, too thin fluid film could lead to the severe wear. The leakage curve reflected that the leakage increasing rate is the slowest when the fluid film thickness is between 5 and 6 µm. It depicted that this fluid film thickness was the most beneficial to mechanical seal. So, the best thickness area for this mechanical seal was from 5 to 6 µm. However, the increasing tendency is obvious when the fluid film thickness was more than 6 µm in the red circle in Figure 15b. It could be concluded that faults have occurred on the mechanical seal. With the development of faults, the leakage curve would still increase. According to Figure 15, a large amount of leakage resulted from too thick fluid film, which was directly related with the extrusion fault. During the simulation operation of this mechanical seal, the extrusion fault was more likely to occur when the fluid film thickness was too large or too small. Dynamic characteristics of fluid film would change when the extrusion fault occurs. Hydrodynamic effects were the typical parameters to measure the sealing performance. As the fluid film thickness beyond the best thickness area marked with red circle in Figure 15c, the hydrodynamic effects of fundamental frequency became weaker, while volume force and inertia force got larger.

**Figure 15.** Leakage under different fluid film.

μ μ μ The thickness of fluid film was set from 1 to 9 µm. The opening force changed from 60 N to 87 N. The rotating speed was set as 2950 r/min. Figure 16 illustrated that with the increasing of fluid film thickness, opening force became larger rapidly. Moreover, the opening force would have a fluctuating downward trend as the fluid film thickness increases. The opening force could be regarded as the important parameter for the seal dynamic characteristics. As shown in the Figure 16, the opening force was small when the fluid film thickness was less than 3 µm. It was because that serve extrusion fault occurs when the thickness was too thin. The mechanical seal surface friction would happen directly when the thickness was too small. Thus, there would be less or no hydrodynamic effect in this situation. A downward trend of opening force was shown in curve when the fluid film thickness was more than 4 µm. Moreover, the opening force was smaller while the thickness is larger. It indicates that serious damage has occurred on the surface of the rotating ring and stationary ring with the increasing of the thickness. The opening force could have relationship with the fluid film lubrication status and seal extrusion fault degree.

Besides, cavitation could be found in Figure 16b,c. With the increase of fluid film thickness, cavitation phenomenon is further strengthened. It is also due to the cavitation phenomenon that the opening force increases firstly and then decreases with the increase of film thickness. The pressure reduction of mechanical seal is mainly due to the local separation caused by narrow gap, micro groove machined on the surface and surface roughness, and vortex caused by micro modeling. In the conventional scale flow, the surface roughness of micro channel, which is often neglected due to its small influence, has a significant influence in the micro channel flow. The micro disturbance caused by the surface roughness often affects the flow at the edge of the fluid film, which is also one of the main reasons for the cavitation of the micro gap fluid film in the mechanical seal.

**Figure 16.** Opening force under different fluid film.

As shown in Figure 17, with the increasing of rotating speed, the leakage grew with many fluctuations. The Figure 17a could reflect the same trend of leakage under different film thickness, and Figure 17b could reflect the difference under each speed. It could be seen that the most obvious growth occurs when the rotating speed *n* from 2500 r/min to 3500 r/min. This was mainly due to that the mechanical seal was from normal operation state to the fault status. The inherent factors could be the opening force and friction torque, which resulted in the rapid growth of leakage. Thus, the leakage could be related to the rotating speed. The fluid film stiffness would not always increase with the rotating speed growth, because the growth of the friction torque was suppressed to a certain extent. Based on the combined action of the hydrodynamic and fluid film stiffness, the seal leakage could reflect the sealing performance and the seal faults.

(**a**) Leakage under different film thickness

(**b**) Leakage at different speeds

**Figure 17.** Leakage under different rotating speeds and fluid film thickness.

#### **6. Conclusions**

The dynamic characteristics of mechanical seal under different fault conditions were processed and analyzed in this research. Several conclusions could be drawn from the results described above.


the hydrodynamic and fluid film stiffness, the seal leakage could reflect the sealing performance and the seal fault. There must be a law between the extrusion degrees and fluid film thickness. This paper researched the law that too thin or too thick fluid film would result in the heavy extrusion fault. The fluid film stiffness, leakage, and opening force are the important parameters, which have recognition capability for the extrusion degrees of the mechanical.

This research work has proved the feasibility that the dynamic characteristics of mechanical seal could be found to reflect the degree of extrusion fault. Besides, the research conclusions could have the reference value of mechanical seal fault detection. Further research should focus on the different type of the mechanical seal and find the accurate correspondence relationship.

**Author Contributions:** Conceptualization, Y.L.; methodology and software, Y.L. and Y.F.; validation, Y.L. and Y.F.; formal analysis, Y.L. and Y.F.; investigation, Y.F. and W.Z.; data curation, Y.L. and Y.F.; writing—original draft preparation, Y.L. and Y.F.; writing—review & editing, Y.L., Y.H. and E.A.; project administration, Y.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 51979127, the Nature Science Foundation of Jiangsu Providence, grant number BK20171403, and German KSB Global Headquarters Research Fund, project number: 1.2018.07.1

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Flow and Noise Characteristics of Centrifugal Fan in Low Pressure Environment**

### **Xilong Zhang, Yongliang Zhang \* and Chenggang Lu**

School of Mechanical and Automotive Engineering, Qingdao University of Technology, Qingdao 266520, China; zhangxilong@qut.edu.cn (X.Z.); papqut@yeah.net (C.L.)

**\*** Correspondence: zhangyongliang@qut.edu.cn

Received: 16 July 2020; Accepted: 12 August 2020; Published: 13 August 2020

**Abstract:** The influence of low-pressure environment on centrifugal fan's flow and noise characteristics was studied experimentally and numerically. A testbed was established to conduct the experimental test on the performance of a centrifugal fan, and the characteristic curve and power consumption curve of the fan under different pressure were obtained. Then the simulation model of the centrifugal fan was established, which was used to simulate the working process of centrifugal fan under different negative pressures. The results showed that the total pressure and static pressure of the fan decrease with the decrease of the ambient pressure. The total and static pressures of the fan under 60 kPa pressure condition decreased by 42.3% and 38.3%, respectively, compared with those of fan under the normal pressure. The main reason for this phenomenon is that the decrease of the environmental pressure leads to the decrease of air density. Besides, with the drop of environmental pressure, the sound pressure and sound power of the fan noise decreases.

**Keywords:** centrifugal fan; noise characteristics; power consumption; negative pressure; sound pressure

### **1. Introduction**

Under the plateau environment, the air is thin and the air pressure is low, which changes the physical properties and flow characteristics of the air. For every 1000 m increment in altitude, the atmospheric pressure drops by about 10 kPa, and the air density also gradually decreases. When the altitude is 5000 m, the air density is 0.7263 kg/m<sup>3</sup> , which is only about half of that in the plain area [1,2]. However, from the sea level to the elevation below 85,000 m, the volume ratio of the main gases such as nitrogen and oxygen is basically the same at each altitude. So, the relative molecular mass of air remains basically unchanged. Density is proportional to atmospheric pressure at a given temperature [3,4]. When the temperature is constant, molecular concentration and air density increase with an increase of pressure [5]. The characteristics of air flow and noise change correspondingly in the application of vehicle fan and air conditioning fan. The flow rate, static pressure, axial power, efficiency, rotational speed, noise, and other performance parameters of the fan are all related to the physical properties of the air, so the flow and noise characteristics of the fan will inevitably be affected by the change of environmental pressure.

At present, many scholars conducted deep research on the flow and noise characteristics of the fan. Lee et al. [6] studied the effect of blade inclination angle, blade thickness, and maximum blade thickness location of the low-speed axial fan on fan efficiency by response surface method. It was concluded that the blade inclination angle had the greatest influence on the fan efficiency. Gholamian et al. [7] used the method of CFD to study the effect of inlet diameter on fan efficiency and flow field. It was found that the size of the inlet diameter had a great influence on the fan performance and efficiency. When the inlet diameter differs by 2 cm, the fan performance and fan efficiency change by more than 14%. In order to improve the performance of the fan, a variety of advanced technologies have been

put forward and a great deal of research has been carried out. For example, fan blade bending sweep technology [8] and blade twisting technology [9] are the most commonly used techniques.

However, the studies on the performance of fans under low pressure environment are still relatively few. The existing studies do not take the impact of environmental pressure on fan efficiency into consideration, which lacks experimental verification. So, it is insufficient to provide theoretical support for the application of fans in low-pressure environment at plateau, and further research is urgently needed.

In the noise characteristics of fans, through theoretical analysis, it was found by Sharland [10] that the aerodynamic noise source of axial fan belonged to dipole source, which was closely related to the pressure pulsation of blade acting on the air passing through the fan. Li [11] regarded the air as an incompressible fluid to calculate the fan performance, which saved the calculation time, and the error was within the allowable scope of the project. By analyzing the influence of metal stamping fan blade thickness on fan performance, it can be known that the greater the blade thickness is, the greater the noise will be generated. Hodgson et al. [12] concluded that the magnitude of fan noise was positively correlated with driving voltage and negatively correlated with outlet flow through the computer cooling fan noise test.

From the existing literature, the studies on the fan were mainly based on one atmosphere. However, the articles that studied the noise characteristics of the fan in the low-pressure environment were not retrieved. So far, it is difficult to reveal the mechanism of low-pressure environment on the flow and noise characteristics of fans. Furthermore, it is difficult to meet the design and calculation requirements of fans in low pressure environment.

#### **2. Analysis of the Influence of Low-Pressure Environment on Noise Characteristics of Centrifugal Fan**

#### *2.1. Simplification of Air Flow in Centrifugal Fan*

In order to study air flow in a centrifugal fan, the flow in the fan is properly simplified as follows:


#### *2.2. Noise Analysis of Centrifugal Fan in Low Pressure Environment*

The noise of the fan blade is mainly determined by the Lighthill fundamental equation [13]. Namely, the sound power (*W*) dominates, which has the following relationship with air flow rate (*u<sup>f</sup>* ):

$$
\rho \mathcal{W} \sim \rho L^2 \frac{\mu\_f^6}{c^3},
\tag{1}
$$

where ρ is air density, *L* is the length along the direction of the flow, and *u<sup>f</sup>* is average flow rate (*u<sup>f</sup>* ).

The formula for calculating the propagation velocity (*c*) of sound in air in the above formula is *c* = (*k*/ρ) 0.5 [14], in which *k* is the volumetric modulus of elasticity of the medium.

According to the theory of fluid molecules, the specific heat of air can be regarded as a constant value in the pressure change range of 0~1 atm, so the *k* value does not change with the ambient pressure. According to Equation (1), when the air temperature is the same, the sound propagation velocity *c* at different pressures can be approximately considered to remain unchanged. Therefore, according to Equation (1), the sound power *W* generated by the fan has a linear relationship with the air density ρ, namely *W*~ρ.

The relationship between sound pressure (*P*) and sound power (*W*) [15] is *P* = p *W*ρ*c*/*A*, in which *A* means effective flow area. So, the sound pressure ratio at different ambient pressures is

$$\frac{P\_H}{P\_0} = \sqrt{\frac{\mathcal{W}\_H \rho\_H}{\mathcal{W}\_0 \rho\_0}},\tag{2}$$

where the subscript *H* and *0* mean that the height above sea level are *H* m and 0 m, respectively.

The above equation can be further simplified as *PH*/*P*<sup>0</sup> = ρ*H*/ρ0, in which it can be seen that the sound pressure *P* is also in a linear relationship with the air density.

From the above analysis, it can be seen that the sound power (*W*) and sound pressure (*P*) of the fan are in a linear relationship with the air density (ρ). While the air density is greatly affected by the environmental pressure, so the noise of the fan is closely related to the environmental pressure.

#### **3. Experimental Study**

#### *3.1. Fan Testbed*

The performance test of centrifugal fan is conducted on a fan testbed, which is composed of air loop, fan, and control unit and auxiliary equipment. The fan testbed is shown in Figure 1. The real experimental device is shown in Figure 2.

**Figure 1.** Centrifugal fan testbed.

**Figure 2.** The real experimental device.

During the experiment, the method of gradually loading a uniform circular plug at the throttling metal mesh is used to simulate the resistance and to form a vacuum in the wind tunnel. The motor's electric parameters of centrifugal fan are measured by using a motor economic operation instrument. The pressures in Point 1 and Point 5 in Figure 1 are measured with a high-precision manometer, which are converted to obtain the flow rate of air [16]. The ambient temperature, humidity, and atmospheric pressure are measured with a thermometer, hygrometer, and atmospheric pressure gauge, respectively.

#### *3.2. Uncertainty in Experiments*

The experimental sample is a centrifugal fan, which is shown in Figure 3. The thickness of the blade is 3 mm, the outlet mounting angle is 60◦ , the outer diameter is 600 mm, and the number of blades is 13.

**Figure 3.** Experimental sample of centrifugal fan.

The rated rotational speed and maximum rotational speed of this centrifugal fan are 5000 r/min and 5994 r/min, respectively. The rated power and maximum power are 78.4 kW and 92 kW, respectively. The rated air volume and rated static pressure are 7.6 m<sup>3</sup> /s and 5427 Pa, respectively.

The uncertainty of the fan testbed is calculated using the method in the literature [17], which is described as:

$$\frac{\delta\_R}{R} = \sqrt{\left(\frac{\partial \mathcal{R}}{\partial v\_1}\right)^2 \left(\frac{\delta v\_1}{v\_1}\right)^2 + \left(\frac{\partial \mathcal{R}}{\partial v\_2}\right)^2 \left(\frac{\delta v\_2}{v\_2}\right)^2 + \dots + \left(\frac{\partial \mathcal{R}}{\partial v\_n}\right)^2 \left(\frac{\delta v\_n}{v\_n}\right)^2} \tag{3}$$

where *R* is a function of the variables *v<sup>i</sup>* (*i* = 1,2, . . . , n), and δ*v<sup>i</sup>* (*i* = 1, 2, . . . *n*) represents the uncertainties of variables *v<sup>i</sup>* . After calculation, the uncertainties of mass flow rate is 0.212%, the fan pressure head is 0.257%, the fan efficiency is 0.97%, and the average uncertainties of the main parameters are less than 1%. The relationship between *R* and *v<sup>i</sup>* is shown in Table 1. In the table, ε is blade displacement coefficient; *D* is blade diameter; *b* is width of blade; *v*<sup>r</sup> is the radial component of the absolute velocity; ρ is air density; *D*<sup>o</sup> is diameter of outlet; *D*<sup>i</sup> is diameter of inlet; *S*<sup>r</sup> is radial blade clearance; *g* is gravitational acceleration; *u* is peripheral speed; *v*<sup>u</sup> is the axial component of the absolute velocity. Therefore, it can be considered that the testbed can meet the requirements of fan test accuracy.

**Table 1.** The relationship between *R* and *v<sup>i</sup>* .


#### *3.3. Experimental Results and Analysis*

As the fan rotates, the air can be considered as an incompressible flow medium; that is, the density of air can be considered as constant. So, the flow characteristics of the fan can be characterized by the volume of medium passing through the fan per unit time and is denoted by *qv*.

As the air passes through the fan, the increase of pressure is the total pressure of the fan and is denoted as *PtF*, i.e.,:

$$P\_{\rm tF} = P\_2 - P\_1 + \rho \frac{u\_{f2}^2 - u\_{f1}^2}{2},\tag{4}$$

where *P*<sup>1</sup> and *P*<sup>2</sup> are the static pressure at the inlet and outlet; *u*f1 and *u*f2 are the average flow rate at the inlet and outlet, and *r* is the rotating radius of the blade.

The increased total power *P<sup>f</sup>* of the medium of the fan in unit time can be expressed as

$$P\_f = q\_\upsilon P\_{t\mathcal{F}\prime} \tag{5}$$

During the experiment, the power consumption curve of the fan is measured by adjusting the speed and volume flow rate of the fan, which is shown in Figure 4.

**Figure 4.** (**a**) Power consumption curve of centrifugal fan; (**b**) Performance characteristic curve of centrifugal fan.

In Figure 4, a total of 6 speed test points are selected in the experiment around the rated speed of 5000 r/min, which are 3000 r/min, 3500 r/min, 4000 r/min, 4500 r/min, 5000 r/min and 5500 r/min, respectively. The relationship of volume flow rate, rotational speed, power, and total pressure of the fan is shown in Figure 4b. During the whole test process, the volume flow rate of the fan varies from 3.5 m<sup>3</sup> /s to 9.4 m<sup>3</sup> /s, and the fan average power varies from 16 kW to 92.8 kW. It can be found that with the same speed of fan and with the decrease of volume flow rate, the total pressure of the fan and the power consumed by the fan increase gradually.

#### **4. Simulation Study**

#### *4.1. Grid Division and Definition of the Boundary*

The finite volume analysis software Ansys Fluent [18] is used to simulate the centrifugal fan, and the flow field model of the fan is shown in Figure 5. The flow field model of fan is divided into three parts, which include inlet extension area, rotating fluid area, and outlet extension area. Besides, a multi-reference frame (MRF) approach is adopted for the rotating fluid area. The medium flowing in the fan channel is air, and the fluid in the fan moving area belongs to turbulence flow. The effect of gravity on the flow field is ignored. The standard *k*-ε model is used. In the simulation process, only the ambient pressure and its corresponding air physical properties are changed, while other boundary conditions are not changed. Given the complexity of the fan simulation model, a mixture of structured and unstructured grids is used to partition the fluid region. In the inlet extension and outlet extension areas, structured grid is adopted [19].

**Figure 5.** Flow field model and grid division of centrifugal fan.

In the rotating fluid areas, the tetrahedral grid is used. The near-wall grids are locally encrypted. To assess the influence of the number of grids on the accuracy of the calculation, the grid independence is calculated. Five grid systems, 0.44 million, 0.72 million, 1 million,1.28 million, and 1.56 million, are tested. It is found that the relative errors of total pressure of fan between the solutions of 1.28 million and 1.56 million are less than 0.07%. It can be considered that the simulation calculation result is independent of the number of grids when the number of grids is encrypted to 1.28 million [20]. The inlet is set as "Velocity inlet", the outlet is set as free "Outflow", and the rotating fluid area is set as "no-slip boundary condition".

#### *4.2. Governing Equations*

The governing equations mainly include mass conservation, momentum conservation, and energy conservation equations, which are

1. Mass conservation equation

$$\frac{\partial \rho}{\partial \tau} + \nabla \cdot (\rho \mathcal{U}) = 0,\tag{6}$$

where *U* is the velocity vector, ρ is the density, and τ is time.

2. Momentum conservation equation

$$\begin{array}{l}\frac{\partial(\rho u\_{1})}{\partial \tau} + div(\rho u\_{1} \cdot \mathcal{U}) = div(\mu grad u\_{1}) - \frac{\partial p}{\partial x} + \mathcal{S}\_{u} \\\frac{\partial(\rho u\_{2})}{\partial \tau} + div(\rho u\_{2} \cdot \mathcal{U}) = div(\mu grad u\_{2}) - \frac{\partial p}{\partial y} + \mathcal{S}\_{v} \\\frac{\partial(\rho u\_{3})}{\partial \tau} + div(\rho u\_{3} \cdot \mathcal{U}) = div(\mu grad u\_{3}) - \frac{\partial p}{\partial z} + \mathcal{S}\_{w} \end{array} \tag{7}$$

where µ is kinetic viscosity, *p* is pressure, *u*1, *u*2, and *u*<sup>3</sup> are the components of the velocity vector in the *X*, *Y*, and *Z* directions, and *Su*, *Sv*, *S<sup>w</sup>* are the generalized source terms.

#### 3. Energy conservation equation

$$\frac{\partial(\rho t)}{\partial t} + div(\rho Ult) = div(\frac{\lambda}{c\_p} \text{grad} T) - \frac{\partial p}{\partial x} + \text{S}\_{\text{T}} \tag{8}$$

where *T* is the temperature, λ is the heat transfer coefficient of the fluid, *c<sup>p</sup>* is the specific heat capacity, and *S<sup>T</sup>* is the viscous dissipation term.

The solver uses the SEGREGATED separate implicit solver. The turbulence energy, turbulence dissipation term, and momentum conservation equations are all discretized using the second-order upwind scheme. The governing equations are solved using the transient-SIMPLE method. Second order upwind scheme is chosen to discretize these governing equations. Due to the strong nonlinear relationship between the variables, the iterative solution uses subrelaxation factors. The inlet and outlet turbulence are all set to 0.5%, and the detection surface of mass flow rate is set at the fan outlet section. For the mass conservation equation, when the results of two adjacent calculations are less than 10−<sup>6</sup> , the numerical simulation results can be considered to converge.

#### *4.3. Simulation Model Verification*

The diagram of the comparison between the simulation calculation results and the experimental results at the rotational speed of 5500 r/min under normal pressure is shown in Figure 6.

**Figure 6.** Comparison between the experimental and the simulation results of total pressure of the fan.

It can be seen from the figure that the error between the experimental value and the calculated value is within 6.8% in the whole experiment process, which indicates that the simulation model is reliable. The deviation between the simulation and the experimental results may be caused to some extent by the simplification of the fan model, measurement errors in the experimental process, and errors inherent in the simulation calculation method. At the same time, the simulation model is mainly used for comparative study. Therefore, the simulation model can be used to simulate the fan performance under different environmental pressures.

#### *4.4. Numerical Calculation Results and Analysis*

#### 4.4.1. Pressure Analysis of the Fan under Different Ambient Pressures

Using this model, the performance of the fan under the same rotational speed (5500 r/min) and different ambient pressure is simulated, and the calculation results are shown in Figure 7.

**Figure 7.** Total pressure of the fan under different ambient pressures.

It can be seen from Figure 7 that the total pressure of the fan decreases with the decrease of flow rate, and the lower ambient pressure indicates that the total pressure of the fan is lower under the same condition. The quantitative analysis shows that compared with the normal pressure condition, the total pressure drop amplitude of the fan is basically the same at each working point under low pressure condition. Taking 60 kPa pressure operating point as an example, the average decrease range of total pressure of the fan is 42.3% compared with the normal pressure condition.

Figure 8 shows the changing curve of the static pressure versus flow rate of the fan under different ambient pressures. It can be seen from the figure that the variation trend of the fan static pressure and the total pressure is similar. The lower the ambient pressure, the lower the fan static pressure. Compared with the normal pressure condition, the static pressure drop extent of the fan is basically the same at each working point under the low-pressure condition. Taking 60 kPa pressure operating point as an example, the average decrease extent of fan static pressure is 38.3% compared with normal pressure operating point. At the same time, the air density at 60 kPa is about 39.5% lower than that at atmospheric pressure. Thus, the decline in the total pressure and static pressure performance of the fan is mainly due to the drop in the air density caused by the reduction of the environmental pressure.

**Figure 8.** Static pressure of the fan under different ambient pressure.

Figure 9 shows the change curve of fan outlet pressure head along with flow rate (*u<sup>f</sup>* ) under different ambient pressure. The figure shows that the outlet pressure head of the fan increases with the increase of flow rate (*u<sup>f</sup>* ) but changes little under different ambient pressure. It can be seen that the environmental pressure has little effect on the outlet pressure head of the fan. Although this change is not obvious, it can be found that the outlet pressure head decreases gradually with the increase of ambient pressure.

**Figure 9.** Fan outlet pressure head under different ambient pressure.

#### 4.4.2. Noise Analysis under Different Ambient Pressure

Fan noise propagation is an unsteady process. Large Eddy Simulation (LES) model is used to simulate the sound field, in which the time step is 0.01 ms and each time step is iterated 20 times. At the same time, two noise monitoring points (as shown in Figure 10) are set along the fan radial direction and axis central direction respectively for real-time monitoring of the sound pressure change of the fan. The location of Point 1 and Point 2 are shown in Figure 10, and all dimensions are in millimetres [21,22]. The noise spectrogram of the prototype fan is shown in Figure 11 when the rotational speed of the fan is 5500 r/min and the flow rate is 8.5 m<sup>3</sup> /s at atmospheric pressure. It can also be seen from Figure 11 that although Point 1 is more than twice as far away from Point 2, the sound pressure level differs little at different frequencies. At the same time, for Point 2 in the axial direction, with the continuous increase of frequency (0–50,000 Hz), the sound pressure level gradually decreases and the trend eases. In the radial direction, with the increase of frequency, the variation of sound pressure level fluctuates greatly, which is the smallest when the frequency reaches 5000 Hz. This is because the change in axial pressure is much smaller than the change in radial pressure.

**Figure 10.** The location of noise measurement Points 1 and 2.

**Figure 11.** The spectrogram of two monitoring points of centrifugal fan under normal pressure.

Changing the ambient pressure and fluid physical properties, the total sound pressure level of the prototype fan under different ambient pressures is obtained as shown in Figure 12. It can be seen from Figure 12 that as the ambient pressure decreases from 100 kPa to 50 kPa, the noise of the fan gradually becomes lower. In measurement Point 1, noise decreases from 103.85 dB to 98.2 dB. In measurement Point 2, it decreases from 121.45 dB to 118.1 dB. So, the sound pressure levels at the two measurement points decrease by about 5.8% and 2.8%, respectively. This means that the sound pressure level at negative pressure is smaller than that at one atmosphere. This also means that the noise produced by the centrifugal fan is less than that produced by one standard atmospheric pressure under the low-pressure environment such as the plateau environment. The general reason mainly lies in that the air density decreases when the air pressure is low, and the pressure variation range in the air during the operation of the centrifugal fan is small.

**Figure 12.** Fan noise under negative pressure.

At the same time, with the decrease of ambient pressure, the sound pressure level of fan noise presents an approximate logarithmic reduction trend. The relationship between sound pressure level and sound power can be transformed into

$$L\_P = 10 \lg \frac{W}{W\_0} - 20 \lg r - 11,\tag{9}$$

where *W*<sup>0</sup> is the reference sound power, and it is generally taken as 10−<sup>2</sup> W. As can be seen from Equation (9), we can obtain that *L<sup>p</sup>* ~ *lgW*. Since the sound power *W* and the air density ρ are in a linear relationship, namely *W* ~ ρ, so *L<sup>p</sup>* ~ *lg*ρ, it can be seen that the logarithmic relationship of the sound pressure level and the air density is linear.

#### **5. Conclusions**

In this paper, the flow and noise characteristics of centrifugal fan in low-pressure environment (50 kPa–100 kPa) are studied; the following conclusions can be obtained:


**Author Contributions:** Conceptualization, X.Z.; methodology, X.Z.; software, Y.Z.; validation, Y.Z.; formal analysis, C.L.; writing—original draft preparation, C.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China, grant number 51806114 and 51874187.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Processes* Editorial Office E-mail: processes@mdpi.com www.mdpi.com/journal/processes

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-1815-2