Next Article in Journal
Applying Reservoir Simulation and Artificial Intelligence Algorithms to Optimize Fracture Characterization and CO2 Enhanced Oil Recovery in Unconventional Reservoirs: A Case Study in the Wolfcamp Formation
Next Article in Special Issue
Effect of Mixing Section Acoustics on Combustion Instability in a Swirl-Stabilized Combustor
Previous Article in Journal
Enhanced Dynamic Performance in Hybrid Power System Using a Designed ALTS-PFPNN Controller
Previous Article in Special Issue
Study on Dynamic Injection Prediction Model of High-Pressure Common Rail Injector under Thermal Effect
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Comparative Study of Fuel Cavitation in Microchannels under Different Turbulence Models

College of Power Engineering, Naval University of Engineering, Wuhan 430033, China
*
Author to whom correspondence should be addressed.
Energies 2022, 15(21), 8265; https://doi.org/10.3390/en15218265
Submission received: 27 September 2022 / Revised: 23 October 2022 / Accepted: 24 October 2022 / Published: 4 November 2022
(This article belongs to the Special Issue Low-Emission Combustion Techniques: Latest Advances and Prospects)

Abstract

:
The fuel injector is a critical component of the internal combustion engine. The diameters of the injector nozzle and the control chamber’s oil inlet and outlet are generally between 0.2 and 0.5 mm, which are typical microchannel structures. During high-pressure injection, the cavitation phenomenon in the channel seriously affects the reliability of the internal combustion engine. The choice of turbulence and cavitation models is the key to investigate the cavitation in the microchannel by using numerical methods. Based on the Winklhofer microchannel fuel experiment, five representative turbulence models were used to construct a microchannel model, and the results were compared and analyzed with the experiment. The results show that the pressure gradient values obtained from the combination of RNG k-ε and ZGB models were similar to the experimental data, with an error of less than 6%. The cavitation distribution calculated from the combination of LES and ZGB models was most consistent with the experimental observation data. The outlet mass flow rate obtained from the LES and ZGB models matched the trend of the experimental data in the pressure difference range of 19 bar to 85 bar, with an error of less than 2%. For the cross-sectional flow rate calculation, the RNG k-ε and ZGB models had the smallest calculation errors, with errors below 11%.

1. Introduction

As one of the main development directions in modern internal combustion engines, high-speed and high-power density internal combustion engines are characterized by high injection pressure and speed and short injection duration. With the increasing injection pressure in a common high-pressure rail system, the cavitation problem inside the injector nozzle holes becomes more and more serious. When the pressure difference ΔP between the inlet and outlet of the injector hole is large, a low-pressure area will be generated near the wall of the flow path inside the injector hole, resulting in the pressure of high-speed fuel in this area rapidly dropping to its saturated vapor pressure value, thus leading to cavitation. Cavitation in the injector hole affects the subsequent spray combustion effect, and long-term cavitation corrosion will lead to spray hole wall peeling, affecting the injection accuracy and reliability of engine operation. Currently, the use of computational fluid dynamics (CFD) methods to validate relevant injector fuel cavitation flow tests and thus study the cavitation in injector holes is one of the hot research topics in the field of internal combustion engines.
The current fuel cavitation tests that are widely used to verify the reliability of numerical model calculations are mainly the two-dimensional injector hole test proposed by Sou et al. [1] and the microchannel test proposed by Winklhofer et al. [2,3]. Because the microchannel test is closer to the actual work of the injector hole, in recent years, scholars at home and abroad have widely used it to validate various new numerical models of the spray hole. Dai et al. [4] used the Winklhofer test for the validation of the reliability of the model of a curved structure nozzle and used the nozzle model to study the effect of nozzle curvature on the transient cavitation characteristics of the nozzle, finding that the nozzle curvature can significantly reduce the cavitation in the hole degree. The greater the arc, the smaller the area of cavitation, and the smaller the average velocity of the spray hole exit, but the mass flow rate did not change significantly. By combining the Winklhofer test data, Sa et al. [5] used the realizable k-ε turbulence model and Schnerr and Sauer (SS) cavitation model in their calculations and verified that the cavitation distribution of the constructed multiphase flow model was consistent with the test, and then quantitatively investigated the effect of the spiral counter-groove structure on the injector needle valve on the turbulent fuel flow in the orifice and the subsequent fuel injection. Cristofaro et al. [6] used the large eddy simulation (LES) model to calculate the cavitation distribution of his multiphase flow model, in agreement with the experimental data, and investigated the effect of liquid fuel viscosity on the mass flow rate, velocity profile and cavitation distribution during the microchannel throttling flow. Zhao et al. [7] used a model with a grid size of 10 μm to verify the cavitation distribution, outlet mass flow rate and cross-sectional velocity distribution of the model similarly to the experiment, and finally investigated the effect of the compressibility of the fuel on the emission coefficient and critical cavitation number of the fuel in the nozzle. Guo et al. [8] used the homogeneous relaxation model and the RNG k-ε turbulence model calculations to verify that the model used was consistent with the experimental outlet mass flow rate and that the error in the cavitation distribution was minimal, thus investigating the effect of needle valve motion on the development of cavitation in the nozzle.
Although numerical calculation results closer to the experimental phenomena can be obtained in the study of the Winklhofer test at present, the influence of the selected computational model on the computational results and experimental validation has not been reported in the above-mentioned papers because the computational models and meshing strategies used are not the same. In recent years, some scholars have conducted related studies based on the influence of cavitation model selection on simulation results. Saha et al. [9] investigated the influence of three cavitation models such as Zwart-Gerber-Belamri (ZGB), SS and Saha-Abu-Ramadan-Li (SAL) on simulation results with the selection of single-fluid and two-fluid methods. Piehl et al. [10] carried out similar research work using three cavitation models such as Original Fire and SS, and Singhal et al. and Kumar et al. [11] combined multiple URANS models with ZGB models to investigate the effect of the selection of Unsteady Reynolds Averaged Navier-Stokes (URANS) models on the computational results and verified them by comparing with the experimental data of Duke et al. [12]. He et al. [13] evaluated the predictive ability of different models for string cavitation in nozzles based on turbulence models such as RNG k-ε, SST k-ω, WALE LES and VLES, and found that LES and VLES obtained more accurate simulation results when studying the above problem. However, with the gradual use of LES and DES in microchannel cavitation simulations, it became more important to compare the differences in the effects of LES, DES, and URANS models on microchannel cavitation flow. The Winklhofer microchannel test mainly has four result indicators: pressure gradient distribution, cavitation distribution, outlet mass flow rate and cross-sectional flow rate distribution. In order to obtain more reliable research conclusions, the simulation data need to be compared and verified comprehensively with the above-proposed test parameters.
In this paper, we used Fluent 2020R2 commercial software to numerically simulate the pressure gradient, cavitation distribution, outlet mass flow rate and cross-sectional flow velocity distribution under various differential pressure conditions in the Winklhofer microchannel test based on five combinations of turbulence and cavitation models in turn, to compare and discuss the differences in numerical results among the model combinations and to evaluate and summarize the influence of model selection on the test validation results.

2. Computational Model Theory

2.1. Turbulence Model

2.1.1. SST k-ω Model

The SST k-ω model is improved from the Baseline k-ω model and inherits the robustness and accuracy of the k-ω class model in the near-wall region and the free-flow independence of the k-ε class model in the far-field region, and respectively defines the transport process of turbulent shear stress in the turbulent viscosity with the transport equations of turbulent kinetic energy and specific dissipation rate, as:
t ( ρ k ) + x i ( ρ k u i ) = x j ( Γ k k x j ) + G k Y k + G b
t ( ρ ω ) + x i ( ρ ω u i ) = x j ( Γ ω ω x j ) + G ω Y ω + D ω + G ω b
where G k is the turbulent energy generation term; G ω is the specific dissipation rate generation term; Γ k and Γ ω are the effective diffusivity of turbulent energy and specific dissipation rate, respectively; Y k and Y ω are the dissipation terms of turbulence-induced turbulent energy and specific dissipation rate, respectively; D ω is the cross-diffusion term; G b is the buoyancy-induced turbulent energy; and G ω b is the buoyancy term in the transport equation. The specific expressions of the above variables are given in ref. [14].

2.1.2. Realizable k-ε Model

Based on the mean-square vorticity fluctuation transport equation, the Realizable k-ε model derives a modified transport equation for the dissipation rate and has excellent computational performance in most flow phenomena, especially separated flows and secondary flows with complex characteristics, and its transport equations for turbulent kinetic energy and dissipation rate are:
t ( ρ k ) + x j ( ρ k u j ) = x j [ ( μ + μ t σ k ) k x j ] + G k + G b ρ ε Y M
t ( ρ ω ) + x i ( ρ ω u i ) = x j ( Γ ω ω x j ) + G ω Y ω + D ω + G ω b
where G k is the turbulent energy due to the mean velocity gradient; G b is the turbulent energy due to buoyancy; Y M is the contribution of pulsating expansion to the total dissipation rate in compressible turbulent flow; μ t is the turbulent viscosity; C 1 is a coefficient term in the equation and has:
C 1 = max [ 0.43 , η η + 5 ]
η = S k ε
S = 2 S i j S i j
where σ k and σ ε are the turbulent Prandtl numbers of turbulent kinetic energy and dissipation rate, respectively, and σ k = 1.0 and σ ε = 1.2 ; C 1 ε and C 2 are the constants in the dissipation rate transport equation, and C 1 ε = 1.44 and C 2 ε = 1.9 . The specific expressions of the above variables are given in ref. [15].

2.1.3. RNG k-ε Model

The RNG k-ε model is based on the Standard k-ε model using the statistical method of reformed group. It is characterized by the improved accuracy of the eddy calculation and the effective calculation of the low Reynolds number effect in the near-wall region. Thus the RNG k-ε model is more accurate and applicable to a wider range of flow cases than the Standard k-ε model, and the transport equations for turbulent kinetic energy and dissipation rate are:
t ( ρ k ) + x i ( ρ k u i ) = x j ( α k μ e f f k x j ) + G k + G b ρ ε Y M
t ( ρ ε ) + x i ( ρ ε u i ) = x j ( α ε μ e f f ε x j ) + C 1 ε ε k ( G k + G 3 ε G b ) C 2 ε ρ ε 2 k R ε
where G k , G b and Y M have the same physical meaning as in the Realizable k - ε model (Equation (3)); μ e f f is the effective viscosity; α k and α ε are the inverse of the effective Prandtl number of the turbulent energy and dissipation rate, respectively; C 1 ε and C 2 ε are the constants in the dissipation rate transport equation, and C 1 ε = 1.42 and C 2 ε = 1.68 ; R ε are additional terms in the dissipation rate transport equation. The specific expressions of the above variables are given in ref. [16].

2.1.4. LES Model

In the LES model, large eddies are solved by direct analysis, while small eddies are solved by constructing a subgrid model. This paper uses the wall-adapting local eddy-viscosity (WALE) model for the large eddy simulation [17]. Compared with the Smagorinsky–Lilly model [18], the WALE model treats the turbulent viscosity of the laminar shear flow as zero, which allows the WALE model to correctly handle the flow in the laminar region of the basin, improving the accuracy of the solution at the near-wall turbulence. The eddy viscosity of the WALE model is
μ t = ρ L S 2 ( S i j d S i j d ) 3 / 2 ( S ¯ i j S ¯ i j ) 5 / 2 + ( S i j d S i j d ) 5 / 4
where S ¯ i j is the solvable scale strain rate tensor; L S and S i j d respectively are
L S = min ( κ d , C w V 1 / 3 )
S i j d = 1 2 ( g ¯ i j 2 + g ¯ j i 2 ) 1 3 δ i j g ¯ k k 2
where, κ is the von Kamen constant, and κ = 0.41 ; d is the distance to the nearest wall; C w is the WALE model constant, and C w = 0.325 ; g ¯ i j is the velocity gradient tensor, and g ¯ i j = u i ¯ x j .

2.1.5. DES Model

The DES model is classified as a hybrid RANS/LES model. It obtains higher accuracy than the URANS model by using the URANS model in the boundary layer, while switching to the LES model in the outer field region where the separated vortices are generated, and avoids the situation where a large amount of computational resources are occupied in the LES model. In this paper, we use the SST k - ω based DES model [19], whose dissipation terms of turbulent kinetic energy are
Y k = ρ β * k ω F D E S
where F D E S is
F D E S = max ( L t C d e s Δ max , 1 )
where L t is the turbulence length scale, and L t = k β * ω ; C D E S is the calibration constant in the DES model, and C D E S = 0.61 ; Δ max is the grid spacing.

2.2. Cavitation Model

2.2.1. Cavitation Model Basic Control Equations

The cavitation model can constitute a multiphase flow cavitation modeling method together with the multiphase flow model for controlling gas–liquid mixtures and the traditional turbulence model, whose control equations are expressed as:
t ( f v ρ ) + ( f v ρ V v ) = ( Γ f v ) + R e R c
where f v is the vapor mass fraction; ρ is the mixture density; V v is the vapor velocity; Γ is the diffusion coefficient; R e is the vapor generation rate; and R c is the vapor condensation rate.

2.2.2. ZGB Model

The ZGB model assumes that the bubbles in the multiphase flow system all have the same size, neglects the effect of the non-condensable gas on the cavitation flow, and considers the interphase mass transfer rate per unit volume using the rate of change of the mass of a single bubble and the bubble number density [20].
When P P v :
R e = F v a p 3 α n u c ( 1 α v ) ρ v B 2 3 P v P ρ l
When P P v :
R c = F c o n d 3 α v ρ v B 2 3 P P v ρ l
where P is the local pressure of the flow field; P v is the saturation vapor pressure; B is the bubble radius, B = 1 × 10 6 m ; α n u c is the volume fraction at the gas nucleus, α n u c = 5 × 10 4 ; F v a p is the evaporation coefficient, F v a p = 50 ; F c o n d is the condensation coefficient, F c o n d = 0.01 ; ρ v and ρ l are the gas phase density and liquid phase density, respectively.

3. Numerical Model Construction and Validation

3.1. Physical Model Construction

In this paper, a “U” shaped pipe from the Winklhofer visualization fuel experiment [3] was selected as the object of our study, and its geometry is shown in Figure 1. The geometrical parameters of the “U” pipe model are: pipe thickness W = 300   μ m , microchannel inlet height H i n = 301   μ m , microchannel outlet height H o u t = 284   μ m , microchannel length L = 1000   μ m , and microchannel inlet corner radius R = 20   μ m . The previous studies [6,7,21] all used microchannel models with additional preset buffer flow channels to verify with the Winklhofer test results for comparison. The purpose was to keep the inlet and outlet pressures of the microchannel in the simulation consistent with the pressure conditions in the test so that the simulation results are closer to the actual test measurements. As shown in Figure 2, the simulation model of Guan [21] is referred to in this paper, and the buffered flow channels of 2000   μ m × 300   μ m × 2312   μ m and 3000   μ m × 300   μ m × 2312   μ m are preset before and after the microchannel model, respectively.

3.2. Numerical Calculation Settings

The fuel thermal property parameters used in the numerical calculations are shown in Table 1 [21], and the fuel is set as a compressible fluid. In addition, the microchannel inlet and outlet are defined as pressure inlet conditions and pressure outlet conditions, respectively, keeping the inlet pressure of 100 bar constant and setting different outlet pressure conditions (19 to 85 bar) for the computational model.
The discretization of the numerical model is based on the finite volume method (FVM), the Mixture model is used to describe the gas–liquid two-phase flow, the PISO algorithm is used for the pressure–velocity coupling method, and the pressure solution format and the gradient solution format are set to PRESTO! The rest of the formats are set to QUICK format. The calculation step is set to 5 × 10−7 s. The specific parameters of each turbulence and cavitation model are shown in Table 2.

3.3. Grid-Independent Verification

For the microchannel basin and the buffer basin near the entrance and exit of the microchannel, local grid encryption was performed using ICEM software to obtain accurate flow parameter results for subsequent model calculations, and the grid division is shown in Figure 3.
For the above-mentioned grid division strategy, the grid irrelevance was verified by using the grids of magnitude 50,000, 200,000, 400,000, 600,000 and 800,000 in turn. The Winklhofer flow experiment was able to observe the supercavitation phenomenon in the microchannel under the condition of 80 bar differential inlet and outlet pressure. Therefore, the grid-independent verification is based on the above pressure conditions, and the mass flow rate at the exit of the microchannel for five grid numbers is shown in Figure 4. The results show that the mass flow rate at the exit of the microchannel tended to a stable value when the number of grids exceeded 600,000. Considering the accuracy of subsequent numerical calculations and the allocation of computational resources, we finally determined that the grid number of the computational model was 600,000, in which the minimum grid size at the near-wall surface was 0.5 μm and the maximum grid size at the central flow channel was about 4 μm.

4. Calculation Results and Discussion

When fuel flows through the microchannel, the flow channel cross-section suddenly becomes smaller, the flow velocity increases and the pressure drops below the saturated vapor pressure of fuel, at which time flow separation and cavitation will occur. These phenomena directly affect the pressure and velocity distribution in the microchannel. They are reflected in the resultant parameters such as outlet mass flow rate, cavitation area distribution, flow velocity distribution and pressure gradient distribution. In this paper, the distribution characteristics of the above physical parameters are studied separately, the differences in results under multiple numerical model application scenarios are compared and analyzed, and finally, the accuracy of the numerical calculation of compressible fuel under multiple model applications is discussed.

4.1. Comparison of Pressure Gradient

In order to elucidate the intrinsic mechanism of cavitation incipient to cavitation development in the channel, this paper first analyzes the pressure gradient variation along the flow direction at the centerline of the model under the ΔP = 70 bar, as shown in Figure 5. Under this pressure condition, the fuel flow in the microchannel reaches the critical cavitation condition. However, the fuel pressure value first decreases rapidly from the channel inlet to the middle of the cavitation region because the cavitation region is not fully developed at this time. It drops to a small value in the middle of the cavitation region with the smallest local basin cross-section, and then the basin cross-section increases and the pressure value rises back. After the fuel fluid flows through the end of the cavitation region and reattaches back to the wall, the basin cross-section gradually decreases and the flow rate gradually increases as the fluid continues to approach the channel exit, resulting in a gradual decrease in pressure as well.
Figure 5 illustrates the differences in pressure gradient profile changes for multiple turbulence and cavitation model applications. The pressure gradient curves obtained from different model combinations almost overlap and are similar to the experimental results as the fuel flows from the upstream buffer basin to the vicinity of the microchannel inlet. As the fuel flows through the cavitation and non-cavitation regions, all five model combinations, except for the RNG k-ε and ZGB model combinations, show a trend where the pressure value first drops to a minimal value, then rises, and finally decreases. The errors between the pressure gradient results of the five models and the experimental data in each zone are given in Table 3. As can be seen from Table 3, except for the model combinations of LES and ZGB, the results obtained in the two zones of −1.00~−0.25 mm and −0.25~0.22 mm are close to the experimental data, and the errors are within 6%. In the 0.22–0.50 mm range, the best simulation results were obtained for the RNG k-ε and ZGB model combination with an error within 6%; in the 0.50–2.00 mm range, the simulation results for the RNG k-ε and ZGB, DES and ZGB model combinations were within 4% error. In general, the error between the pressure gradient simulation results obtained from the combination of RNG k-ε and ZGB models and the experimental data is the smallest. From Equation (9), it can be seen that the RNG k-ε model is based on the Standard k-ε model with the addition of terms that can improve the accuracy for fast strain flow [16]. Therefore, the RNG k-ε model shows superior simulation performance in the two main zones of fast strain flow in microchannels, −0.25 to 0.22 mm and 0.22 to 0.50 mm.

4.2. Comparison of Cavitation Distribution

Figure 6 compares the cavitation distribution results between the experimental results and the simulation of various model combinations. Due to the ΔP between the front and rear buffer regions and the increase in fuel flow rate in the microchannel, the fuel pressure in the local region decreases to the fuel evaporation pressure, and the fuel in the region instantly evaporates into the gas phase, and the cavitation occurs.
From the vapor volume fraction cloud diagram, it can be seen that there are significant differences between the cloud diagrams obtained from different model combinations and the experimental results. Due to the presence of non-condensable gases in the actual fuel flow, the fuel liquid is more prone to gas nucleation during the pressure drop, and cavitation is more likely to occur. At 60 bar ΔP, the cavitation at the channel’s entrance was simulated by all five model combinations, consistent with the experimental results. For 70 bar ΔP, the cavitation clouds obtained by the LES and ZGB model combinations were close to the experimental observations. In contrast, the cavitation domain lengths obtained by the other five model combinations were smaller than those in the experimental results. For the pressure difference of 80 bar, the simulation results of the LES and ZGB and DES and ZGB model combinations were able to better reflect the experimental data, while the cavitation domain lengths obtained from the SST k - ω and ZGB model combinations were slightly smaller than the experimental values.
It can be seen from the experimental results that as the ΔP increased from 60 to 80 bar, the gas phase fuel was no longer confined to the reflux region at the entrance of the microchannel, but gradually extended along the wall to the exit. However, the cavitation domains observed in the microchannels in the experiments were still larger than those obtained in the simulations, and the cavitation distributions obtained from different turbulence models were different. Yu et al. [22] showed that different cavitation distributions were observed using different cavitation models for the same microchannel structure. Altimira et al. [23] also found that the cavitation domain was larger in the experiment than in the simulation and pointed out that fuel compressibility affects the vapor condensation rate, which may be the reason for the small cavitation domain in the simulation. Meanwhile, he further found that using the different densities of grid distribution and adjusting the size of control parameters affecting cavitation bubble development in the cavitation model did not significantly affect the cavitation distribution at the near-wall surface. Giannadakis et al. [24] showed that the cavitation domain development did not depend on the initial diameter of the cavitation bubble.
In summary, the cavitation domain attached to the microchannel wall developed as the ΔP between inlet and outlet increased, gradually extending along the wall to the outlet of the channel. The reason for the variability in the cavitation distribution between different model combinations at the same inlet and outlet pressure difference was the difference in the mathematical formulation of the different turbulence models because the cavitation domain was mainly composed of multiple large vortices clustered at the macroscopic level (one can see the wake of the vortices attached from the wall in the cavitation cloud at 60 bar pressure difference from the Winklhofer experiment). The above error analysis process also reflects that the LES model directly resolves the large vortices in a more suitable way for the simulation calculation of the cavitation domain than the way the remaining four turbulence models deal with the flow.

4.3. Comparison of Outlet Mass Flow Rates

Figure 7 reveals the microchannel outlet mass flow rate’s variation characteristics and the differences in the simulation results obtained from five model combinations. In this paper, the outlet mass flow rate is studied for thirteen ΔP conditions (ΔP1 = 19 bar, ΔP2 = 45 bar, ΔP3 = 58 bar, ΔP4 = 60 bar, ΔP5 = 63 bar, ΔP6 = 65 bar, ΔP7 = 67 bar, ΔP8 = 69 bar, ΔP9 = 70 bar, ΔP10 = 71 bar, ΔP11 = 75 bar, ΔP12 = 80 bar, ΔP13 = 85 bar), which are entirely consistent with the conditions used in the experiments.
From the experimental results, it can be seen that in the process of gradually increasing the ΔP from 19 bar to 85 bar, the outlet mass flow rate firstly increases with the increase in ΔP between the inlet and outlet. Then, the mass flow rate curve shows an inflection point at 70 bar ΔP, after which the outlet mass flow rate maintains a certain value and no longer changes with the increase in ΔP. As shown in Figure 7, numerical results obtained by all five model combinations were in good agreement with the experimental data in the range of 19 bar to 70 bar ΔP. As can also be seen from Table 4, by dividing the model results into two zones (19 to 70 bar and 70 to 85 bar) with different variation trends of outlet mass flow rate according to the experimental data, the error between the results obtained from the five models and the experiments were within 6%. Among them, the simulation error of the LES and ZGB model combination was the smallest, and the error was within 2% in both zones. After a pressure difference of 70 bar, the results of the four model combinations, except for the SST k - ω and ZGB model combinations, were slightly larger than the experimental data, but the trend of no further increase in the outlet mass flow rate could be reproduced. The LES (WALE) model assumes zero turbulent viscosity for the laminar shear flow, an assumption that allows the model to treat the laminar region in the domain correctly [18]. As a result, the LES(WALE) model yielded simulation results that were in general agreement with the experimental data during the process in which the microchannel undergoes cavitation priming at 60 bar ΔP and subsequently progresses to sheet cavitation at 70 bar ΔP.

4.4. Comparison of Velocity Distribution

The velocity distribution along the y-axis for two cross-sections (x1 = 0.053 mm; x2 = 0.17 mm) within the microchannel at 55 bar and 67 bar ΔP are shown in Figure 8. Since the gasified fuel occupies the cavitation region when the ΔP is 67 bar, it is not easy to measure the velocity variation at the near-wall surface during the experiment [5]. According to the trend of cross-sectional flow velocity distribution in the Winklhofer test, the simulation results of cross-sections x1 and x2 under two pressure differences can be divided into sections to calculate the average value of flow velocity and tabulated and summarized, as shown in Table 5, Table 6, Table 7 and Table 8.
Under the 55 bar ΔP condition, the simulation results at x1 can somewhat reflect the velocity variation trend in the experimental data, as shown in Figure 8b,d. On the other hand, at x2, all model combinations cannot reproduce the complex flow velocity variation in the range of about 50~100 μm away from the upper and lower walls of the microchannel. As can be seen from Table 5 and Table 7, the simulation results are still quite different from the experimental data in the near-wall basin at a distance of about 0~50 μm from the upper and lower walls. When the y-axis position point is about 50~250 μm (i.e., the central basin of cross-sections x1 and x2), the results obtained from the RNG k - ε and ZGB models are closest to the experimental data, i.e., the calculation errors in both cross-sections are within 11%. The calculation errors of the DES and ZGB models are the largest, and the calculation errors in both cross-sections are over 14%.
Combined with Figure 8c,e, the results obtained from the combination of the five models all match the trend of the test data at section x1 under the 67 bar ΔP condition, but this still makes it difficult to reproduce the flow velocity distribution in the range of about 50–100 μm away from the wall at x2. As can be seen from Table 6 and Table 8, the computational error of the flow velocity distribution near the wall at the 67 bar ΔP condition is slightly improved compared to that at the 55 bar ΔP condition. However, the simulation results are still somewhat different from the experimental data. For the central basin of cross-sections x1 and x2, the calculation errors of RNG k - ε and ZGB models are still the smallest, i.e., the calculation errors at both cross-sections are within 8%. In contrast, DES and ZGB have the largest computational errors, which exceed 19% in both cross-sections.
Since the RNG k - ε model can construct an effective viscosity term based on the RNG theory that takes into account the low Reynolds number effect [16], this means that the RNG k - ε model has superior performance in dealing with the near-wall flow or microchannel flow where the wall surface has a strong influence on the flow. Furthermore, the RNG k - ε model retains the computational accuracy of the k - ε type model for the high Reynolds number flow region. Therefore, we believe that the model treatment described above is the main reason why the RNG k - ε and ZGB models eventually yield the simulation results closest to the experimental data.

5. Conclusions

In this paper, the cavitation phenomena in the experimental model of Winklhofer microchannel under different ΔP conditions were investigated numerically using the SST k-ω, Realizable k-ε, RNG k - ε , LES and DES turbulence models combined with the ZGB cavitation model, respectively. The following conclusions were obtained:
(1)
At 70 bar ΔP, although the combination of RNG k - ε and ZGB models do not reflect the trend where the pressure along the flow direction at the centerline of the microchannel drops to a very small value, then rises, and finally decreases, the RNG k - ε and ZGB model combination has the best agreement with the experimental value compared with the other four model combinations, and its error is controlled within 6% in each zone.
(2)
By comparing the experimental data with the simulation results, the cavitation domain at the near wall develops continuously during the increase in the inlet and outlet pressure difference and gradually extends along the wall to the microchannel outlet. Under the same differential pressure condition, the cavitation domain length obtained from simulation is smaller than that obtained from experimental observation. The cavitation distribution varies among different model combinations. The LES and ZGB model combinations can obtain simulation results similar to the cavitation distribution of the Winklhofer test at 60 bar, 70 bar and 80 bar ΔP conditions. In comparison, the simulation results of the remaining five model combinations at 70 bar and 80 bar ΔP conditions are still somewhat different from the experimental data.
(3)
For the numerical calculation of the outlet mass flow rate, all five model combinations can obtain numerical results in the range of pressure difference from 19 bar to 70 bar, in good agreement with the experimental data. After the differential pressure exceeds 70 bar, the calculation results of the remaining four model combinations can reflect the trend of the experimental data, except for the SST k - ω and ZGB model combinations. Among them, the LES and ZGB model combinations have the best simulation results, and the errors in both study sections are within 2%.
(4)
For the numerical calculation of the flow velocity distribution along the y-axis, the results obtained from the simulation of the five model combinations still differed from the experimental data in the near-wall basin at a distance of about 0~50 μm from the upper and lower walls. In the central basin of cross-sections x1 and x2, the trends of the calculated results of RNG k - ε and ZGB models are in the best agreement with the experimental data, and the calculated errors are below 11% and 8% for the two differential pressure conditions of 55 bar and 67 bar, respectively.
(5)
Compared with the Winklhofer test data, the RNG k - ε and ZGB models are suitable for calculating the pressure gradient variation at the centerline of the microchannel and the flow velocity distribution in a certain cross-section of the microchannel when cavitation occurs. The LES and ZGB models are suitable for calculating the injector microchannel’s cavitation domain development and outlet mass flow rate.

Author Contributions

Conceptualization, Z.L. (Ziming Li), Z.L. (Zhenming Liu) and P.C.; methodology, Z.L. (Ziming Li), Z.L. (Zhenming Liu) and P.C.; software, Z.L. (Ziming Li); validation, Z.L. (Ziming Li) and Z.L. (Zhenming Liu); formal analysis, Z.L. (Ziming Li); resources, J.L. and J.W.; writing—original draft preparation, Z.L. (Ziming Li) and Z.L. (Zhenming Liu); writing—review and editing, Z.L. (Ziming Li), Z.L. (Zhenming Liu), P.C., J.L. and J.W.; supervision, Z.L. (Zhenming Liu); funding acquisition, Z.L. (Zhenming Liu). 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 No. 51879269.

Data Availability Statement

Some data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

C 1 A coefficient term of Realizable k-ε model
C 1 ε A constant in the dissipation rate transport equation
C 2 ε A constant in the dissipation rate transport equation
C w WALE model constant
C D E S Calibration constant in the DES model
d Distance to the nearest wall (m)
D ω Cross-diffusion term
f v Vapor mass fraction
F c o n d Condensation coefficient
F v a p Evaporation coefficient
g ¯ i j Velocity gradient tensor term
G b Buoyancy-induced turbulent energy term
G k Turbulent energy generation term
G ω Specific dissipation rate generation term
G ω b Buoyancy term in the transport equation
H i n Microchannel inlet height (μm)
H o u t Microchannel outlet height (μm)
L Microchannel length (μm)
L t Turbulence length scale (m)
P Pressure (Pa)
P v Saturation vapor pressure (Pa)
ΔPDifferential pressure (Pa)
R Microchannel inlet corner radius (μm)
R c Vapor condensation rate term
R e Vapor generation rate term
R ε Additional terms in the dissipation rate transport equation
S ¯ i j Solvable scale strain rate tensor term
V v Vapor velocity (m/s)
Y k Dissipation terms of turbulence-induced turbulent energy rate term
Y ω Dissipation terms of turbulence-induced specific dissipation rate term
Y M Contribution of pulsating expansion to the total dissipation rate term
α n u c Volume fraction at the gas nucleus
σ k Turbulent Prandtl numbers of turbulent kinetic energy term
σ ε Turbulent Prandtl numbers of dissipation rate term
κ Von Kamen constant
μ e f f Effective viscosity (Pa·s)
Γ Diffusion coefficient
Γ k Effective diffusivity of turbulent energy term
Γ ω Specific dissipation rate term

Abbreviation

CFDComputational fluid dynamics
DESDetached Eddy Simulation
FVMFinite volume method
LESLarge Eddy Simulation
RNGRenormalization-Group
SALSaha-Abu-Ramadan-Li
SS Schnerr and Sauer
SST Shear Stress Transport
URANS Unsteady Reynolds Averaged Navier-Stokes
VLESVery Large Eddy Simulation
ZGBZwart-Gerber-Belamri

References

  1. Sou, A.; Tomiyama, A.; Hosokawa, S.; Nigorikawa, S.; Maeda, T. Cavitation in a two-dimensional nozzle and liquid jet atomization (LDV measurement of liquid velocity in a nozzle). JSME Int. J. Ser. B Fluids Therm. Eng. 2006, 49, 1253–1259. [Google Scholar] [CrossRef] [Green Version]
  2. Winklhofer, E.; Philipp, H.; Hirsch, A.; Morozov, A. Cavitation and spray formation in diesel flow situations. In Proceedings of the ILASS-Europe 2000, Darmstadt, Germany, 11–13 September 2000; pp. 11–13. [Google Scholar] [CrossRef]
  3. Winklhofer, E.; Kull, E.; Kelz, E.; Morozov, A. Comprehensive hydraulic and flow field documentation in model throttle experiments under cavitation conditions. In Proceedings of the ILASS-Europe 2001, Zurich, Switzerland, 2–6 September 2001; pp. 574–579. [Google Scholar] [CrossRef]
  4. Dai, Y.; Zhang, X.; Zhang, G.; Cai, M.; Zhou, C.; Ni, Z. Numerical analysis of influence of cavitation characteristics in nozzle holes of curved diesel engines. Flow Meas. Instrum. 2022, 85, 102172. [Google Scholar] [CrossRef]
  5. Sa, B.; Klyus, O.; Markov, V.; Kamaltdinov, V. A numerical study of the effect of spiral counter grooves on a needle on flow turbulence in a diesel injector. Fuel 2021, 290, 120013. [Google Scholar] [CrossRef]
  6. Cristofaro, M.; Edelbauer, W.; Koukouvinis, P.; Gavaises, M. Influence of Diesel Fuel Viscosity on Cavitating Throttle Flow Simulations under Erosive Operation Conditions. ACS Omega 2020, 5, 7182–7192. [Google Scholar] [CrossRef] [PubMed]
  7. Zhao, J.; Liu, W.; Zhao, J.; Grekhov, L. Numerical investigation of gas/liquid two-phase flow in nozzle holes considering the fuel compressibility. Int. J. Heat Mass Transf. 2020, 147, 118991. [Google Scholar] [CrossRef]
  8. Guo, G.; He, Z.; Wang, Q.; Lai, M.-C.; Zhong, W.; Guan, W.; Wang, J. Numerical investigation of transient hole-to-hole variation in cavitation regimes inside a multi-hole diesel nozzle. Fuel 2021, 287, 119457. [Google Scholar] [CrossRef]
  9. Saha, K.; Battistoni, M.; Som, S.; Li, X. Modeling of Cavitation in Fuel Injectors with Single- and Two-Fluid Approaches. In Two-Phase Flow for Automotive and Power Generation Sectors; Springer: Singapore, 2019; pp. 185–201. [Google Scholar] [CrossRef]
  10. Piehl, J.; Bravo, L. Assessment of Cavitation Models for Computational Fluid Dynamics Analysis of Erosion Risk in a Hydrocarbon-Fueled Nozzle; US Army Research Laboratory: Adelphi, MD, USA, 2018. [Google Scholar]
  11. Kumar, A.; Ghobadian, A.; Nouri, J.M. Assessment of Cavitation Models for Compressible Flows Inside a Nozzle. Fluids 2020, 5, 134. [Google Scholar] [CrossRef]
  12. Duke, D.J.; Kastengren, A.L.; Tilocco, F.Z.; Powell, C.F. Synchrotron X-ray measurements of cavitation. In Proceedings of the ILASS-Americas 25th Annual Conference on Liquid Atomization and Spray Systems, Pittsburgh, PA, USA, 5–8 May 2013; pp. 21–23. [Google Scholar]
  13. He, Z.; Guan, W.; Wang, C.; Guo, G.; Zhang, L.; Gavaises, M. Assessment of turbulence and cavitation models in prediction of vortex induced cavitating flow in fuel injector nozzles. Int. J. Multiph. Flow 2022, 157, 104251. [Google Scholar] [CrossRef]
  14. Menter, F.R. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 1994, 32, 1598–1605. [Google Scholar] [CrossRef] [Green Version]
  15. Shih, T.-H.; Liou, W.W.; Shabbir, A.; Yang, Z.; Zhu, J. A new k-ϵ eddy viscosity model for high reynolds number turbulent flows. Comput. Fluids 1995, 24, 227–238. [Google Scholar] [CrossRef]
  16. Yakhot, V.; Orszag, S.A. Renormalization group analysis of turbulence. I. Basic theory. J. Sci. Comput. 1986, 1, 3–51. [Google Scholar] [CrossRef]
  17. Nicoud, F.; Ducros, F. Subgrid-Scale Stress Modelling Based on the Square of the Velocity Gradient Tensor. Flow Turbul. Combust. 1999, 62, 183–200. [Google Scholar] [CrossRef]
  18. Smagorinsky, J. General circulation experiments with the primitive equations: I. The basic experiment. Mon. Weather. Rev. 1963, 91, 99–164. [Google Scholar] [CrossRef]
  19. Menter, F.R.; Kuntz, M.; Langtry, R. Ten years of industrial experience with the SST turbulence model. Turbul. Heat Mass Transf. 2003, 4, 625–632. [Google Scholar]
  20. Zwart, P.J.; Gerber, A.G.; Belamri, T. A two-phase flow model for predicting cavitation dynamics. In Proceedings of the ICMF 2004, 5th International Conference on Multiphase Flow, Yokohama, Japan, 30 May–3 June 2004. [Google Scholar]
  21. Guan, W.; He, Z.; Zhang, L.; El-Seesy, A.I.; Wen, L.; Zhang, Q.; Yang, L. Effect of asymmetric structural characteristics of multi-hole marine diesel injectors on internal cavitation patterns and flow characteristics: A numerical study. Fuel 2021, 283, 119324. [Google Scholar] [CrossRef]
  22. Yu, H.; Goldsworthy, L.; Brandner, P.; Garaniya, V. Development of a compressible multiphase cavitation approach for diesel spray modelling. Appl. Math. Model. 2017, 45, 705–727. [Google Scholar] [CrossRef]
  23. Altimira, M.; Fuchs, L. Numerical investigation of throttle flow under cavitating conditions. Int. J. Multiph. Flow 2015, 75, 124–136. [Google Scholar] [CrossRef]
  24. Giannadakis, E.; Gavaises, M.; Arcoumanis, C. Modelling of cavitation in diesel injector nozzles. J. Fluid Mech. 2008, 616, 153–193. [Google Scholar] [CrossRef]
Figure 1. The microchannel geometry of the Winklhofer.
Figure 1. The microchannel geometry of the Winklhofer.
Energies 15 08265 g001
Figure 2. The simulation calculation model of the microchannel structure.
Figure 2. The simulation calculation model of the microchannel structure.
Energies 15 08265 g002
Figure 3. The simulation mesh model of the microchannel structure.
Figure 3. The simulation mesh model of the microchannel structure.
Energies 15 08265 g003
Figure 4. Verification of grid independence based on 80 bar ΔP conditions.
Figure 4. Verification of grid independence based on 80 bar ΔP conditions.
Energies 15 08265 g004
Figure 5. Comparison of the pressure gradient results at the centerline of the model under different turbulence models.
Figure 5. Comparison of the pressure gradient results at the centerline of the model under different turbulence models.
Energies 15 08265 g005
Figure 6. Comparison of cavitation distribution regions at microchannels under different turbulence models.
Figure 6. Comparison of cavitation distribution regions at microchannels under different turbulence models.
Energies 15 08265 g006
Figure 7. Comparison of mass flow rate results at the exit of microchannels under different turbulence models.
Figure 7. Comparison of mass flow rate results at the exit of microchannels under different turbulence models.
Energies 15 08265 g007
Figure 8. Comparison of flow velocity distribution in microchannel cross-section under different turbulence models. (a) The schematic diagram of the location of the two cross-sections within the microchannel; (b) Velocity distribution of section x1 in 55 bar ΔP; (c) Velocity distribution of section x1 in 67 bar ΔP; (d) Velocity distribution of section x2 in 55 bar ΔP; (e) Velocity distribution of section x2 in 67 bar ΔP.
Figure 8. Comparison of flow velocity distribution in microchannel cross-section under different turbulence models. (a) The schematic diagram of the location of the two cross-sections within the microchannel; (b) Velocity distribution of section x1 in 55 bar ΔP; (c) Velocity distribution of section x1 in 67 bar ΔP; (d) Velocity distribution of section x2 in 55 bar ΔP; (e) Velocity distribution of section x2 in 67 bar ΔP.
Energies 15 08265 g008aEnergies 15 08265 g008b
Table 1. ISO4113 fuel thermophysical properties.
Table 1. ISO4113 fuel thermophysical properties.
ParametersValue
Liquid phase density (kg/m3)830
Liquid phase kinetic viscosity (kg/(m·s))0.0024
Saturated vapor pressure (Pa)2000
Vapor phase density (kg/m3)0.029
Gas phase kinetic viscosity (kg/(m·s))3.1 × 10−6
Table 2. The parameter settings for each model during numerical calculations.
Table 2. The parameter settings for each model during numerical calculations.
SST   k - ω Realizable   k - ε RNG   k - ε LESDESZGB
σ k , 1 1.176 σ k 1.0 C 1 ε 1.42 C w 0.325 C d e s 0.61 F v a p 50
σ ω , 1 2.0
σ k , 2 1.0 σ ε 1.2 C 2 ε 1.68 F c o n d 0.01
σ ω , 2 1.168
Table 3. Comparison of the error of the simulated average value of pressure gradient for each model combination with the experimental data.
Table 3. Comparison of the error of the simulated average value of pressure gradient for each model combination with the experimental data.
Model CombinationsZone/mm
−1.00~−0.25−0.25~0.220.22~0.500.50~2.00
SST k-omega&ZGB2.55%2.22%15.14%18.61%
Realizable k-epsilon&ZGB2.13%0.79%11.44%14.51%
RNG k-epsilon&ZGB2.38%4.18%5.68%3.48%
LES&ZGB1.20%16.44%64.68%22.27%
DES&ZGB2.16%5.26%27.02%3.03%
Table 4. Comparison of the error of the simulated average value of outlet mass flow rate for each model combination with the experimental data.
Table 4. Comparison of the error of the simulated average value of outlet mass flow rate for each model combination with the experimental data.
Model CombinationsZone/bar
19~7070~85
SST k-omega&ZGB3.14%5.20%
Realizable k-epsilon&ZGB0.60%3.13%
RNG k-epsilon&ZGB0.44%3.15%
LES&ZGB0.53%1.14%
DES&ZGB3.67%5.85%
Table 5. The error of average cross-sectional x1 flow rate at 55 bar differential pressure compared to test data.
Table 5. The error of average cross-sectional x1 flow rate at 55 bar differential pressure compared to test data.
Model CombinationsZone/mm
0~4848~150150~266266~300
SST k-omega&ZGB38.27%12.02%15.27%30.61%
Realizable k-epsilon&ZGB43.14%8.02%11.95%32.34%
RNG k-epsilon&ZGB41.11%1.93%5.87%31.63%
LES&ZGB42.25%14.07%15.44%32.17%
DES&ZGB34.04%14.81%18.78%40.44%
Table 6. The error of average cross-sectional x1 flow rate at 67 bar differential pressure compared to test data.
Table 6. The error of average cross-sectional x1 flow rate at 67 bar differential pressure compared to test data.
Model CombinationsZone/mm
0~4646~150150~260260~300
SST k-omega&ZGB18.88%13.98%9.33%24.82%
Realizable k-epsilon&ZGB20.57%9.15%5.01%27.61%
RNG k-epsilon&ZGB17.59%7.10%2.87%23.75%
LES&ZGB13.57%20.69%15.76%27.93%
DES&ZGB17.86%14.62%10.36%29.28%
Table 7. The error of average cross-sectional x2 flow rate at 55 bar differential pressure compared to test data.
Table 7. The error of average cross-sectional x2 flow rate at 55 bar differential pressure compared to test data.
Model CombinationsZone/mm
0~5555~150150~254254~300
SST k-omega&ZGB69.51%23.58%20.86%18.89%
Realizable k-epsilon&ZGB54.83%18.93%14.89%12.15%
RNG k-epsilon&ZGB34.17%10.12%6.42%4.11%
LES&ZGB78.30%28.96%20.97%43.61%
DES&ZGB39.17%28.37%23.28%2.70%
Table 8. The error of average cross-sectional x2 flow rate at 67 bar differential pressure compared to test data.
Table 8. The error of average cross-sectional x2 flow rate at 67 bar differential pressure compared to test data.
Model CombinationsZone/mm
0~6363~150150~256256~300
SST k-omega&ZGB17.11%19.05%14.55%31.18%
Realizable k-epsilon&ZGB22.72%10.99%6.94%35.13%
RNG k-epsilon&ZGB22.71%7.79%4.92%32.78%
LES&ZGB12.05%29.22%19.18%47.27%
DES&ZGB18.31%19.65%14.39%35.17%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Li, Z.; Liu, Z.; Chen, P.; Liu, J.; Wu, J. Numerical Comparative Study of Fuel Cavitation in Microchannels under Different Turbulence Models. Energies 2022, 15, 8265. https://doi.org/10.3390/en15218265

AMA Style

Li Z, Liu Z, Chen P, Liu J, Wu J. Numerical Comparative Study of Fuel Cavitation in Microchannels under Different Turbulence Models. Energies. 2022; 15(21):8265. https://doi.org/10.3390/en15218265

Chicago/Turabian Style

Li, Ziming, Zhenming Liu, Ping Chen, Jingbin Liu, and Jiechang Wu. 2022. "Numerical Comparative Study of Fuel Cavitation in Microchannels under Different Turbulence Models" Energies 15, no. 21: 8265. https://doi.org/10.3390/en15218265

APA Style

Li, Z., Liu, Z., Chen, P., Liu, J., & Wu, J. (2022). Numerical Comparative Study of Fuel Cavitation in Microchannels under Different Turbulence Models. Energies, 15(21), 8265. https://doi.org/10.3390/en15218265

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop