In this section, the numerical simulation results together with the experimental results are presented. First, the validity of steady numerical solvers was confirmed by comparing the mass flow rate obtained from the compressible RANS simulations of the axisymmetric critical nozzle with the measured ones. Then, the compressible steady-state RANS simulation was performed on the 3D critical nozzle, and the unsteady compressible LES was performed using the RANS results as an initial condition. The validity of the LES was confirmed by comparing the predicted surface pressure fluctuations with the measured ones. Finally, the wavenumber-frequency analysis was carried out and the acoustic power due to the critical nozzle flow was assessed.
4.1. Two-Dimensional Axisymmetric Steady Simulation
The compressible RANS simulation was performed for a two-dimensional axisymmetric critical nozzle and the accuracy of the simulation was verified through the comparison of the predicted mass flow rate with measured data. The critical nozzle targeted is the toroidal-throat Venturi nozzle as specified in ISO 9300. Its geometry and the numerical mesh with the specified boundary conditions are shown in
Figure 5.
The simulations were carried out for the two types of nozzle of which the throat diameter was 0.4 mm and 0.56 mm, respectively. The inlet and outlet pressure were specified on the inlet and outlet boundaries, respectively, the adiabatic and no-slip wall conditions were set on the wall, and the axisymmetric condition was applied along the centerline axis of the venturi nozzle. The grids were constructed by using a structured grid system, and y
+ was set to be around five, and the total number of grids was 350 × 100. The efficiency of the critical nozzle can be determined by observing the variation of mass flow rate according to the change of pressure difference between the inlet and the outlet. In this study, the pressure of the outlet was varied to control the pressure difference between the inlet and outlet of the nozzles. The applied inlet and outlet pressure conditions are shown in
Table 1 [
16].
Figure 6 shows the predicted iso-contours of the Mach number inside the nozzles of d = 0.4 m at the pressure ratio of 0.666 and d = 0.56 m at the pressure ratio of 0.778. It is seen that shock waves are generated around the throat in both cases, which implies that the choked flow appears.
Figure 7 compares the predicted ones with the measured data in terms of the variation of the non-dimensional mass flow rate according to the pressure ratio between the inlet and outlet. The
x-axis represents the ratio of the outlet static pressure to the total inlet pressure, and the
y-axis represents the ratio of the mass flow rate(
qm) to the mass flow rate(
qm, ref) when the outlet pressure is equal to the atmospheric pressure. In the case of the nozzle flow, as the pressure difference between the inlet and outlet becomes larger, the flow is more accelerated up to the nozzle throat and the mass flow rate increases. However, when the outlet pressure becomes equal to or less than the critical pressure at which the chocked flow starts to form, the Mach number in the throat reaches one and does not exceed one anymore, and thus the mass flow rate is kept constant. In
Figure 6, because the outlet pressure in both nozzles is lower than the critical pressure, the Mach number is one in the nozzle throat. Then the flow is accelerated to supersonic in the diffusion region through the nozzle throat, and shock waves are generated. This phenomenon corresponds to a typical process for creating a supersonic flow in an internal flow. In
Figure 7, it can be seen that the mass flow rate through the nozzle increases as the outlet pressure decreases to the critical pressure. As the pressure ratio is lower than the critical ratio, the choked flow occurs and the mass flow rate is kept constant. The numerical results are in excellent agreement with the experimental results [
16].
4.2. Three-Dimensional Steady Simulation
In the pipe system used in the offshore plant, high-temperature and high-pressure gas is transferred, and a pressure relief valve is used to prevent damage to the pipe which may be caused by strong pressure waves. However, the sudden pressure variation in the pressure relief valve emits high acoustic energy in the downstream direction of the pipe. To simulate this sudden pressure variation, the numerical simulations using the prescribed numerical methods were performed on the three-dimensional critical nozzle as shown in
Figure 8. The computation grids are shown in
Figure 9, and the boundary conditions were set to be the same as those used in the preceding two-dimensional axisymmetric problem. The grids used for the flow simulation were structured grids, the number of grids was 6,470,000, and the diameter of the nozzle throat was 8.5 mm. The overall grid quality was confirmed using the commercial grid generation tool ICEM CFD and the grid quality index was greater than 0.61, which confirms the good quality of the grids. In addition, y
+ based on the inflow condition was kept below 1. To construct the initial condition for the unsteady LES, three-dimensional steady-state compressible RANS simulation was performed beforehand.
Figure 10 presents the numerical results of pressure and temperature distributions inside the pipe and the nozzle. It can be seen that the shock wave is formed inside the critical nozzle due to the high pressure difference in the pipe inlet and outlet, and the sudden pressure variation occurs before and after the shock wave. The shock inside the nozzle meets the boundary layer of the wall surface and the boundary layer separation occurs at this location. This complex interaction generates a strong coherent turbulent flow which in turn causes intense aerodynamic noise. The corresponding unsteady phenomena can be identified in
Figure 11, which is discussed in detail in the next section. The similar phenomenon was also observed and discussed in the preceding study by Kim et al. [
8].
4.3. Three-Dimensional Unsteady Simulation
LES is preferred for high accuracy unsteady simulation but requires many grids and a long simulation time. Also, it requires proper initial conditions for stable numerical simulation. Therefore, three-dimensional unsteady LES was performed with the initial condition of the preceding RANS simulation results. The total number of the grid was the same as the RANS simulation. The simulation time interval was set to be 5 × 10
−5 s.
Figure 11 shows the iso-contours of pressure and vorticity magnitudes in the zoomed plot as well as in the entire computational domain. It can be found that a sudden pressure drop occurs at the nozzle throat, and strong vortices are formed downstream just behind the throat, which causes a large amount of fluctuating turbulent disturbance components to be convected in the downstream direction of the pipe. This turbulent disturbance generated acoustic pressure inside the pipe.
In the ISO 9300 standards, it is required to measure the pressure at the positions 2D and 4D downstream away from the end of the critical nozzle, where D is the pipe diameter. The pressure fluctuations at the same locations denoted as Sensor 1 and Sensor 2 in
Figure 3 and
Figure 11 are measured using the experimental setup described in
Section 3. The experimental devices and methods follow the ISO 9300 and the experimental conditions are the same as the simulation conditions.
Figure 12 compares the 1/3 octave band spectrum of pressure between the measured and the predicted ones at each sensor position. There is good agreement between two results at the position of Sensor 2 while there is a relatively large discrepancy between two results at the position of Sensor 1 where the numerical results over-predict in comparison with the measured data. This seems to be due to the fact that the Sensor 1 is located in the region where the strong vortex is observed in the simulation while the location of Sensor 2 is far enough away from the strong turbulence region. However, the predicted pressure levels at both locations become much lower than the measured one in the high-frequency region above 2000 Hz, which seems to be due to poor resolution of the grids in the corresponding frequency range. More precise numerical simulations can be performed by using more resolved grids and fewer time intervals but at a much more expensive numerical cost. Since the AIV phenomenon is known to occur in the frequency range of 500 Hz to 2000 Hz, the current numerical resolution was determined to be sufficient for that purpose.
4.4. Wavenumber-Frequency Analysis
Based on the unsteady LES results, the wavenumber-frequency analysis was performed by using the Fourier transform in the cylindrical coordinate system as described in
Section 2.4 with the in-house code, of which the validity was confirmed in the previous studies [
8,
13]. The acoustic power caused by the sudden pressure variation inside the nozzle provides the source of the AIV of the pipe wall. Therefore, for more precise assessment of AIV problems, it is necessary to extract the acoustic wave components propagating in the downstream direction from the fluctuating pressure field in the pipe. Therefore, the pressure extraction area is selected to include the prescribed monitoring position as shown in
Figure 13.
The simulation time interval and the frequency interval used for the wavenumber-frequency analysis were 5 × 10
−5 s and 3.23 Hz, respectively, and the grid spacing was (
Δθ,
Δz) = (2.5°, 0.0015 m). The average speed of sound and Mach number in each region were 343 m/s, 0.064 in Region 1 and 0.011 in Region 2. The results of the wavenumber-frequency analysis are shown in
Figure 14.
In
Figure 14, the black lines represent the modal cut-off lines in the pipe, and the circumferential direction mode is represented by
m and the radial direction mode is represented by
n. The dashed line at the right bottom represents the wave components propagating at the constant speed V = 21.91 m/s in Region 1 and V = 3.87 m/s in Region 2, which is the incompressible pressure components convected at the mean flow velocity in the downstream direction and associated with turbulence fluctuation acting as an aerodynamic noise source. Since the acoustic wave generated by the aerodynamic noise sources propagates at the speed of sound, the compressible pressure fluctuation (plane wave components in the pipe) propagates in the downstream direction at the speed of
, which is the sum of the sound velocity and the flow velocity, and at the speed of
in the upstream direction. In the case of other non-planar waves, the phase speed is higher than
and
in the upstream and downstream directions, respectively. Therefore, the pressure wave components located in the region between the two lines of which the phase speed is
and
denotes those propagating at a speed higher than the speed of sound and thus are compressible acoustic components. The cut-off line for each duct mode can be defined as Equation (10) and is depicted as a curve in
Figure 14. In the case of the incompressible pressure fluctuation, the components propagating at the average mean flow velocity are not clearly visible, which is different from the results shown in the study by Kim et al. [
8]. The reason for this is because the vortices of various sizes are formed and propagate at various speeds due to an abrupt change in the cross-sectional area between the nozzle and the downstream pipe. In the case of the compressible pressure fluctuation component, it can be seen that the spectral density is prominently distributed in the vicinity of the cut-off line of each mode and becomes gradually weaker as the order becomes higher. As mentioned above, the predicted pressure level assures up to 2000 Hz, and high grid resolution is required for higher frequency range. Therefore, the high-frequency components above 2000 Hz are rapidly attenuated.
The application of inverse wavenumber-frequency transform on the region of wavenumber-frequency diagram including only the compressible pressure fluctuations leads to the acoustic power spectral density, eliminating the incompressible pressure fluctuations.
The acoustic power spectrum of the internal aerodynamic noise obtained in this way is shown in
Figure 15. It shows the power spectrum of three components: the incompressible pressure fluctuation part due to turbulence fluctuation, the compressible pressure fluctuation part due to acoustic waves, and the total pressure fluctuation including the incompressible and compressible parts. From the spectrum and overall values of all parts, the power level in Region 1, which is the region more upstream, is greater than in Region 2. The influence of the incompressible part due to flow is dominant in the low frequency region below 1000 Hz, but the influence of the compressible part due to the acoustic wave is dominant in the region of 1000 Hz to 2000 Hz. However, in the higher frequency range, it is rapidly dissipated due to poor time and spatial resolution. It can be seen that the incompressible part becomes weaker from Region 1 to Region 2 because the vortex rapidly dissipates as shown in
Figure 11b.
For the incompressible part related to the turbulence fluctuation, during the convection of vortical waves downstream, its energy is transferred to smaller-scale eddies oscillating with the high frequency, that is so-called energy cascade phenomena occurring in the inertial subrange. The decay rate of its spectrum is known to follow the −5/3 power law. The even smaller eddy scale, known as the Kolmogorov Scale, belongs to the dissipation range where energy is dissipated due to heat or sound, and is known to follow the −7 power law [
17]. It can be found from
Figure 15 that, in both Region 1 and Region 2, the inertial subrange exists in the region between 100 Hz and 1000 Hz, and a dissipation range exists in the frequency region higher than 1000 Hz. Especially in Region 1, where larger coherent turbulent fluctuations exist and are followed by its energy cascade phenomena, the inertial subrange and dissipation range are more prominent and distinct. Note that the decay rate of the spectrum in the frequency range higher than around 2000 Hz to 3000 Hz is due to numerical damping. The good agreement of the decay rate of the predicted spectrum confirm the validity of the LES simulation and the wavenumber-frequency analysis. It is noteworthy that the magnitude of compressible components is higher than that of incompressible ones in the frequency range belonging to the dissipation rage. This fact ensures that the energy of turbulent fluctuations is transferred to the compressible acoustic waves in the dissipation range. It can be also seen from the overall power level that both incompressible and compressible components are attenuated downstream, but the incompressible part is further reduced compared to the compressible part due to more rapid dissipation of the vortex.
Overall, these results show that the current numerical methods predict the acoustic power due to the valve flow in a pipe more precisely than the existing ones relying on the quasi-empirical equations by Carucci and Muller, considering more complex geometry and conditions.