Mathematical model
The flow area of the electronic expansion valve experiences abrupt contraction and expansion, leading to a pronounced throttling effect. As a result, the valve exhibits notable characteristics such as phase transition, vortex generation, and cavitation evolution. To accurately describe the phase transition process caused by throttling, the Reynolds‒averaged Navier‒Stokes (RANS) equation method was employed. In this study, a simplified three‒dimensional (3D) model of the flow region of the electronic expansion valve was developed, and Fluent software was utilized to incorporate various models, including turbulence model, multiphase model, cavitation model, and noise model for cavitation two – phase flow based on the RANS method. These models were allowed to capture the intricate details of the cavitation flow.
Turbulence model
The two-equation model is well-suited for flows driven by differential pressure. In view of the intricate internal geometry of the electronic expansion valve, post-valve throttling leads to an elevation in flow velocity and a transition to a high-Reynolds-number turbulent state. Thus, the k – ε turbulence model proves to be a suitable choice. Considering the presence of a pressure gradient throughout the flow process, the impacts of viscous interactions among fluid molecules on the flow dynamics cannot be neglected. Therefore, the Realizable k – ε turbulence model is selected.
$$\frac{\partial }{{\partial t}}\left( {\rho k} \right)+\frac{\partial }{{\partial {x_j}}}\left( {\rho k{u_j}} \right)=\frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu +\frac{{{\mu _t}}}{{{\sigma _k}}}} \right)\frac{{\partial k}}{{\partial {x_j}}}} \right]+{G_k}+{G_b} – \rho \varepsilon – {Y_M}+{S_k}$$
(1)
$$\frac{\partial }{{\partial t}}\left( {\rho \varepsilon } \right)+\frac{\partial }{{\partial {x_j}}}\left( {\rho \varepsilon {u_j}} \right)=\frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu +\frac{{{\mu _t}}}{{{\sigma _\varepsilon }}}} \right)\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \right]+\rho {C_1}S\varepsilon – \rho {C_2}\frac{{{\varepsilon ^2}}}{{k+\sqrt {\nu \varepsilon } }}+{C_{1\varepsilon }}\frac{\varepsilon }{k}{C_{3\varepsilon }}{G_b}+{S_\varepsilon }$$
(2)
where.
\({C_1}=\hbox{max} \left[ {0.43,\frac{\eta }{{\eta +5}}} \right],\;\eta =S\frac{k}{\varepsilon },\;S=\sqrt {2{S_{ij}}{S_{ij}}}\)
where \(\rho\) is the flow density, kg/m2; k is the turbulence kinetic energy, J; \({x_j}\) is the coordinate positions; \({u_j}\) is the velocity of\({x_j}\), m/s; \(\mu\) is the coefficient of viscosity, Pa·s; \(\varepsilon\) is the turbulent dissipation rate;\({G_k}\) is the generation of turbulence kinetic energy due to the mean velocity gradients, J; \({G_b}\)is the generation of turbulence kinetic energy due to buoyancy, J; \({Y_M}\)is the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate; \({\sigma _k}\)and\({\sigma _\varepsilon }\) are the turbulent Prandtl numbers for k and \(\varepsilon\) respectively; \({S_k}\) and \({S_\varepsilon }\) are user-defined source terms; and \({C_{1\varepsilon }}\), \({C_2}\) and \({C_{3\varepsilon }}\) are constants.
The turbulent (or eddy) viscosity \(\mu {}_{t}\) is computed by combining k and \(\varepsilon\) as follows:
$${\mu _t}=\rho {C_\mu }\frac{{{k^2}}}{\varepsilon }$$
(3)
where
$${C_\mu }=\frac{1}{{{A_0}+{A_S}\frac{{k{U^*}}}{\varepsilon }}}$$
(4)
and
$${U^*}=\sqrt {{S_{ij}}{S_{ij}}+{{\tilde {\Omega }}_{ij}}{{\tilde {\Omega }}_{ij}}}$$
(5)
\({\tilde {\Omega }_{ij}}={\Omega _{ij}} – 2{\varepsilon _{ijk}}{\omega _k}\)
\({\Omega _{ij}}=\overline {{{\Omega _{ij}}}} – {\varepsilon _{ijk}}{\omega _k}\)
where \(\overline {{{\Omega _{ij}}}}\) is the mean rate of rotation tensor viewed in a moving reference frame with the angular velocity \({\omega _k}\). The model constants \({A_0}\) and \({A_S}\) are given by
$${A_0}=4.04,{A_S}=\sqrt 6 \cos \varphi$$
(6)
where
$$\begin{gathered} \varphi = \frac{1}{3}\cos ^{{ – 1}} (\sqrt 6 W),W = \frac{{S_{{ij}} S_{{jk}} S_{{ki}} }}{{\tilde{S}^{3} }}, \hfill \\ \tilde{S} = \sqrt {S_{{ij}} S_{{ij}} } ,S_{{ij}} = \frac{1}{2}\left( {\frac{{\partial u_{j} }}{{\partial x_{i} }} + \frac{{\partial u_{i} }}{{\partial x_{j} }}} \right) \hfill \\ \end{gathered}$$
(7)
The model constants \({C_{1\varepsilon }}\),\({C_2}\) ,\({\sigma _k}\) ,\({\sigma _\varepsilon }\)and\({\rho _\varepsilon }\)have the following default values:
\({C_{1\varepsilon }}=1.44,\;{C_2}=1.9,\;{\sigma _k}=1.0,\;{\sigma _\varepsilon }=1.2.\)
Multiphase flow model
In the field of fluid dynamics, the Eulerian, mixture, and VOF (Volume of Fluid) models are extensively employed multiphase-flow models. Eulerian models are predominantly applied to address complex problems such as those in bubble columns and updrafts39. Considering the relative simplicity of the research objectives in the present study, either the mixture model or the VOF model will be selected. It should be noted that the mixture model permits inter – phase interpenetration. Since this study involves the interpenetration between the liquid – phase refrigerant and the gas – phase refrigerant, the mixture model is adopted to depict the cavitation flow in an electronic expansion valve, which can be expressed as:
$$\frac{{\partial \rho }}{{\partial t}}+\frac{{\partial \rho {u_i}}}{{\partial {x_i}}}=0$$
(8)
$$\frac{{\partial \rho {u_i}}}{{\partial t}}+\frac{{\partial \rho {u_i}{u_j}}}{{\partial {x_j}}}= – \frac{{\partial p}}{{\partial {x_i}}}+\left[ {\left( {\mu +{\mu _t}} \right)+\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}}+\frac{{\partial {u_j}}}{{\partial {x_i}}} – \frac{2}{3}{\delta _{ij}}\frac{{\partial {u_k}}}{{\partial {x_k}}}} \right)} \right]$$
(9)
where δij is the Kronecker tensor symbol. The mixture density of the multiphase flow model ρ is:
$$\rho =\alpha {\rho _v}+\left( {1 – \alpha } \right){\rho _l}$$
(10)
where ρv and ρl are the vapor and liquid density respectively, kg/m³.The mixture viscosity of the multiphase flow model µ is:
$$\mu =\alpha {\mu _v}+\left( {1 – \alpha } \right){\mu _l}$$
(11)
Where µv and µl are the vapor and liquid kinetic viscosity respectively, Pa·s. The vapor volume fraction α is computed by: \(\alpha =f\frac{\rho }{{{\rho _v}}}\).
where f is the vapor mass fraction.
where \(\alpha\) is the vapor volume fraction; \({\rho _v}\) is the vapor density, kg/m2; \(\overrightarrow {{V_v}}\) is the vapor phase velocity, m/s; Re is the vapor generation rate; Rc is the vapor condensation rate; \({\Re _B}\) is the bubble radius, m; n is the bubble number densities; P is the liquid pressure, Pa; \({P_v}\) is the saturated vapor pressure at that state, Pa; Fvap and Fcond are evaporation and condensation correction factors, respectively; \({\rho _v}\) and \({\rho _l}\) a is the vapor density and liquid density, respectively, kg/m3.
Cavitation model
Cavitation occurs when the static pressure of the fluid drops below the vapor pressure at a certain temperature. After a comprehensive comparison between the Schnerr–Sauer cavitation model and the Zwart–Gerber–Belamri cavitation model, the Schnerr–Sauer model is chosen due to its proficiency in accurately depicting the cavitation phenomenon. This selection is based on the fact that the Schnerr–Sauer model can effectively account for the phase-transition process taking place between the gas and liquid phases. Additionally, it takes into consideration the influence of turbulent pressure fluctuations on the cavitation pressure40. The corresponding mathematical control equations are as follows10:
$$\frac{\partial }{{\partial t}}\left( {\alpha {\rho _v}} \right)+\nabla \cdot \left( {\alpha {\rho _v}\overrightarrow {{V_v}} } \right)={R_e} – {R_c}$$
(12)
$${\Re _B}={\left( {\frac{\alpha }{{1 – \alpha }}\frac{3}{{4\pi }}\frac{1}{n}} \right)^{\frac{1}{3}}}$$
(13)
when P ≤ Pv,
$${R_e}={F_{vap}}\frac{{{\rho _v}{\rho _l}}}{\rho }\alpha \left( {1 – \alpha } \right)\frac{3}{{{\Re _B}}}\sqrt {\frac{2}{3}\frac{{\left( {{P_v} – P} \right)}}{{{\rho _l}}}}$$
(14)
when P ≥ Pv,
$$R_{c} = F_{{cond}} \frac{{\rho _{v} \rho _{l} }}{\rho }\alpha \left( {1 – \alpha } \right)\frac{3}{{\Re _{B} }}\sqrt {\frac{2}{3}\frac{{\left( {P – P_{v} } \right)}}{{\rho _{l} }}}$$
(15)
Noise model
A remarkable characteristic of the broadband noise source modeling lies in its substantial improvement in accuracy when a substantial number of noise sources coexist. In the domain of fluid dynamics research, particularly when taking into account that cavitation noise predominantly exists within the broadband spectrum, the broadband noise model (encompassing broadband noise sources) is opted for conducting in – depth noise analysis of the cavitation noise generated in electronic expansion valves. This selection is of utmost importance for precisely characterizing the acoustic properties related to cavitation phenomena. The governing equations of this model are presented as follows34:
$${P_A}=\lambda {\rho _0}\left( {\frac{{{U^3}}}{l}} \right)\frac{{{U^5}}}{{\alpha _{0}^{5}}}$$
(16)
where \({P_A}\) is the sound power, W/m2; \(\lambda\) is the model constant; U is the turbulent velocity, m/s; l is the turbulent characteristic length, m; and \({\alpha _0}\) is the sound velocity, m/s.
The focus of the study
In the study, modeling of electronic expansion valve for air conditioning was carried out using 3D modeling soft SolidWorks 2018, as depicted in (Fig. 1a). The valve opening is determined by the position of the valve stem, which is influenced by the spring and fluid forces. The up‒and‒down movement of the valve stem serves to open or close the electronic expansion valve. A physical model of the valve was constructed using ANSYS SpaceClaim 2020 modeling software, based on the existing dimensions of the valve. The fluid domain within the valve was derived through reverse modeling, as illustrated in (Fig. 1b).
The geometry model of the electronic expansion valve. (a) the three-dimensional geometric section (b) the computational fluid domain.
To ensure uniform inlet flow rate distribution and mitigate the impact of backflow on calculation results, the horizontal and vertical lengths of the inlet and outlet are appropriately extended by ANSYS SpaceClaim 2020, as illustrated in (Fig. 2). Generally, the horizontal and vertical dimensions should be at least three times the diameter of the flow channel to enhance the accuracy of flow field calculations. When fluid flows from the left end of the pipeline into the lower end and exits with a positive flow direction, the air conditioning system is in a cooling condition. Conversely, when fluid flows from the lower end of the pipeline into the left end and exits with a reverse flow direction, the air conditioning system is in a heating condition.
Electronic expansion valve computational fluid domain.
The opening pulse interval of the electronic expansion valve ranges from 0 to 500. When the number of pulses is 0, the valve is in a closed state with no refrigerant flow. Conversely, when the number of pulses reaches 500, the electronic expansion valve is fully open. To avoid using an overly abstract parameter like pulses and to simplify calculations related to fluid flow cross-section, we use the rising height of the spool as an alternative measure for valve opening. The opening of the electronic expansion valve at various pulse numbers was determined by using a pulse number of 0 as a reference point to describe the valve rise height, as presented in (Table 1).
Simulation parameter setting
Due to its good refrigeration performance, R134a refrigerant is a common refrigerant in air conditioning. Therefore, this study used R134a as the refrigerant. Establish monitoring points at the inlet and outlet of the valve core during numerical calculations to evaluate the convergence of the calculation results. The location of the monitoring points is shown in (Fig. 3). Given that the electronic expansion valve possesses a stable structure and its internal flow dynamics are solely determined by the inlet and outlet boundary conditions, pressure‒inlet and pressure‒outlet conditions are employed for the calculations while assuming an adiabatic wall surface and no-slip condition.
The location of monitoring points.
In the numerical simulation process, an initial valve opening is set, and the valve is driven to open and close by a user-defined function (UDF) in the software. The calculation is carried out in a transient state, with a time step of 10 − 6 s. The inlet and outlet boundary conditions are set as the inlet and outlet pressures, respectively. The inlet is the liquid phase R134a refrigerant at 5 °C. As the refrigerant flows through the valve, the flow cross-sectional area decreases while the flow rate remains unchanged, causing the flow velocity to increase abruptly and the local pressure of the liquid to decrease. This results in cavitation and the formation of a two-phase flow, leading to a gas-liquid mixture of media being discharged from the electronic expansion valve’s outlet. The physical parameters of the 5 °C R134a refrigerant are shown in (Table 2), which are taken from the physical properties parameter table of the R134a refrigerant, and the refrigerant purity is one of the main factors affecting experimental uncertainty.
Dynamic grid setting
In numerical simulation studies, dynamic grid models are frequently employed to simulate fluid flow scenarios where the shape of the flow field boundary varies over time due to boundary movement. In an electronic expansion valve, a combined force from the fluid and the spring imparts an upward or downward force on the spool. Consequently, the boundary layer moves in response to this force, either upward or downward. Once the forces are balanced, and there is no longer any upward or downward velocity, the flow field stabilizes. The force acting on the spool is illustrated in (Fig. 4).
Schematic diagram of the spool force for the electronic expansion valve.
The spring force is expressed as F1 = K·Δx, where K is a constant of 6 N/mm. The value of Δx is related to the length of the spring in its natural state, which can be measured as x1 = 7.52 mm in the actual model. When the electronic expansion valve is not in operation, the valve port is completely closed, and the spring is in a state of compression with a length of x2 = 5.3 mm, resulting in a spring deformation of Δx = 2.22 mm.
When fluid flows through the valve, it exerts a force on the spool that opposes the spring’s force. Under the combined action of these forces, the spring deformation changes. The force acting on the spool can be calculated as Force = F2–F1, where F2 is the force exerted by the fluid and F1 is the spring force. If Force > 0, the valve stem moves upward, causing an increase in the deformation variable Δx. If Force < 0, the valve stem moves downward, resulting in a decrease in Δx. If Force = 0, the valve stem remains stationary, and Δx remains unchanged. Once the deformation variable Δx no longer changes, the valve opening reaches a stable state, and the internal flow field of the valve stabilizes.
If the entire computational fluid domain of the electronic expansion valve is taken as a whole, then the spool will be accompanied by grid reconstruction during the rising process, and the pulling action of the grid nodes will lead to the deformation of the sharp corners in the red box in (Fig. 5). This irregularity may lead to changes in the internal flow field, particularly the pressure distribution, and prolong the grid reconstruction time. Additionally, an excessive number of fine grids can introduce negative volumes. To address these issues, a straightforward enhancement to the valve model involves adding an interface connecting the surface between the vertical and horizontal motion faces, as depicted in (Fig. 6). This modification ensures both computational speed and accuracy during grid reconstruction.
Some of the custom functions are shown below:
\({\text{force}}\,=\,{\text{f}}\_{\text{glob}}\left[ {\text{1}} \right] – {\text{f}}\_{\text{spring}};\)
\({\text{dv}}\,=\,{\text{dtime}}*{\text{force}}/{\text{m}};\)
\({\text{vel}}\left[ {\text{y}} \right]\,=\,{\text{v}}\_{\text{prev}}\,+\,=\,{\text{dv}};\)
where force resultant force, f_glob[1] is extracted fluid force, f_spring is the elasticity of the spring, vel[y] is valve core up and down movement speed.
Grids and independency verification
The study employs ANSYS FLUENT 2020 software to perform numerical simulations of the cavitation characteristics and fluidic noise associated with electronic expansion valves. ANSYS FLUENT 2020 utilizes the finite volume method for numerical calculations, necessitating the division of the fluid calculation domain into a finite number of grids. The parameters of these cells are then replaced with parameter values at a specific point, and calculations are solved for all grids in a logical sequence. Prior to the numerical simulation, it is essential to divide the fluid calculation domain into a finite number of grid cells and replace the parameters within them with parameter values at a certain point. All grids must then be solved according to a specific logical order, as the number of grids and their density are directly related to the calculation results. Due to the intricate structure of the electronic expansion valve’s throat and its small volume, which necessitates high computational requirements, a tetrahedral grid is selected. Grid encryption is applied to the orifice position and the cavitation area around the spool of the flow path based on the refrigerant’s flow characteristics to ensure accurate computations. To verify whether the number of grids meets the calculation accuracy requirements, the grid is performed according to the above criteria, and grid encryption is applied in the throat area. There are five grid numbers modeled as 1.75 × 10 5, 3.72 × 10 5, 6.02 × 10 5, 9.38 × 10 5, and 1.31 × 10 6. Figure 7 shows the number of grids of 1.75 × 10 5 and 1.31 × 10 6.
Deformation of the computational domain during grid reconstruction.
Interface between the vertical and horizontal motion face.
Schematic diagram of the number of grids. (a) 1.75 × 105, (b) 1.75 × 105.
Grid division is a pivotal factor in numerical simulations, with the number of grids directly impacting computational accuracy and determining the convergence and correctness of residual results. Consequently, to guarantee calculation precision and enhance efficiency by reducing computation time, it is imperative to assess the grid density’s relevance prior to the commencement of calculations. The grid independence test was conducted for a steady simulation, utilizing constant inlet and outlet pressure boundary conditions of 1.5 MPa and 0.7 MPa, respectively, where the opening of the valve is 30%. The mass flow rate through the outlet and acoustic pressure level 15 mm horizontally from the post-throttling monitoring point are employed as a criterion for verifying grid independence. Figure 8 shows the comparison results of the mass flow rate and acoustic pressure level under different grid quantities. The mass flow rate and acoustic pressure level decreas gradually as the quantity of the grid increases. When the grid quantity is increased to 6.02 × 105 with the first boundary layer set to 5 × 10− 5 m, the error of the mass flow relative to the maximum grid quantity (1.31 × 106) is 0.84%, and the error of the acoustic pressure level relative to the maximum grid quantity (1.31 × 106) is 0.97%, which has been reduced to less than 1%. Moreover, the maximum y + value calculated using this setting is 28.6, which satisfies the conditions of the standard wall function, with the y + value required to be in the range of 0–50.
The quality of the grids is directly proportional to the accuracy of the calculation. However, an excessive number of grids can impede the computational process. Considering both the mass flow rate and y + value, a grid scheme with a grid quantity of 6.02 × 105 and a first boundary layer of 5 × 10− 5m is chosen.
The comparison results under different grid quantity. (a) The mass flow rate through the outlet. (b) The acoustic pressure level 15 mm from the post-throttling monitoring point.
Experimental platform
To validate the realism of the simulation, a noise measurement platform was constructed, and the operational noise was measured. The schematic diagram of the experiment platform is shown in (Fig. 9). The individual components in (Fig. 16) are as follows: 1-Compressor; 2-Oil separator; 3-Temperature sensor; 4-Pressure sensor; 5-Condenser; 6-Reservoir; 7-Subcooler; 8-Filter; 9-Mass flow meter; 10-Front liquid sight; 11-Ball valve; 12-Microphone; 13-Electric expansion valve; 14- Rear liquid sight; 15- Evaporator; 16- Gas‒liquid separator (Zhou et al., 2023).
The schematic of the experimental platform for measuring the flow-induced noise of the electronic expansion valve.
Experimental comparison and analysis
To ensure the accuracy of the results, validation of the current numerical model was essential. Comparative analyses between experimental and simulated flow rate tests were performed for valve openings ranging from 10 to 80%, with increments of 10%. The inlet and outlet pressures were maintained at 1.5 MPa and 0.7 MPa, respectively. The inlet flow rate was set to 0.105 kg/s, and the time step was set to 10− 6 s. The results are illustrated in (Fig. 10). As observed from the figure, the flow rates obtained by the simulation are in good agreement with the experimental data. The maximum error is less than 5%. It has been verified that the numerical model used in this study is acceptable.
Comparison between numerical results and experimental data for different openings of the valve.
Furthermore, uncertainties in testing can lead to discrepancies between test outcomes and simulation values. The purity of the refrigerant plays a crucial role in determining test results, as impurities can impact the air conditioning system’s flow and noise levels. Additionally, the noise insulation of the test chamber is a pivotal factor that influences test outcomes, while poor sealing can compromise data authenticity. Test equipment accuracy and sensitivity are also essential factors that affect test results; using low-frequency and accurate equipment can result in significant errors.