Wavenumber-Frequency Analysis of Internal Aerodynamic Noise in Constriction-Expansion Pipe

High-pressure gas is produced during the oil production process at offshore plants, and pressure relief devices, such as valves, are widely used to protect related systems from it. The high-pressure gas in the pipes connected to the flare head is burned at the flare stack, or, if it is nontoxic, is vented to the atmosphere. During this process, excessive noise is generated by the pressure relief valves that are used to quickly discharge the high-pressure gas to the atmosphere. This noise sometimes causes severe acoustic-induced vibration in the pipe wall. This study estimated the internal aerodynamic noise due to valve flow in a simple constriction-expansion pipe, by combining the large eddy simulation technique with a wavenumber-frequency analysis, which made it possible to decompose the fluctuating pressure into the incompressible hydrodynamic pressure and compressible acoustic pressure. First, the steady-state flow was numerically simulated, and the result was compared with a quasi-one-dimensional theoretical solution, which confirmed the validity of the current numerical method. Then, an unsteady simulation analysis was performed to predict the fluctuating pressure inside a pipe. Finally, the acoustic pressure modes in a pipe were extracted by applying the wavenumber-frequency transform to the total pressure field. The results showed that the acoustic pressure fluctuations in a pipe could be separated from the incompressible ones. This made it possible to obtain accurate information about the acoustic power, which could be used to assess the likelihood of a piping system failure due to acoustic-induced vibration, along with information about the acoustic power spectrum of each acoustic mode, which could be used to facilitate the systematic mitigation of the potential acoustic-induced vibration in piping systems.


Introduction
High-temperature and high-pressure gas is produced by various processes in chemical, petrochemical, and offshore oil production plants.Excessive noise and vibration can be generated when such gases are discharged, usually as they pass through a valve in the piping system used to protect the system and equipment from being damaged by such high-pressure gases.In the 1980s, pipe failures due to acoustic-induced vibration (AIV) were reported.Thus, the assessment and control of noise and vibration in piping systems have become essential parts of the design stage.Carucci et al. [1] first investigated pipe failure using nine pipe failure cases and 27 normal pipe cases.They suggested a criteria curve for the sound power and pipe diameter.This curve has been used as the basic screening method for the AIV of pipes.Eisinger [2] proposed simple piping design criteria by including the pipe wall thickness as an additional parameter.These were adopted in the Norsok Standard [3].Recently, Bruce et al. [4] suggested new guidelines for the pipe AIV problem, in which the maximum vibration level is proportional to a nondimensional parameter that is given by the pipe diameter divided by the pipe wall thickness.The Energy Institute [5] published guidelines for the pipe AIV problem based on the likelihood of failure (LOF), which depends on parameters, such as the branch connection, main pipe diameter ratio, and material type.
However, all of these studies assumed that the main noise sources in a piping system are the pressure relief valves and estimated their acoustic power levels using an empirical formula that was based on API521 [6] and IEC 60534-8-3 [7].In addition, they used an empirical criteria curve, the main parameters of which were the sound power and pipe dimensions obtained from the available pipe failure cases.Although these methods are very efficient at quickly screening a large number of pipes, it is difficult to apply efficient mitigation measures, because these tools rely on a limited number of input parameters to quantify the LOF, without any understanding of the relevant physical mechanisms.Recently, Agar et al. [8] proposed a finite element methodology based on fluid-structure coupling, which makes it possible to predict the dynamic stress in a complex piping system.However, even though this method could be used to quantitatively and effectively assess mitigation actions, the acoustic energy that was generated by a pressure-reducing device was estimated using a simple empirical formula, which inevitably limited its predictive accuracy with significant uncertainty.
In this study, to overcome this limitation, a systematic first-principles-based method was developed by combining high-order computational fluid dynamics (CFD) techniques with a wavenumber-frequency analysis.This method makes it possible to separate the acoustic pressure waves from the total pressure fluctuations in a pipe flow induced by valves, and to compute the acoustic power that is generated by the valves.First, the valve flow is predicted using high-resolution CFD techniques, in which the valve is modeled as a simple constriction-expansion pipe for simplicity.The validity of the calculated steady pressure is confirmed by comparing it to the theoretical solution obtained from quasi-one-dimensional (1D) inviscid compressible flow equations.Then, an unsteady simulation is performed to predict the fluctuating total pressure inside a pipe.Finally, the acoustic pressure modes in the pipe are extracted by applying the wavenumber-frequency transform to the total pressure field.The results show that the acoustic pressure fluctuations can be separated from the incompressible ones; this can provide accurate information about the acoustic power spectrum, which is used as an important parameter to determine the LOF of a piping system due to AIV.
The remainder of this paper is organized as follows.Section 2 introduces the theoretical background related to valve noise in piping systems, including the wavenumber-frequency analysis.Section 3 describes the governing equations and the numerical methods in detail.Section 4 presents the numerical results, including the steady and unsteady flow fields.Finally, the acoustic pressure modes in a pipe are extracted from the total pressure field using the wavenumber-frequency analysis, which makes it possible to predict the acoustic power spectrum.Section 5 concludes this paper.

Valve Flow and Its Quasi-1D Approximation
Commercially available valves may be categorized into globe valves, butterfly valves, pressure relief valves, and so on.A schematic diagram of a typical pressure relief valve is shown in Figure 1.
A relief valve is designed to open at a set pressure to protect a pressure vessel or other equipment.When the inlet pressure exceeds this set pressure, the relief valve opens and a portion of the gas is diverted through the valve system.As the fluid begins to move faster, its compressibility plays an increasingly important role, and it undergoes a significant change in density.If the pressure ratio across the constriction area exceeds a critical value, a choked flow with shock waves occurs [9], which causes additional aerodynamic noise.
The flow passage through a valve can be considered as a simple constriction and expansion pipe in terms of its cross-sectional area variation, whereby the area is reduced to a minimum at the throat of the valve and is then increased.An orifice plate works by essentially the same principle [10], in that it reduces the local pressure by providing a mechanical constriction in the pipe, although a valve is much more sophisticated than an orifice.
Appl.Sci.2017, 7, 1137 3 of 16 throat of the valve and is then increased.An orifice plate works by essentially the same principle [10], in that it reduces the local pressure by providing a mechanical constriction in the pipe, although a valve is much more sophisticated than an orifice.The high-speed flow through a pipe with an orifice can be simplified as a quasi-1D compressible flow if the cross-sectional area is smoothly varied: the velocity components normal to the axial direction of the pipe (in the x and y directions) are small when compared to the velocity in the axial (z) direction, such that the flow field can be assumed to be a function of z only.In other words, it can be assumed that the flow is uniform across the cross-sectional plane of the pipe, which is called a quasi-1D flow.In this quasi one-dimensional flow, internal flow can be predicted analytically by solving the relevant equations presented by White [9].
As the pressure difference between the upstream and downstream regions increases, the mass flow rate also increases until the system reaches a critical condition, under which the mass flow no longer increases with an increase in the pressure difference.If the pressure difference increases further, the flow velocity in the throat remains constant at Mach 1.This is also called a choked condition [11].When the outlet pressure is below the critical pressure, normal shock waves occur in the pipe, which are generated perpendicular to the flow direction.Generally, normal shock waves occur in a very thin region on the scale of 10 −7 m, and are characterized by sudden changes in the flow quantities [12].In a steady inviscid adiabatic 1D normal shock wave problem, the flow quantities downstream of the shock wave can be calculated from the upstream flow properties, which are used as the inlet conditions.The changes in the pressure and velocity across the shock wave are expressed as follows:   where p1 and p2 are the pressures before and after shock wave, respectively.Because the shock wave is an adiabatic and irreversible process, an energy loss occurs and the entropy increases.The entropy variation across the shock wave is expressed as follows: The high-speed flow through a pipe with an orifice can be simplified as a quasi-1D compressible flow if the cross-sectional area is smoothly varied: the velocity components normal to the axial direction of the pipe (in the x and y directions) are small when compared to the velocity in the axial (z) direction, such that the flow field can be assumed to be a function of z only.In other words, it can be assumed that the flow is uniform across the cross-sectional plane of the pipe, which is called a quasi-1D flow.In this quasi one-dimensional flow, internal flow can be predicted analytically by solving the relevant equations presented by White [9].
As the pressure difference between the upstream and downstream regions increases, the mass flow rate also increases until the system reaches a critical condition, under which the mass flow no longer increases with an increase in the pressure difference.If the pressure difference increases further, the flow velocity in the throat remains constant at Mach 1.This is also called a choked condition [11].When the outlet pressure is below the critical pressure, normal shock waves occur in the pipe, which are generated perpendicular to the flow direction.Generally, normal shock waves occur in a very thin region on the scale of 10 −7 m, and are characterized by sudden changes in the flow quantities [12].In a steady inviscid adiabatic 1D normal shock wave problem, the flow quantities downstream of the shock wave can be calculated from the upstream flow properties, which are used as the inlet conditions.The changes in the pressure and velocity across the shock wave are expressed as follows: where p 1 and p 2 are the pressures before and after shock wave, respectively.Because the shock wave is an adiabatic and irreversible process, an energy loss occurs and the entropy increases.The entropy variation across the shock wave is expressed as follows: The entropy variation can be expressed using the total pressure change as follows: In the case of a gas with T < 1000 K, the specific heat can be assumed to be constant.A perfect gas with a specific heat constant is referred to as a calorically perfect gas.In this paper, the flow through a pipe can be assumed to be of a calorically perfect gas.Therefore, it can be assumed that the specific heat (c p ) and gas constant (R) are constant.
Because the mass flow rate through the pipe should be constant, the relationship between the cross-sectional area and Mach number in the pipe can be derived using the isentropic flow relation and mass conservation.
where A is the area of the pipe and A* is the critical area, which is the area of the throat.In this paper, because the geometry of the pipe is known, the theoretical solution derived from the above equations is compared with the numerical solution to confirm the validity of the numerical method.

Wavenumber-Frequency Analysis for Internal Aerodynamic Noise
AIV is known to originate with the high acoustic energy that is generated by pressure-reducing devices, which causes pipe shell vibration and thus induces excessive dynamic stress.Acoustic energy propagates at the speed of sound, i.e., the energy travels at the same speed as the compressible pressure waves.The total pressure fluctuation induced by valve flow in a pipe consists of compressible and incompressible pressure fluctuations.The incompressible part is the hydrodynamic pressure of a turbulent flow field with a relatively high energy content, whereas the compressible pressure fluctuation propagates at the speed of sound, which is determined by the compressibility of the fluid.The incompressible pressure fluctuation propagates at the convection speed of the mean flow.Therefore, it is essential to separate the compressible and incompressible acoustic pressure fluctuations to accurately assess the AIV.This can be achieved using a wavenumber-frequency analysis, which is described in detail below [13].
Because the cross-sectional shape of a pipe is circular, it is better to perform a wavenumber-frequency analysis in a cylindrical coordinate system.The analytical solution of the three-dimensional (3D) wave equation in a circular duct can be written in the following form: where r, θ, and z are the coordinates in the radial, circumferential, and axial directions, respectively, and t is time.J m is the 1st kind of Bessel function of the m-th order in the circumferential direction and ω is the angular frequency, k z is the wavenumber in z-direction.C is the unknown constant and the subscript + and − mean the upstream and downstream parts, respectively.Using the Fourier transform for the z-direction propagated wave, the spectral power density is expressed as follows: Appl.Sci.2017, 7, 1137 where k r,m,n is the radial wavenumber of the m-th circumferential mode and n-th radial mode, and k z is the wavenumber in the axial direction.The dispersion relation of the internal acoustic waves to a constant mean flow of Mach number M can be expressed in the following form: Thus, the axial wavenumber in the positive and negative axial directions can expressed as follows: Because the acoustic modes affecting the pipe wall vibration are the θ and z modes, the Fourier transform is considered only in the θ and z directions, whereas the coordinate in the r direction is treated as a constant.For the cut-on higher-order mode, the following relation must be satisfied.
In this paper, the discrete Fourier transform is performed using the Hanning window, which can be written in the following form:

Acoustic Power Level
The formula for the acoustic power of a fluid in a pipe is briefly introduced, because most of the existing empirical approaches consider the acoustic power level to be a critical parameter.The acoustic particle velocity in a duct with a circular cross-section can be described as follows [14]: From Equations ( 6) and ( 12), the acoustic power is defined as follows: where I z is the intensity of the propagated acoustic wave, A is the cross-sectional area, and I is the time-averaged acoustic intensity.After some algebra, the formula for the acoustic power can be written as follows: Note that this is valid when the effects of the mean flow are negligible.

Governing Equations
As the governing equations for the fluid flow in a pipe, the compressible Navier-Stokes equations are used in the following form: where ρ is the fluid density, → v is a velocity vector, p is the static pressure, = τ is the viscous stress tensor, and → F is the external body force.The stress tensor, = τ, is found as follows: where µ is the dynamic viscosity, and I is the unit tensor.The energy conservation is described by the following equation: where E is the energy per unit mass, h j is the heat transfer coefficient, and J j and S h are the diffusion flux and heat source, respectively.

Numerical Methods
First, a steady-state simulation was performed.The Reynolds-averaged Navier-Stokes equations (RANS) were used as the governing equations in the following form: Equation ( 19) could be solved using the k-ω SST turbulence model, where k represents the turbulent kinetic energy, and ω denotes the specific turbulent dissipation rate.The k-ω SST model uses the k-ω model near the wall and a k-ε model in the region far from the boundary for accurate and stable calculations.This turbulence model is known to be the most suitable model for a supersonic internal flow16.The equations are written in the following form: where G k is the production term for the turbulent kinetic energy k due to the velocity gradient, and G ω is the generation term for the specific dissipation rate ω, Γ and Y are the effective diffusivity and dissipation terms, respectively.The k-ω SST model requires a cross-diffusion term D ω , which contributes to the conversion relationship between the k-ω model and the k-ε model.S k and S ω are source terms.The density-based implicit solver was used to solve the RANS equation numerically.
The flow was discretized using a second-order upwind scheme, and the turbulent kinetic energy and dissipation rate were found using a first-order upwind scheme.
Then, the large eddy simulation (LES) method with a Smagorinsky-Lilly model was used for the unsteady simulation to predict the fluctuating pressure waves inside a pipe with a high resolution.The turbulent viscosity of the Smagorinsky-Lilly model was calculated as follows.
where S ij and L s are the rate of the strain tensor and the mixing length, respectively.These are expressed in the following form: where κ is the von Karman constant, d is the distance to the nearest wall, C s is the Smagorinsky constant, and ∆ is the local grid size, which is defined as follows: Although Lilly proposed C s = 0.23 for a homogeneous isotropic turbulence, because it is known that excessive damping occurs for large pressure fluctuation components, C s ≈ 0.1 was used to obtain an accurate unsteady pressure in the current numerical simulation.The second-order upwind scheme was used as the spatial discretization scheme.The maximum Courant number was set to 1 and the time step was 0.05 ms, for which the Nyquist frequency was 10,000 Hz.The nonreflective boundary was applied to the outlet to accurately capture the acoustic pressure waves in the computation domain, by eliminating any possible contamination due to reflected waves from the outflow boundary.The simulations were performed with the commercial CFD code Fluent, which is based on finite volume methods.A wavenumber-frequency analysis was performed and the results were combined with the LES results to separate the compressible acoustic pressure waves from the incompressible hydrodynamic pressure fluctuations.

Numerical Results
First, the steady-state results were used to confirm the validity of the current numerical method through a comparison of its prediction result with those of the analytical solution, and then the unsteady LES results were investigated.The total pressure and temperature values were prescribed on the inlet and outlet surfaces of the pipe as boundary conditions.The boundary values of the pressure and temperature were determined based on the design values of the piping system that were obtained using the software HYSYS, which is widely used for process simulation.

Axisymmetric Steady Simulation
The valve flow was simplified as a quasi-1D flow in a pipe consisting of uniform inlet and outlet ducts with the same diameter (0.292 m), and an orifice with a minimum diameter of 0.144 m, as shown in Figure 2.These dimensions were based on the data for a real valve in a pipe that is used in an existing ocean plant.
The inlet and outlet flow conditions were prescribed using the following design data.The inlet total pressure was 9.0 bar, inlet static pressure was 7.0 bar, inlet temperature was 523 K, outlet pressure was 4.4 bar, and outlet temperature was 433 K.An axisymmetric boundary condition was applied along the center line of the pipe, and the no-slip boundary condition was applied on the inner pipe wall.The number of grid points in a typical calculation was 200,000.
The results of the axisymmetric steady simulation were compared with those that were obtained from the quasi-1D theory.Figure 3 shows the convergence time of the simulation in terms of the inlet and outlet mass flows.The mass flow rates are seen to converge quickly after 5000 iterations to the same value (25.7 kg/s) at the inlet and outlet.Because the flow is an inflow at the inlet, the mass flow rate has a positive value, and vice versa at the outlet.All of the subsequent analyses were performed using the converged numerical results.

Axisymmetric Steady Simulation
The valve flow was simplified as a quasi-1D flow in a pipe consisting of uniform inlet and outlet ducts with the same diameter (0.292 m), and an orifice with a minimum diameter of 0.144 m, as shown in Figure 2.These dimensions were based on the data for a real valve in a pipe that is used in an existing ocean plant.The inlet and outlet flow conditions were prescribed using the following design data.The inlet total pressure was 9.0 bar, inlet static pressure was 7.0 bar, inlet temperature was 523 K, outlet pressure was 4.4 bar, and outlet temperature was 433 K.An axisymmetric boundary condition was applied along the center line of the pipe, and the no-slip boundary condition was applied on the inner pipe wall.The number of grid points in a typical calculation was 200,000.
The results of the axisymmetric steady simulation were compared with those that were obtained from the quasi-1D theory.Figure 3 shows the convergence time of the simulation in terms of the inlet and outlet mass flows.The mass flow rates are seen to converge quickly after 5000 iterations to the same value (25.7 kg/s) at the inlet and outlet.Because the flow is an inflow at the inlet, the mass flow rate has a positive value, and vice versa at the outlet.All of the subsequent analyses were performed using the converged numerical results.Figure 4 shows the iso-contours of the static pressure, temperature, and velocity magnitude in the vicinity of the orifice.Their variations on the center line (y = 0) of the pipe are also shown in Figure 5, where the velocity results are described in terms of the Mach number.Note that the orifice is located at x = 1 m.Upstream from the orifice, the total pressure is approximately 9.0 bar, and it continuously decreases as the flow approaches the orifice.Downstream of the orifice, the static pressure increases sharply up to 4.4 bar and remains almost constant far downstream.It is found that the shock wave is located at approximately x = 1.1 m.Other variables show similar behaviors.To verify the numerical results, these distributions along the center line were compared with those of the theoretical solutions.To improve the understanding of the numerical solutions, three conditions  The inlet and outlet flow conditions were prescribed using the following design data.The inlet total pressure was 9.0 bar, inlet static pressure was 7.0 bar, inlet temperature was 523 K, outlet pressure was 4.4 bar, and outlet temperature was 433 K.An axisymmetric boundary condition was applied along the center line of the pipe, and the no-slip boundary condition was applied on the inner pipe wall.The number of grid points in a typical calculation was 200,000.
The results of the axisymmetric steady simulation were compared with those that were obtained from the quasi-1D theory.Figure 3 shows the convergence time of the simulation in terms of the inlet and outlet mass flows.The mass flow rates are seen to converge quickly after 5000 iterations to the same value (25.7 kg/s) at the inlet and outlet.Because the flow is an inflow at the inlet, the mass flow rate has a positive value, and vice versa at the outlet.All of the subsequent analyses were performed using the converged numerical results.Figure 4 shows the iso-contours of the static pressure, temperature, and velocity magnitude in the vicinity of the orifice.Their variations on the center line (y = 0) of the pipe are also shown in Figure 5, where the velocity results are described in terms of the Mach number.Note that the orifice is located at x = 1 m.Upstream from the orifice, the total pressure is approximately 9.0 bar, and it continuously decreases as the flow approaches the orifice.Downstream of the orifice, the static pressure increases sharply up to 4.4 bar and remains almost constant far downstream.It is found that the shock wave is located at approximately x = 1.1 m.Other variables show similar behaviors.To verify the numerical results, these distributions along the center line were compared with those of the theoretical solutions.To improve the understanding of the numerical solutions, three conditions Figure 4 shows the iso-contours of the static pressure, temperature, and velocity magnitude in the vicinity of the orifice.Their variations on the center line (y = 0) of the pipe are also shown in Figure 5, where the velocity results are described in terms of the Mach number.Note that the orifice is located at x = 1 m.Upstream from the orifice, the total pressure is approximately 9.0 bar, and it continuously decreases as the flow approaches the orifice.Downstream of the orifice, the static pressure increases sharply up to 4.4 bar and remains almost constant far downstream.It is found that the shock wave is located at approximately x = 1.1 m.Other variables show similar behaviors.To verify the numerical results, these distributions along the center line were compared with those of the theoretical solutions.To improve the understanding of the numerical solutions, three conditions for the numerical results were compared with the theoretical ones: steady with friction, steady without friction, and unsteady with friction.Because the analytical solution was based on the 1D inviscid compressible Euler equations, the numerical results without viscous friction showed the closest agreement with the theoretical ones.The steady and unsteady RANS results showed similar behaviors, except for the oscillation of the steady results in the vicinity of the shock waves.
However, there was excellent agreement between the numerical and theoretical results, except for the region just after the shock wave, where the complexity of the interaction between the shock wave and boundary layer invalidated the 1D approximation based on the theoretical approach.These results ensured the validity of the current numerical method.
for the numerical results were compared with the theoretical ones: steady with friction, steady without friction, and unsteady with friction.Because the analytical solution was based on the 1D inviscid compressible Euler equations, the numerical results without viscous friction showed the closest agreement with the theoretical ones.The steady and unsteady RANS results showed similar behaviors, except for the oscillation of the steady results in the vicinity of the shock waves.However, there was excellent agreement between the numerical and theoretical results, except for the region just after the shock wave, where the complexity of the interaction between the shock wave and boundary layer invalidated the 1D approximation based on the theoretical approach.These results ensured the validity of the current numerical method.

Three-Dimensional Steady Simulation
To obtain the initial values for the unsteady LES simulation, a 3D steady RANS simulation was performed using the same inlet and outlet boundary conditions as the axisymmetric simulation.The geometry and grids are shown in Figure 6.The number of grids used was 2,700,000.To obtain more reliable unsteady results by reducing the possible reflection waves from the outlet, the length of the computational domain was extended so that the outlet was 20.5D away from the orifice in the downstream direction.for the numerical results were compared with the theoretical ones: steady with friction, steady without friction, and unsteady with friction.Because the analytical solution was based on the 1D inviscid compressible Euler equations, the numerical results without viscous friction showed the closest agreement with the theoretical ones.The steady and unsteady RANS results showed similar behaviors, except for the oscillation of the steady results in the vicinity of the shock waves.However, there was excellent agreement between the numerical and theoretical results, except for the region just after the shock wave, where the complexity of the interaction between the shock wave and boundary layer invalidated the 1D approximation based on the theoretical approach.These results ensured the validity of the current numerical method.

Three-Dimensional Steady Simulation
To obtain the initial values for the unsteady LES simulation, a 3D steady RANS simulation was performed using the same inlet and outlet boundary conditions as the axisymmetric simulation.The geometry and grids are shown in Figure 6.The number of grids used was 2,700,000.To obtain more reliable unsteady results by reducing the possible reflection waves from the outlet, the length of the computational domain was extended so that the outlet was 20.5D away from the orifice in the downstream direction.

Three-Dimensional Steady Simulation
To obtain the initial values for the unsteady LES simulation, a 3D steady RANS simulation was performed using the same inlet and outlet boundary conditions as the axisymmetric simulation.The geometry and grids are shown in Figure 6.The number of grids used was 2,700,000.To obtain more reliable unsteady results by reducing the possible reflection waves from the outlet, the length of the computational domain was extended so that the outlet was 20.5D away from the orifice in the downstream direction.
closest agreement with the theoretical ones.The steady and unsteady RANS results showed similar behaviors, except for the oscillation of the steady results in the vicinity of the shock waves.However, there was excellent agreement between the numerical and theoretical results, except for the region just after the shock wave, where the complexity of the interaction between the shock wave and boundary layer invalidated the 1D approximation based on the theoretical approach.These results ensured the validity of the current numerical method.

Three-Dimensional Steady Simulation
To obtain the initial values for the unsteady LES simulation, a 3D steady RANS simulation was performed using the same inlet and outlet boundary conditions as the axisymmetric simulation.The geometry and grids are shown in Figure 6.The number of grids used was 2,700,000.To obtain more reliable unsteady results by reducing the possible reflection waves from the outlet, the length of the computational domain was extended so that the outlet was 20.5D away from the orifice in the downstream direction.To verify the results, the mass flow rate predicted is shown in Figure 7.It can be seen that the mass flow rate converges at the inlet and outlet.The iso-contour distributions of the pressure, temperature, and Mach number on the central cross-sectional plane in the vicinity of the orifice are shown in Figure 8.These results are almost the same as those in the axisymmetric case.The spatial variations of these variables along the pipe axis are shown in Figure 9. Again, there is an excellent agreement between the numerical and theoretical results, with the exception of the transient region.
To verify the results, the mass flow rate predicted is shown in Figure 7.It can be seen that the mass flow rate converges at the inlet and outlet.The iso-contour distributions of the pressure, temperature, and Mach number on the central cross-sectional plane in the vicinity of the orifice are shown in Figure 8.These results are almost the same as those in the axisymmetric case.The spatial variations of these variables along the pipe axis are shown in Figure 9. Again, there is an excellent agreement between the numerical and theoretical results, with the exception of the transient region.To verify the results, the mass flow rate predicted is shown in Figure 7.It can be seen that the mass flow rate converges at the inlet and outlet.The iso-contour distributions of the pressure, temperature, and Mach number on the central cross-sectional plane in the vicinity of the orifice are shown in Figure 8.These results are almost the same as those in the axisymmetric case.The spatial variations of these variables along the pipe axis are shown in Figure 9. Again, there is an excellent agreement between the numerical and theoretical results, with the exception of the transient region.To verify the results, the mass flow rate predicted is shown in Figure 7.It can be seen that the mass flow rate converges at the inlet and outlet.The iso-contour distributions of the pressure, temperature, and Mach number on the central cross-sectional plane in the vicinity of the orifice are shown in Figure 8.These results are almost the same as those in the axisymmetric case.The spatial variations of these variables along the pipe axis are shown in Figure 9. Again, there is an excellent agreement between the numerical and theoretical results, with the exception of the transient region.As shown in Figure 10, a complicated flow structure is formed in the region downstream of the shock wave as a result of the complex interaction between the shock wave and the boundary layer at the wall surface.These complex interactions generate a strong turbulent flow, and thus cause aerodynamic noise.The following section shows how this was investigated in detail.
As shown in Figure 10, a complicated flow structure is formed in the region downstream of the shock wave as a result of the complex interaction between the shock wave and the boundary layer at the wall surface.These complex interactions generate a strong turbulent flow, and thus cause aerodynamic noise.The following section shows how this was investigated in detail.

Three-Dimensional Unsteady Simulation
This section discusses an unsteady simulation that was performed using the LES to capture the compressible acoustic field, as well as the incompressible hydrodynamic flow field.For efficient computation, the 3D RANS steady simulation result was used as an initial condition for the LES simulation.The unsteady simulation was carried out with a time step of 0.05 ms and 8,220,000 grids.In the case of the unsteady LES analysis, the unsteady flow characteristics should be captured by the numerical method combined with its mesh.The important unsteady phenomena in the CAA simulation are two-fold: one is aerodynamic noise source mechanism and the other is acoustic wave propagation.Therefore, the quality of grid for CAA simulation is generally determined by the highest frequency of the concerned frequency range of the problem.In other words, as the grid-resolution increases, the resolved frequency range also increases.An important aerodynamic noise source mechanism in the current problem is closely related to the turbulence in the boundary layer, especially, downstream of the orifice.The important length scale related to the aeroacoustic source is y+ so that the minimum grid length is determined, such that y+ = 10.The length scale related to acoustic wave propagation is the wavelength of the highest frequency component, which is set to be 10,000 Hz in the current problem.The maximum grid length is matched to be ∆ /8, where is the wavelength of the acoustic wave of 10,000 Hz.The geometry and boundary conditions were the same as those for the 3D RANS simulation, as shown in Figure 6.
Figure 11 shows the variations in the mass flow rates at the inlet and outlet in the LES simulations.Because the 3D RANS result was used for the initial condition, a transient period is seen before the converged numerical solutions occur after 1500 iterations.Figure 12 shows the distributions of the static pressure, turbulent viscosity, and vorticity on the central cross-sectional plane and in an enlarged area downstream of the orifice.The strong turbulent fluctuation due to the complex interaction between the shock and the separated boundary layer is observed in the region just downstream of the shock wave.

Three-Dimensional Unsteady Simulation
This section discusses an unsteady simulation that was performed using the LES to capture the compressible acoustic field, as well as the incompressible hydrodynamic flow field.For efficient computation, the 3D RANS steady simulation result was used as an initial condition for the LES simulation.The unsteady simulation was carried out with a time step of 0.05 ms and 8,220,000 grids.In the case of the unsteady LES analysis, the unsteady flow characteristics should be captured by the numerical method combined with its mesh.The important unsteady phenomena in the CAA simulation are two-fold: one is aerodynamic noise source mechanism and the other is acoustic wave propagation.Therefore, the quality of grid for CAA simulation is generally determined by the highest frequency of the concerned frequency range of the problem.In other words, as the grid-resolution increases, the resolved frequency range also increases.An important aerodynamic noise source mechanism in the current problem is closely related to the turbulence in the boundary layer, especially, downstream of the orifice.The important length scale related to the aeroacoustic source is y+ so that the minimum grid length is determined, such that y+ = 10.The length scale related to acoustic wave propagation is the wavelength of the highest frequency component, which is set to be 10,000 Hz in the current problem.The maximum grid length is matched to be ∆x max = λ min /8, where λ min is the wavelength of the acoustic wave of 10,000 Hz.The geometry and boundary conditions were the same as those for the 3D RANS simulation, as shown in Figure 6.
Figure 11 shows the variations in the mass flow rates at the inlet and outlet in the LES simulations.Because the 3D RANS result was used for the initial condition, a transient period is seen before the converged numerical solutions occur after 1500 iterations.Figure 12 shows the distributions of the static pressure, turbulent viscosity, and vorticity on the central cross-sectional plane and in an enlarged area downstream of the orifice.The strong turbulent fluctuation due to the complex interaction between the shock and the separated boundary layer is observed in the region just downstream of the shock wave.
There are two types of excitation sources for the vibration of a pipe wall: the compressible pressure fluctuation that is associated with acoustic waves, and the incompressible pressure fluctuation, which is associated with a turbulent flow.As shown in Figure 12c, the incompressible pressure fluctuations associated with a strong vortex convected with the mean flow can induce pipe vibration.This phenomenon is called flow-induced vibration.In addition, the acoustic waves generated by the turbulent fluctuation and its interaction with the solid surface propagate at the velocity of sound plus the mean flow velocity.The compressible pressure fluctuation associated with the acoustic waves also causes pipe vibration, which is called AIV.However, because the compressible pressure fluctuation propagates at the speed of sound, it is difficult to directly observe it from the flow field.Therefore, the compressible pressure fluctuation and incompressible pressure fluctuation were separated using a wavenumber-frequency analysis.There are two types of excitation sources for the vibration of a pipe wall: the compressible pressure fluctuation that is associated with acoustic waves, and the incompressible pressure fluctuation, which is associated with a turbulent flow.As shown in Figure 12c, the incompressible pressure fluctuations associated with a strong vortex convected with the mean flow can induce pipe vibration.This phenomenon is called flow-induced vibration.In addition, the acoustic waves generated by the turbulent fluctuation and its interaction with the solid surface propagate at the velocity of sound plus the mean flow velocity.The compressible pressure fluctuation associated with the acoustic waves also causes pipe vibration, which is called AIV.However, because the compressible pressure fluctuation propagates at the speed of sound, it is difficult to directly observe it from the flow field.Therefore, the compressible pressure fluctuation and incompressible pressure fluctuation were separated using a wavenumber-frequency analysis.

Wavenumber-Frequency Analysis
The wavenumber-frequency analysis was performed using the LES results presented in Section 4.3.Figure 13 shows the wavenumber-frequency transform of the static pressure on a cylindrical inner duct wall, i.e., r = d/2 for a length of 1 m from the pipe outlet.The frequency and time intervals are  , respectively.The average speed of sound and mean   There are two types of excitation sources for the vibration of a pipe wall: the compressible pressure fluctuation that is associated with acoustic waves, and the incompressible pressure fluctuation, which is associated with a turbulent flow.As shown in Figure 12c, the incompressible pressure fluctuations associated with a strong vortex convected with the mean flow can induce pipe vibration.This phenomenon is called flow-induced vibration.In addition, the acoustic waves generated by the turbulent fluctuation and its interaction with the solid surface propagate at the velocity of sound plus the mean flow velocity.The compressible pressure fluctuation associated with the acoustic waves also causes pipe vibration, which is called AIV.However, because the compressible pressure fluctuation propagates at the speed of sound, it is difficult to directly observe it from the flow field.Therefore, the compressible pressure fluctuation and incompressible pressure fluctuation were separated using a wavenumber-frequency analysis.

Wavenumber-Frequency Analysis
The wavenumber-frequency analysis was performed using the LES results presented in Section 4.3.Figure 13 shows the wavenumber-frequency transform of the static pressure on a cylindrical inner duct wall, i.e., r = d/2 for a length of 1 m from the pipe outlet.The frequency and time intervals are  , respectively.The average speed of sound and mean

Wavenumber-Frequency Analysis
The wavenumber-frequency analysis was performed using the LES results presented in Section 4.3.Figure 13 shows the wavenumber-frequency transform of the static pressure on a cylindrical inner duct wall, i.e., r = d/2 for a length of 1 m from the pipe outlet.The frequency and time intervals are ∆ f = 10 Hz and ∆t = 5 × 10 −5 s, respectively.The spatial and wavenumber intervals are ∆z = 3.844 × 10 −3 m, ∆θ = 10 • , and ∆k z = 1.004, respectively.The average speed of sound and mean flow Mach number over the target region are c 0 = 317 m/s and M = 0.296, respectively.The speed of sound is considered to be the smallest value in the entire flow field.
The black cut-off lines corresponding to the modes of (m, n) = (0:4, 0:3) are also drawn to make it easier to understand the behavior of the pressure waves inside the pipe in the wavenumber-frequency domains.The m-mode is the circumferential mode, and the n-mode is the radial mode.The dashed line represents the incompressible pressure wave convecting at a constant mean flow speed of 0.88U 0 .The strong compressible acoustic waves propagating at speeds of c + U 0 (downstream propagating waves) and −c + U 0 (upstream propagating waves) can be identified in the wavenumber-frequency diagram of the m = 0 and n = 0 modes.The condition for the higher-order modes is shown in Equation (10).From Equation ( 9), the variation in the z-direction wavenumber, k z , for the higher-order modes is described in Figure 13 in the parabola form, according to the n-modes [14].It can be seen from Figure 13 that the n = 0 and n = 1 modes contributes dominantly, whereas the higher-order modes with n = 2 or higher are relatively negligible.Cowling mentioned that pipe failure due to acoustic waves occurs in the range of 500-2000 Hz [8,15].In this study, the higher-order modes above n = 2 occurred at frequencies greater than 2000 Hz.Thus, they did not belong to the pipeline failure region due to acoustic waves.The line of the n = 0 mode shows the velocity of the acoustic wave propagating upstream and downstream, where the area inside the line is the compressible pressure fluctuation component, and the area outside the line is the incompressible pressure fluctuation component.Based on the results of the wavenumber-frequency analysis, it can be seen that the compressible and incompressible pressure fluctuations can be separated from the entire flow field, as shown in Figure 14.(upstream propagating waves) can be identified in the wavenumber-frequency diagram of the m = 0 and n = 0 modes.The condition for the higher-order modes is shown in Equation (10).From Equation ( 9), the variation in the z-direction wavenumber, kz, for the higher-order modes is described in Figure 13 in the parabola form, according to the n-modes [14].It can be seen from Figure 13 that the n = 0 and n = 1 modes contributes dominantly, whereas the higher-order modes with n = 2 or higher are relatively negligible.Cowling mentioned that pipe failure due to acoustic waves occurs in the range of 500-2000 Hz [8,15].In this study, the higher-order modes above n = 2 occurred at frequencies greater than 2000 Hz.Thus, they did not belong to the pipeline failure region due to acoustic waves.The line of the n = 0 mode shows the velocity of the acoustic wave propagating upstream and downstream, where the area inside the line is the compressible pressure fluctuation component, and the area outside the line is the incompressible pressure fluctuation component.Based on the results of the wavenumber-frequency analysis, it can be seen that the compressible and incompressible pressure fluctuations can be separated from the entire flow field, as shown in Figure 14.  Figure 15 compares the acoustic power level spectrum that was obtained from the inverse wavenumber-frequency transform of the compressible and total pressure fluctuation components.It shows that the acoustic power spectrum can be separated from that due to the total pressure fluctuation.Figure 16 shows the acoustic power spectrum of each circumferential mode.These were obtained from the inverse wavenumber-frequency transform of the wavenumber-frequency diagram, as shown in Figure 14.The contribution of each circumferential acoustic mode to the total acoustic power can be quantitatively assessed based on these results.This information can be utilized because the higher modes (m > 0) are known to be more critical in causing flow-induced vibration in a pipe.Carucci et al. [1,4] presented an empirical formula for predicting the acoustic power level in pipes based on flow field information.Eisinger et al., Riegel et al., and CSTI proposed standard values for the acoustic power level based on the pipe dimensions [4], and recommended designing a pipe with a value smaller than the standard value.Table 1 compares the overall acoustic power levels that are predicted using the current approach with those obtained using these Figure 15 compares the acoustic power level spectrum that was obtained from the inverse wavenumber-frequency transform of the compressible and total pressure fluctuation components.It shows that the acoustic power spectrum can be separated from that due to the total pressure fluctuation.Figure 16 shows the acoustic power spectrum of each circumferential mode.These were obtained from the inverse wavenumber-frequency transform of the wavenumber-frequency diagram, as shown in Figure 14.The contribution of each circumferential acoustic mode to the total acoustic power can be quantitatively assessed based on these results.This information can be utilized because the higher modes (m > 0) are known to be more critical in causing flow-induced vibration in a pipe.Carucci et al. [1,4] presented an empirical formula for predicting the acoustic power level in pipes based on flow field information.Eisinger et al., Riegel et al. and CSTI proposed standard values for the acoustic power level based on the pipe dimensions [4], and recommended designing a pipe with a value smaller than the standard value.Table 1 compares the overall acoustic power levels that are predicted using the current approach with those obtained using these empirical methods.Note that the empirical formulas cannot consider the influence of the mean flow, because these models were developed using a jet noise model that does not include the mean flow effects.Therefore, in this study, the acoustic power was calculated without considering the influence of the mean flow for a fair comparison.It can be seen that the empirical approaches over-estimate the acoustic power levels.The reason for this is that the empirical formula were developed using the data obtained from the real valves with more complex geometries than the present model.However, the empirical methods cannot provide detailed information about the spectral distribution of acoustic power levels that the present method can do.
Figure 15 compares the acoustic power level spectrum that was obtained from the inverse wavenumber-frequency transform of the compressible and total pressure fluctuation components.It shows that the acoustic power spectrum can be separated from that due to the total pressure fluctuation.Figure 16 shows the acoustic power spectrum of each circumferential mode.These were obtained from the inverse wavenumber-frequency transform of the wavenumber-frequency diagram, as shown in Figure 14.The contribution of each circumferential acoustic mode to the total acoustic power can be quantitatively assessed based on these results.This information can be utilized because the higher modes (m > 0) are known to be more critical in causing flow-induced vibration in a pipe.Carucci et al. [1,4] presented an empirical formula for predicting the acoustic power level in pipes based on flow field information.Eisinger et al., Riegel et al., and CSTI proposed standard values for the acoustic power level based on the pipe dimensions [4], and recommended designing a pipe with a value smaller than the standard value.Table 1 compares the overall acoustic power levels that are predicted using the current approach with those obtained using these empirical methods.Note that the empirical formulas cannot consider the influence of the mean flow, because these models were developed using a jet noise model that does not include the mean flow effects.Therefore, in this study, the acoustic power was calculated without considering the influence of the mean flow for a fair comparison.It can be seen that the empirical approaches over-estimate the acoustic power levels.The reason for this is that the empirical formula were developed using the data obtained from the real valves with more complex geometries than the present model.However, the empirical methods cannot provide detailed information about the spectral distribution of acoustic power levels that the present method can do.

Conclusions
The LES technique and a wavenumber-frequency analysis were combined to estimate the internal compressible acoustic waves and incompressible hydrodynamic pressure waves in a simple constriction-expansion pipe.First, the steady-state flow was numerically simulated, and the result was compared with a quasi-1D theoretical solution, which confirmed the validity of the current numerical method.Then, an unsteady simulation was performed to predict the unsteady pressure fluctuation inside a pipe.Next, the compressible acoustic pressure waves and incompressible

Conclusions
The LES technique and a wavenumber-frequency analysis were combined to estimate the internal compressible acoustic waves and incompressible hydrodynamic pressure waves in a simple constriction-expansion pipe.First, the steady-state flow was numerically simulated, and the result was compared with a quasi-1D theoretical solution, which confirmed the validity of the current numerical method.Then, an unsteady simulation was performed to predict the unsteady pressure fluctuation inside a pipe.Next, the compressible acoustic pressure waves and incompressible hydrodynamic pressure waves in the pipe were separated by applying a wavenumber-frequency transform to the total pressure field.Finally, the inverse wavenumber-frequency transform facilitated the estimation of the acoustic power spectrum by separating the compressible and incompressible pressure fluctuations.The estimated acoustic power level was compared with those that were obtained using the existing empirical models, which highlighted the fact that the current method can provide more detailed information about the acoustic source for the so-called AIV of a piping system, and thus facilitate systematic measures for the mitigation of this AIV.

Figure 1 .
Figure 1.Schematic diagram of typical pressure relief valve.

Figure 1 .
Figure 1.Schematic diagram of typical pressure relief valve.

Figure 2 .
Figure 2. Grids in physical domain for pipe flow for steady and unsteady computational fluid dynamics (CFD) simulations (a) pipe inlet; (b) pipe throat; (c) pipe outlet.

Figure 3 .
Figure 3.Time histories of mass flow rates.

Figure 2 .
Figure 2. Grids in physical domain for pipe flow for steady and unsteady computational fluid dynamics (CFD) simulations (a) pipe inlet; (b) pipe throat; (c) pipe outlet.

Figure 2 .
Figure 2. Grids in physical domain for pipe flow for steady and unsteady computational fluid dynamics (CFD) simulations (a) pipe inlet; (b) pipe throat; (c) pipe outlet.

Figure 3 .
Figure 3.Time histories of mass flow rates.

Figure 3 .
Figure 3.Time histories of mass flow rates.

Figure 5 .
Figure 5. Spatial distributions of (a) static pressure; (b) temperature; and, (c) Mach number along axis of pipe.

Figure 5 .
Figure 5. Spatial distributions of (a) static pressure; (b) temperature; and, (c) Mach number along axis of pipe.

Figure 5 .
Figure 5. Spatial distributions of (a) static pressure; (b) temperature; and, (c) Mach number along axis of pipe.

Figure 5 .
Figure 5. Spatial distributions of (a) static pressure; (b) temperature; and, (c) Mach number along axis of pipe.

Figure 6 .
Figure 6.Grids in physical domain for pipe flow for three-dimensional (3D) steady CFD simulations with applied boundary conditions.

Figure 7 .
Figure 7. Time histories of mass flow rates.

Figure 7 .
Figure 7. Time histories of mass flow rates.

Figure 7 .
Figure 7. Time histories of mass flow rates.

Figure 7 .
Figure 7. Time histories of mass flow rates.

Figure 10 .
Figure 10.Vorticity contours near throat: (a) separation location due to shock-boundary layer interaction; (b) shock location at center line.

Figure 10 .
Figure 10.Vorticity contours near throat: (a) separation location due to shock-boundary layer interaction; (b) shock location at center line.

Figure 11 .
Figure 11.Time histories of mass flow rates for large eddy simulation (LES).
s, respectively.The spatial and wavenumber intervals are

Figure 11 .
Figure 11.Time histories of mass flow rates for large eddy simulation (LES).

Figure 11 .
Figure 11.Time histories of mass flow rates for large eddy simulation (LES).
s, respectively.The spatial and wavenumber intervals are
radial mode.The dashed line represents the incompressible pressure wave convecting at a constant mean flow speed of 0 88 .0 U .The strong compressible acoustic waves propagating at speeds of

Figure 13 .
Figure 13.Wavenumber-frequency diagram of surface static pressure on cylindrical plane r = d/2.

Figure 14 .
Figure 14.Region of compressible pressure fluctuation component for m = 0 mode.

Figure 15 .
Figure 15.Power level spectrum of total and acoustic component.

Figure 15 .
Figure 15.Power level spectrum of total and acoustic component.

Figure 16 .
Figure 16.Power level spectrum of each mode.

Figure 16 .
Figure 16.Power level spectrum of each mode.

Table 1 .
Overall acoustic power levels found using empirical formulas.

Table 1 .
Overall acoustic power levels found using empirical formulas.