- Email: [email protected]
A frequency domain linearized Navier–Stokes method including acoustic damping by eddy viscosity using RANS
Journal of Sound and Vibration 346 (2015) 229–247
Contents lists available at ScienceDirect
Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi
A frequency domain linearized Navier–Stokes method including acoustic damping by eddy viscosity using RANS Andreas Holmberg a,b,n, Axel Kierkegaard b, Chenyang Weng b,c a b c
KTH, CCGEx – the Competence Center for Gas Excange, Teknikringen 8, 100 44 Stockholm, Sweden KTH, Marcus Wallenberg Laboratory for Sound and Vibration Research, Teknikringen 8, 100 44 Stockholm, Sweden KTH, Linné Flow Centre, Teknikringen 8, 100 44 Stockholm, Sweden
a r t i c l e i n f o
abstract
Article history: Received 30 July 2013 Received in revised form 6 February 2015 Accepted 16 February 2015 Handling Editor: P. Joseph Available online 11 March 2015
In this paper, a method for including damping of acoustic energy in regions of strong turbulence is derived for a linearized Navier–Stokes method in the frequency domain. The proposed method is validated and analyzed in 2D only, although the formulation is fully presented in 3D. The result is applied in a study of the linear interaction between the acoustic and the hydrodynamic field in a 2D T-junction, subject to grazing flow at Mach 0.1. Part of the acoustic energy at the upstream edge of the junction is shed as harmonically oscillating disturbances, which are conveyed across the shear layer over the junction, where they interact with the acoustic field. As the acoustic waves travel in regions of strong shear, there is a need to include the interaction between the background turbulence and the acoustic field. For this purpose, the oscillation of the background turbulence Reynold's stress, due to the acoustic field, is modeled using an eddy Newtonian model assumption. The time averaged flow is first solved for using RANS along with a k-ε turbulence model. The spatially varying turbulent eddy viscosity is then added to the spatially invariant kinematic viscosity in the acoustic set of equations. The response of the 2D T-junction to an incident acoustic field is analyzed via a plane wave scattering matrix model, and the result is compared to experimental data for a T-junction of rectangular ducts. A strong improvement in the agreement between calculation and experimental data is found when the modification proposed in this paper is implemented. Discrepancies remaining are likely due to inaccuracies in the selected turbulence model, which is known to produce large errors e.g. for flows with significant rotation, which the grazing flow across the T-junction certainly is. A natural next step is therefore to test the proposed methodology together with more sophisticated turbulence models. & 2015 Elsevier Ltd. All rights reserved.
1. Introduction Duct junctions and sudden area changes are common elements in pipe networks. While sometimes added to scatter sound, e.g. expansion chambers in exhaust system configurations, especially T-junctions and wall corrugations are known to cause severe resonances, see e.g. [1–3]. The governing aeroacoustic phenomena are generally well known, but the prediction of frequency and magnitude of the acoustic field in the design stage is often beyond reach. Unfortunately, building and
n
Corresponding author at: KTH, CCGEx – the Competence Center for Gas Excange, Teknikringen 8, 100 44 Stockholm, Sweden. E-mail address: [email protected] (A. Holmberg).
http://dx.doi.org/10.1016/j.jsv.2015.02.030 0022-460X/& 2015 Elsevier Ltd. All rights reserved.
230
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
testing prototypes is a time-consuming and expensive process. There is therefore a demand for more efficient and accurate prediction tools. In this paper, we describe a numerically efficient method for calculating linear acoustic propagation, scattering and hydrodynamic interaction. The numerical method is a modification of a linearized Navier–Stokes equations method for the frequency domain, described by Kierkegaard et al. [4]. It is a hybrid approach, in which the flow is divided into a mean and a time dependent but coherent perturbation part. As the mean flow is calculated by RANS using a k–ε turbulence model and the perturbation equation set are linear, the method is very time efficient. The modification we propose begins with substituting the flow decomposition used in [4] by a triple decomposition, where the flow is divided into a time average, a phase average and a random part. This choice stems from the need to include the oscillation of the turbulence Reynold's stress, due to the acoustic field. This oscillation induces an additional acoustic dissipation process, which in general cannot be neglected [5]. The method of using a triple decomposition in order to include the turbulence damping in the linearized Navier–Stokes equations (LNSE) was first introduced in the Ph.D. thesis by Holmberg [6], which is the basis of this paper. Recently, Gikadi et al. [7] presented a paper similar to this, in which 3D LNSE calculations were performed, also using the triple decomposition to include the turbulence effect on the acoustic field. The background flow was in their work calculated using LES, while here RANS is used. We apply our method in the study of a T-junction of rectangular ducts, subject to grazing flow of Mach 0.1 across the opening, for which experimental results exist. The interaction between the acoustic field, the geometry and the shear layer across the junction is studied via a scattering matrix. From this matrix, the frequency ranges of amplification and attenuation can be predicted [8], as well as the frequencies of whistling [9]. Previously, the scattering matrix for several grazing flow Mach numbers was determined by Karlsson and Åbom [8] for cylindrical ducts. A numerical solution has also been obtained by Föller et al. [10], who superimposed plane waves on a large eddy simulation, for a grazing flow Mach number of 0.1. The only significant deviation between experimental and simulation results, was found in the scattering matrix element describing acoustic transmission from the upstream to the downstream main pipe. Föller et al. [10] argued that the much stronger transmission found in the experiments was unphysical. However, this large transmission is also found in a new set of experiments conducted for the rectangular duct junction studied here. We will show the Mach 0.1 grazing flow case, and compare to our numerical results. An interesting feature of the paper by Gikadi et al. [7] is that they also studied a T-junction – in their case a junction of circular pipes, hence the need for 3D simulations. The simulation results in [7] are compared to the experimental results of Karlsson and Åbom [8] – just as in this paper – for the case of Mach 0.1 grazing flow. In the conclusion section, their results will be discussed in light of the results presented in this paper. This paper is outlined as follows; in Section 2, the interaction between the hydrodynamic and acoustic fields in the Tjunction is described. The modification of the governing equations is derived in Section 3. The methods to analyze the calculated results are described in Section 4, while the numerical implementation of the equations and the analysis methods are discussed in Section 5. Results are shown and discussed in Section 6, while a summary of the paper and an outlook is presented in Section 7. 2. Theory of the T-junction When the acoustic particle velocity turns into the sidebranch at the upstream edge of the T-junction, part of the acoustic energy is shed as a small disturbance [11]. In the case of a mean flow in the system, a shear layer is formed across the Tjunction opening, and the disturbance is conveyed by the mean flow across this layer. Due to the instability of the shear layer, the disturbance grows as it travels across the junction [12]. The interaction between this disturbance, the mean flow and the acoustic field can be studied e.g. by the vortex-sound theory by Howe [13]. He expresses the transfer of power from the hydrodynamic field to the acoustic field via the integral Z Π ¼ ρ0 ω v u dV (1) where V is the volume of integration, ω is the distribution of vorticity and v and u are the incompressible and the acoustic velocity fields. As the integral is zero when u and v are parallel, the case of plane wave propagation in a straight duct yields no energy transfer. In the T-junction however, the angle between the particle velocity and the convection velocity is nonzero, and thus, energy transfer can occur, where the direction of transfer depends on the phase of the acoustic field. Since the disturbances are shed when the acoustic particle velocity turns into the sidebranch, the integral is positive during the first half of the acoustic period, thus resulting in a negative sound power contribution. During the second period, the acoustic particle velocity changes sign, resulting in a positive sound power contribution. A net increase is possible since the magnitude of the disturbance, shed by the acoustic wave grows as it is conveyed across the opening. Depending on the ratio between the acoustic and the hydrodynamic timescales, i.e. on the Storuhal number, there will be a net increase or decrease of the acoustic energy. For low acoustic velocity amplitudes, i.e. of the order of O( 3) and below when normalized by the mean flow velocity [14], the magnitude of the shed disturbance is invariant on (although correlated with) the incoming acoustic amplitude, which, according to Eq. (1) yields a linear relation between the incoming and amplified acoustic fields [14]. The rate of growth of the shed disturbance depends on the shear layer as well as the hydrodynamic wavelength of the disturbance. Michalke [12] investigated this for free inviscid parallel flows having a hyperbolic-tangent profile, using modal
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
231
stability analysis. His results are unfortunately not directly applicable to the T-junction, especially not for mixed grazing-bias flows, that is, flow in all three branches. Previous studies also clearly show that the power amplification depends strongly on e.g. the edges of the junction [15], the amount of grazing and bias flow [16,17] and the cross sections of the ducts [17]. It is not surprising therefore, that no reasonable versatile engineering models of the T-junction exist. 3. Derivation of the governing equations In this section, we present the derivation of a modified version of the frequency domain formulation of the linearized Navier–Stokes equations presented by Kierkegaard et al. [4]. As mentioned in the introduction, we propose this modification in order to include dissipation of the acoustic field in the presence of strong vorticity. As in [4], the starting point is the Navier–Stokes equations [18], here written as
ρ where
Dρ ∂u þρ i ¼ 0 Dt ∂xi
(2)
Dui ∂p ∂τij ¼ þ þ ρF i ∂xi ∂xj Dt
(3)
∂ui ∂uj 2 ∂uk þ δij ∂xj ∂xi 3 ∂xk
τij ¼ μ
(4)
Here, ρ, u and p, are the density, velocity and pressure fields. F is a volume force, μ the dynamic viscosity and δ the Kronecker delta function. By assuming an isentropic process, we can relate the density to the pressure, and decouple the equation system from the energy equation. Now, instead of dividing the flow field into a mean and a coherent perturbed field, as done in [4], we adopt the triple decomposition used by Reynolds and Hussain [5]. The velocity and pressure fields are here divided into three parts e i ðx; t Þ þ u0i ðx; t Þ ui ðx; t Þ ¼ ui ðxÞ þ u
(5)
0
e ðx; t Þ þ pðx; t Þ pðx; t Þ ¼ pðxÞ þ p
(6)
The last term vanishes when a phase average is calculated, and hence it is the turbulent part of the flow. The second term vanishes when a time average is calculated, and is therefore the coherent perturbed field, which includes the acoustic waves. The first term is simply the time averaged flow field. For the density we assume a simpler decomposition, namely
ρðx; t Þ ¼ ρ þ ρe ðx; t Þ
(7)
The turbulent part of the density is thus assumed negligible, and the time averaged density is assumed constant throughout the flow. In other words, the flow is approximated as a harmonic field superposed on to an incompressible flow field. This is likely to be accurate for low Mach numbers (below Mach 0.3 say). The mean static pressure is allowed to vary in
Fig. 1. A T-junction; definitions of coordinate systems, positive directions for the Mach numbers and numbering of side branches used in the paper.
Table 1 Physical size of the computational domain.
Duct length [mm] RANS Duct length [mm] LNSE Duct height/width [mm] Wave decomposition zone, distance from junction [mm] Buffer zone length [mm]
Branch I
Branch II
Branch III
1500 1500 25 300—1200 500
1500 1500 25 400–1200 500
1000 1500 25 400–1400 500
232
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
space due to velocity gradients, while the mean density is assumed constant throughout the flow. Now, applying Eqs. (5) and (7) to the continuity equation, we get e i þ u0i e ui þ u e ∂ ρþρ dρ þ ¼0 (8) ∂xi dt The solution we seek is the coherent perturbed field. Therefore we take the phase average of Eq. (8), linearize the result and introduce the assumption of monochromatic waves with ending up with e iωt time dependence. The result, transformed to the frequency domain, is e e ∂ρ ∂ui ∂u e ¼ ρ i ui þ iω ρ (9) ∂xi ∂xi ∂xi Before continuing with the momentum equation, we note that the turbulent part of Eq. (8) can be written as
ρ þ ρe
∂u0i e ∂ρ ¼ u0i ∂xi ∂xi
(10)
We now apply the flow decomposition, i.e. Eqs. (5)–(7) on the momentum equation, i.e. Eq. (3), and rearrange the result using Eq. (10). The result is ∂u þ u e i þu0i ei d u ∂u0i ∂u0i u0j e ∂ρ i e j þ u0j ej ρ þ ρ e e þ uj þ u ρ þ ρe þ ρþρ þu0i u0j þ uj þ u dt ∂xj ∂xj ∂xj ∂xj 0 1 0 e j þ u0j ∂ uj þ u e k þ u0k e i þ u0i e þp ∂ ρþρ ∂ @∂ u i þ u 2 ∂ uk þ u ¼ þ μ þ δij A þ ρ þ ρe F i (11) ∂xj 3 ∂xi ∂xj ∂xk ∂xi Now, nothing has been said about the volume force Fi yet. In this work we will implement it as a harmonically fluctuating e Fe i is zero, and we therefore volume force multiplied by the mean density. Hence, in the computations, the nonlinear term ρ e F i ¼ ρFe i . After introducing this implementation, we continue by subtracting the time average from the phase have ρ þ ρ
Fig. 2. Time average flow field: (a) velocity magnitude [m/s] and (b) eddy viscosity (νT) [m2/s].
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
233
average of Eq. (11). Linearizing the result thus yields e ∂e uj ∂e ρ uk ∂e ui ∂ 2 ∂e ρ ddteu i þ ue j ∂∂xeuji þ uj ∂∂xujj þ ρe uj ∂∂xeuji ¼ ∂x þ μ þ δ ∂xj 3 ∂xk ij ∂xj ∂xi i þ ρFei þ ρ
∂ u0i u0j 〈u0i u0j 〉 ∂xj
ρe
∂ 〈u0i u0j 〉
ρe
u0i u0j ∂ ∂xj ∂xj
! e ρ
∂〈u0i u0j 〉 ∂xj
(12)
We now have a closure problem due to the unknown terms ∂〈u0i u0j 〉 u0i u0j ∂e r ij x¼ ∂xj ∂xj
(13)
e e ρ ∂ρ ρe e ∂〈u0i u0j 〉 e ρ (14) r ij ¼ 〈u0i u0j 〉 þ u0i u0j ∂ ∂xj ∂xj ∂xj e ρ r ij . However, while this is true for small which must be modeled. To start with, we consider e r ij to be negligible compared to e e ρ e i =∂xj Þ þ perturbations, it should be kept in mind that e r ij is not necessarily smaller than the terms ð∂=∂xj Þμ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij Þ, which we have decided to keep. The reason why we keep the latter terms is that, depending on ð∂u how e r ij is modeled, they may not be negligible everywhere in the computational domain. Before continuing, it is interesting to note the physical interpretation of e r ij . Since e r ij is expressed as the difference between the time average and the phase average of the turbulent Reynold's stresses, it can be interpreted as the oscillation of the turbulent Reynolds's stress due to the acoustic field, as is also discussed in [5]. Now, in order to model e r ij , we start with the time average term. According to the Boussinesq hypothesis [19], the deviatoric Reynolds stress is proportional to the mean rate of strain 2 ∂ui ∂uj þ (15) ¼ 2νT Sij u0i u0j þ kδij ¼ νT 3 ∂xj ∂xi where νT is the eddy viscosity, and k is the turbulent kinetic energy. Now, for the phase average term, the idea of a Newtonian eddy model introduced in Reynolds and Hussain [5] could be used
Fig. 3. Perturbation density for sidebranch excitation: (a) 690 Hz and (b) 1370 Hz.
234
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
ei ej ek ∂ ui þ u ∂ uj þ u 2 ∂ uk þ u 2 þ δij þ kδij 3 3 ∂xj ∂xi ∂xk e e e ∂u ∂ u ∂ui ∂u i 2 ∂u k 2 j j ¼ νT þ þ δ þ þ kδij 3 ∂xj ∂xi ∂xj ∂xi 3 ∂xk ij 2 ¼ 2νT Sij þ e S ij þ kδij 3
〈u0i u0j 〉 ¼ νT
(16)
where the divergence of the time average velocity is set to zero, since it become zero if the time average of the continuity equation is linearized. Now, the above equation may be considered as improper, since the Boussinesq closure is an equilibrium expression that assumes the anisotropy tensor to depend only on the local instantaneous value of the mean strain rate tensor. The acoustic shear on the other hand, is a time dependent quantity, and thus the time response (phase) between the instantaneous value of the mean shear and of the acoustic shear should be taken into account. This is however not done here. Instead we assume that the effect of this is negligible. Thus, we get e j 2 ∂u e i ∂u ek ∂u e þ δij (17) r ij ¼ νT ∂xj ∂xi 3 ∂xk The eddy viscosity is readily obtained from the RANS calculations. Now remains the coupling between the pressure and the density. As previously mentioned, we here assume an isentropic process. The equation of state can thus be approximated by [21] Dp Dρ ¼ c2 Dt Dt
(18)
In all published work concerning the frequency domain linearized Navier–Stokes equations method [4, 22,23], the quiescent medium approximation of Eq. (18) is used, seemingly without any noticeable errors induced. We here chose to continue in the line of the previous work, due to the increased speed of the computations this approximation benefits from. Thus we write
ρe ¼ c2 ρe
e (x-component) for sidebranch excitation: (a) 690 Hz and (b) 1370 Hz. Fig. 4. Perturbation velocity u
(19)
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
The governing equations are finally e j 2 ∂u ∂u ei e e i ∂u ek e ∂u du ∂u ∂u ∂ρ ∂ e j i þ uj j þ ρ e uj i ¼ c2 þ þu ρ μ þ ρνT þ δij þ ρFei dt ∂xj ∂xj ∂xj ∂xi ∂xj ∂xj ∂xi 3 ∂xk ui
e e ∂ρ ∂ui ∂u e ¼ ρ i þ iω ρ ∂xi ∂xi ∂xi
235
(20)
(21)
Comparing these equations with those derived by Kierkegaard et al. [4], where the flow decomposition consisted of a division into a time average part and a coherent perturbed part, we can note that the only difference is the introduced eddy viscosity, which is added to the kinematic viscosity. 4. Analysis method In the case of a duct system where the gas flow only have linear effects on the acoustic field, N-port networks [20] can be used to calculate plane wave propagation and scattering. Each element in the network is then treated as an N-port – a linear time invariant system relating an input to an output. For T-junctions, orifice plates, shallow cavities and other elements in which strong shear layers are present in the regions of wave propagation, nonlinear coupling between the acoustic and the hydrodynamic fields may occur. Using an N-port network to study the installation effects of these elements is thus a questionable approach. However, recently N-ports determined in the low level – and thus linear – regime have been used to predict whistling frequencies of T-junctions [9] and orifice plates [23]. For planar waves, the N-port applied on the Tjunction is referred to as the 3-port scattering matrix, which is defined accordingly [8] 0 I 1 2 30 I 1 0 I 1s pþ pþ ρI τI;II τI;III p B II C 6 7B pII C B II C s B p þ C ¼ 4 τII;I B C ρ τ (22) II;III 5@ A þ II @ A @ p þ A or p þ ¼ Sp þ p þ pIII þ
τIII;I τIII;II
ρIII
pIII
pIII þ
e (y-component) for sidebranch excitation: (a) 690 Hz and (b) 1370 Hz. Fig. 5. Perturbation velocity v
236
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
where S is the scattering matrix, p- is the incoming acoustic waves, p þ is the outgoing waves, and psþ the sound that is generated independent on the incoming sound waves, see Fig. 1. The diagonal elements of the scattering matrix thus represent the reflection coefficients towards the T-junction, while the other elements correspond to transmission coefficients between the three branches. Another intrinsic property of the T-junction that can be studied using via the scattering matrix is the power balance between outgoing and incoming energy [8]. The scattering matrix can also be used to determine if a (linear) reflection connected to the T-junction is sufficient to generate a nonlinear response, resulting in a whistle, as was experimentally shown by Karlsson and Åbom [9]. Kierkegaard et al. [23] attempted to predict whistling using numerically obtained scattering matrices of orifice plates subject to bias flow. They found out however, that the predictions where quite sensitive to errors (3 percent) in the scattering matrices. Therefore a sensitivity analysis of the stability analysis should always be made. Due to its many uses it is clear that the scattering matrix is a powerful tool if calculated correctly. We will consider the error in the scattering matrix when the success of the suggested turbulence damping is discussed.
5. Implementation The calculation procedure consists of three main steps; first the time averaged part of the flow is solved for, using an incompressible, steady-state RANS method. Then the coherent perturbation equations are solved, and finally, the propagating waves are extracted and used to calculate the scattering matrix. In this work, we used Comsol 4.2 with a MatLab link. Due to computational cost, the computational realm is set to 2D, both for the RANS solution and the acoustic solution. The physical size of the T-junction is shown in Table 1. The meshes used in this work are described in detail in Appendix A.
5.1. RANS implementation The mean flow is calculated using incompressible RANS equations with k–ε turbulence modeling [19]. The reason for using the k–ε turbulence model is that it is one of the most common and simplest turbulence models found in commercial software. The turbulent viscosity is obtained from the turbulent kinetic energy k and the dissipation rate of turbulence 2 energy ε, using the model constant C μ ¼ 0:09, via the equation νT ¼ C μ ðk =εÞ.
Fig. 6. Magnitude of the scattering matrix elements at zero bias flow; black: experiments, red: calculations e r ij ¼ 0, blue: calculations e e i =∂xj Þþ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij . (For interpretation of the references to color in this figure legend, the reader is referred to the web r ij ¼ νT ð∂u version of this article.)
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
237
The time averaged solution used as input to the acoustic calculations in the next section is obtained without using wall functions. To achieve convergence however, wall functions [24] are used in two consecutive pre-steps; first, starting from a uniform inflow velocity profile, the flow is calculated in a duct of 30 m in order to obtain the outlet flow velocity, pressure, turbulence kinetic energy k, and the dissipation rate of turbulence energy ε. The outlet properties are then used as inlet condition in the main branch I, see Fig. 1. The second step involving wall-functions is a T-junction mean flow calculation, where the objective is to generate an initial condition of the flow, used to setup calculations without wall functions.
Fig. 7. Phase of the scattering matrix elements at zero bias flow; black: experiments, red: calculations e r ij ¼ 0, blue: calculations e i =∂xj Þ þ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij . (For interpretation of the references to color in this figure legend, the reader is referred to the web e r ij ¼ νT ð∂u version of this article.)
Fig. 8. Perturbation density for sidebranch excitation 1370 Hz (e r ij ¼ 0Þ.
238
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
5.2. NSE implementation The sound fields are excited by applying the source terms F in Eq. (28) to one duct branch at a time. This is done in order to determine the scattering matrix, which can only be done if there are at least three linearly independent sound fields. Hence, the perturbed field is solved for three times. The maximum for the continuous function for F is set 1500 mm from the junction in each branch, and the function is continuous and non-zero in a region of 0.05 m in both directions. The acoustic boundary layers are assumed to have negligible effects on the end result. Hence, we implement a slip condition for the perturbation velocity on all walls. In theory, the results obtained using the methodology suggested here, will not be influenced by the end reflections in the duct. However, for numerical stability reasons, the reflections should be kept to a minimum. We therefore apply non-reflecting boundary conditions at in- and outflow boundaries, which based upon the Helmholtz equations are formulated as e ¼ ikρ e; n^ U∇ρ
e ¼ iku e; n^ U ∇u
e ¼ ikv e n^ U∇v
(23a,b,c)
e is zero for vertical boundaries and u e is zero for horizontal boundaries. However, these boundary conditions break where v down in the presence of vorticity. For this reason, buffer zones with artificial viscosity are placed before each duct end. The mean flow is not calculated for the buffer zones instead the values of the mean flow quantities at the duct ends are extended to hold for the entire buffer zone as well. This is done in order to avoid acoustic scattering at the interface.
5.3. Wave decomposition To obtain the scattering matrix, the upstream and downstream propagating waves must be extracted from the solved perturbation fields. We here choose the acoustic density instead of acoustic pressure, but the choice is arbitrary, as the constant factor between the two quantities imposed by the isentropic relation disappears when calculating the scattering matrix. As we have implemented a wall slip condition for the perturbation velocity, there is no acoustic boundary layer, and assuming the acoustic damping due to bulk viscosity to be negligible we define the wavenumbers as being real-valued. We can then setup a nonlinear equation system with the real valued phase, wavenumber and amplitude of each wave as unknowns
ρe xj ¼ ρe þ eiðφ þ þ k þ xj Þ þ ρe eiðφ k xj Þ
(24)
The coordinate axis x is positive outwards and negative into the T-junction in each branch, and the scattering matrix is determined at the position xj ¼0. To solve the equation we apply a nonlinear curve-fitting algorithm, were for each xj, the input is the cross section mean value of the acoustic density. To suppress errors due to vorticity, the cross section mean is calculated at 30 points for each cross section, while the number of cross sections j is set to 150 per duct, in the wave decomposition zones shown in Table 1. In order to obtain a wave decomposition zone of 1000 mm for the sidebranch, it is made 500 longer for the acoustic domain than it is for the RANS calculations. The mean flow quantities taken at the cross section a distance of 1000 mm away from the junction, where the RANS domain ends, extended to hold for the extra duct length in the acoustic domain, just as was done for the buffer zones.
Fig. 9. Perturbation density for sidebranch excitation 1370 Hz (e r ij ¼ 0Þ, zoomed out view.
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
239
6. Results and discussion As already mentioned, the mean flow configuration we consider in this paper is that of a pure grazing flow of Mach 0.1, implying a Reynolds number of approximately 47600 based upon the duct height. To study the effects of the proposed modification of the numerical method, results are shown both with and without the modification, and are compared to experimental results. We also include a grid convergence study for both calculation methods and for the RANS calculations in Appendix A. The scattering matrix data is compared to the experimental results found in [17], for a T-junction consisting of rectangular ducts with a cross section of 120 25 mm2, subject to the same flow case as is studied here. While not being true 2D, the rectangular junction will to a large extent behave as a 2D edge in terms of coherent vortex shedding and acoustic propagation times. Concerning the cross-sectional dimensions of the experimental rig, it is argued in the experimental paper [17] that based on the findings by Alfredsson and Johansson [25], the width to height ratio is large enough for the rig to exhibit fully developed turbulent channel flows upstream of the T-junction. Another support for comparing 2D calculations with the rectangular rig experiments is found in the paper by Salt et al. [26], who performed PIV-measurements, acoustic measurements and FEM modelling on a T-junction with a width/height ratio of 117/70 in the sidebranch and 127/70 in the main branches. They conclude that neglecting the fluctuating component in the third dimension only have minor effects on the flow-acoustic coupling mechanism. As the width/height ratio is smaller in the paper by Salt et al. [26] than for the case compared to here [17], the experimental results are expected to essentially behave as if the rig was 2D. The magnitude of the velocity field and the turbulent viscosity, obtained from the final RANS calculation are shown in Fig. 2, where the distinct shear layer across the opening can be observed. There is also a clear recirculation in the sidebranch. It is interesting as well to note that even in the center of the duct, the turbulent viscosity is two orders of magnitude larger than the kinematic viscosity. A few samples of typical coherent perturbation fields are given in Figs. 3–5, where the perturbation density and velocity components are shown at 690 Hz (a Strouhal number of 0.50) and at 1370 Hz (a Strouhal number of 1.0), for sidebranch excitation (branch III). In these figures, we have used the closure assumption e e i =∂xj Þ þ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij . The resulting scattering matrices from this assumption, as well as data r ij ¼ νT ð∂u obtained using e r ij ¼ 0, are shown together with experimental data in Figs. 6–7. When approaching the Strouhal number of unity, a strong peak can be observed in the magnitude of the calculated e i =∂xj Þ þ ð∂u e j =∂xi Þ scattering matrix obtained setting e r ij ¼ 0. This peak disappear using the model e r ij ¼ νT ð∂u e k =∂xk Þ δij Þ. We can also note that at Strouhal numbers less than 0.5 say, the results using ð2=3Þð∂u e e i =∂xj Þ þ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij is better than when setting e r ij ¼ νT ð∂u r ij ¼ 0, while there is a midrange of Strouhal numbers where setting e r ij ¼ 0 results in almost perfect agreement (of the magnitude) to the experimental data. In order to further study the peak at the Strouhal number of approximately 1.0, we present two figures capturing the density perturbations at 1370 Hz for sidebranch excitation and the case where e r ij ¼ 0, see Figs. 8 and 9. In these figures strong instabilities are clearly observed which are convected with the time averaged flow, while at the same time forming a wavelike pattern. Their presence stems from an unbalance in perturbation energy, due to the lack of a proper dissipation term in the perturbation equations. Although this effect could probably be treated to some extent by artificially increasing the kinematic viscosity, it is more physically relevant to include the effect of the Reynolds stresses via the eddy viscosity. Mathematically the difference is that the eddy viscosity is spatially varying, which implies that the turbulence damping can be zero also when the shear of the perturbed quantities is not, and vice versa. 7. Conclusion Propagation, scattering and hydrodynamic interaction of acoustic waves have been calculated using a linearized Navier– Stokes method in the frequency domain. The decomposition of the flow into time average, phase average and random parts yields an additional term in the linearized equation set. This term was modeled by the Newtonian eddy model [5], i.e. proportional to the turbulent viscosity. The importance of retaining and modeling this term is studied for the case of a rectangular T-junction subject to grazing flow of Mach 0.1. It should be noted that the validation and analysis of the proposed method is in this paper limited in 2D, although the formulation is fully presented in 3D. Comparing to experimental data, the prediction of the acoustic-hydrodynamic interaction is several factors stronger, unless the modification of the equations suggested in this paper is implemented. As the turbulent viscosity is already solved for during the RANS calculations and simply added to the kinematic viscosity, there is no additional computation time involved. This modified version of the method presented by Kierkegaard et al. [4] is therefore very promising. A natural next step is to try other viscosity models, since the k-ε model is known to be inaccurate for e.g. flows with normal-stress anisotropy, high strain rates and curvatures [27]. Comparing with the results presented by Gikadi et al. [7], it is clear that their results support the observations made here; i.e. that the inclusion of the turbulence damping on the perturbation quantities yields stable and physical relevant results, while neglecting the turbulence effect may cause severe over predictions (mostly) of the acoustic amplification by the hydrodynamic shear flow. The results by Gikadi et al. [7] also shows frequencies where a small error is induced by including the turbulence damping, although at a very few number of frequency points compared to what is found here. It is not surprising that their results are more accurate however, since they solve the background flow using an LES in 3D, while we use 2D RANS with k-ε. The method to choose depends on the accuracy required. We note though that the method
240
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
presented here is easily extended in 3D, which would be faster than the LES method presented by Gikadi et al. [7]. It would therefore be a very interesting case to compute, preferably with several different turbulence models.
Appendix A In this appendix, the main results of conducted grid convergence studies are presented. For the RANS calculations, three different meshes were tested, see Table A1. In all three meshes, quadrilateral elements were used below yoh/25, where y is the wall normal direction and h¼0.025 m is the duct height. The height of the wall-adjacent element was kept constant as the mesh was refined. Everywhere else in the computational domain, triangular elements were used. The convergence was studied by plotting the x-component of the time averaged velocity at three different lines, all with the length of the duct height, see Fig. A1. The velocity profiles at the three lines are shown in Figs. A2–A4. In Fig. A5, the magnitude of the calculated scattering matrices are shown for each RANS mesh, where the acoustic closure e i =∂xj Þ þð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij is used. The finest acoustic mesh, see Table A2, was applied in hypothesis e r ij ¼ νT ð∂u this study. The velocity profiles obtained using the medium and fine RANS meshes are almost identical, while the course RANS mesh yields slightly different profiles. The influence of this on the acoustic result is in Fig. A5 shown to be tremendous. There is also a slight variation between the results obtained by the medium and the fine RANS mesh, which is seen in the highest and the lowest Strouhal number ranges. The fine RANS mesh solution was applied in the end, and was also used in the acoustic mesh convergence study. A closeup of the medium acoustic mesh, which was used for the final acoustic results in this paper, is shown in Fig. A6. The mesh consists of three main parts. At the downstream edge, a quadrilateral mesh was applied in a rectangle up to a wall normal direction of yoh/20, and with a streamwise length of h/2. Far from the T-junction, in the plane wave propagation regions of the ducts, a quadrilateral mesh was also applied. In between these two regions, a triangular mesh was used. The number of triangular elements across the duct opening, the number of wall normal quadrilateral elements, the maximum growth rate and the maximum element size was specified for three different meshes, see Table A2. In the plane wave propagation regions (far field), the wall normal size ratio between the largest element and the smallest, wall adjacent element was set to 3. A variant of the medium mesh was also tested, by setting this size ratio to 1, see Fig. A7. The result of the convergence study (for the size ratio set to 3) is shown Figs. A8 and A9, in terms of the magnitude of the scattering matrix elements at five frequencies, viz. 100, 400, 700, 1000 and 1400 Hz. The sensitivity to the wall normal size ratio is investigated by computing the scattering matrix in the full frequency range of interested. The result of this study is shown in Fig. A10 for the case without turbulent damping, and in Fig. A11 for the case with this damping implemented.
Table A1 RANS mesh properties. Mesh
Coarse
Medium
Fine
Triangular elements Quadrilateral elements Total number of elements Max. element size Min. element size Max. element growth Number of elements at yo h/25
445,676 329,840 775,516 h/25 h/2500 1.025 40
258,344 164,920 423,264 h/25 h/2500 1.05 20
156,301 491,60 205,461 h/15 h/2500 1.05 10
Fig. A1. Figure showing the lines at which the x-component of the time average velocity is plotted. .
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
241
Fig. A2. Time averaged velocity component in the x-direction at cut line A. Blue dashed line: coarse RANS mesh, red line: medium RANS mesh, green line: fine RANS mesh. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
Fig. A3. Time averaged velocity component in the x-direction at cut line B. Blue dashed line: coarse RANS mesh, red line: medium RANS mesh, green line: fine RANS mesh. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
Fig. A4. Time averaged velocity component in the x-direction at cut line C. Blue dashed line: coarse RANS mesh, red line: medium RANS mesh, Green line: fine RANS mesh. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
242
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
e i =∂xj Þþ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij , black: coarse RANS mesh, red: Fig. A5. Magnitude of calculated scattering matrix elements, using e r ij ¼ νT ð∂u medium RANS mesh, blue: fine RANS mesh. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
Table A2 Acoustic mesh properties. Mesh
Coarse
Medium
Fine
Triangular elements Quadrilateral elements Max. element size Max. growth rate Elements across opening No. of wall normal elements at downstream edge No. of wall normal elements in plane wave region Approximate computation time per frequency point on standard laptop PC
2042 6005 h/5 1.2 25 3 5 4
7668 23,900 h/10 1.1 50 5 10 5
29,689 95,600 h/20 1.05 100 10 20 45
Fig. A6. Close up of the medium mesh used for acoustic calculations. .
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
243
Fig. A7. Close up of the medium mesh without wall normal grid stretching in the far field. .
e i =∂xj Þ þ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij , black: coarse acoustic mesh, red: Fig. A8. Magnitude of calculated scattering matrix elements, using e r ij ¼ νT ð∂u medium acoustic mesh, blue: fine acoustic mesh. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
244
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
Fig. A9. Magnitude of calculated scattering matrix elements, using e r ij ¼ 0, black: coarse acoustic mesh, red: medium acoustic mesh, blue: fine acoustic mesh. For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
245
Fig. A10. Magnitude of calculated scattering matrix elements, using e r ij ¼ 0, black: experimental results, red: no wall normal grid stretching in the far field, see Fig. A7, blue: wall normal grid stretching in the far field applied, see Fig. A6. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
246
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
e i =∂xj Þ þ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij , black: experimental results, red: no Fig. A11. Magnitude of calculated scattering matrix elements, using e r ij ¼ νT ð∂u wall normal grid stretching in the far field, see Fig. A7, blue: wall normal grid stretching in the far field applied, see Fig. A6. Note that the blue curve lies almost exactly atop the red curve. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
References [1] S. Ziada, Flow-excited acoustic resonance in industry, Journal of Pressure Vessel Technology 132 (2010) 1–9. [2] P. Lafon, S. Caillaud, J.P. Devos, C. Lambert, Aeroacoustical coupling in a ducted shallow cavity and fluid/structure effects on a steam line, Journal of Fluids and Structures 18 (2003) 695–713. [3] S.P.C. Belfroid, D.P. Shatto, R.M.C.A.M. Peters, Flow induced pulsations caused by corrugated pipes, Proceedings of ASME Pressure Vessel and Piping Division Conference, 2007, San Antonio, Texas. [4] A. Kierkegaard, S. Boij, G. Efraimsson, A frequency domain linearized Navier–Stokes equations approach to acoustic propagation in flow ducts with sharp edges, Journal of Acoustic Society of America 127 (2) (2010) 710–719. [5] W.C. Reynolds, A.K.M.F. Hussain, The mechanics of an organized wave in turbulent shear flow. Part 3. Theoretical models and comparisons with experiments, Journal of Fluid Mechanics 54 (2) (1972) 263–288. [6] Holmberg, A., Aeroacoustic Characterization using Multiport Methods, PhD Thesis, Royal Institute of Technology, Stockholm, 2012. [7] J. Gikadi, S. Föller, T. Sattelmayer, Impact of turbulence on the prediction of linear aeroacoustic interactions: acoustic response of a turbulent shear layer, Journal of Sound and Vibration 333 (24) (2014) 6548–6559, http://dx.doi.org/10.1016/j.jsv.2014.06.033. [8] M. Karlsson, M. Åbom, Aeroacoustics of T-junctions – an experimental investigation, Journal of Sound and Vibration 329 (2010) 1793–1808. [9] M. Karlsson, M. Åbom, On the use of linear acoustic multiports to predict whistling in confined flows, Acta Acustica United with Acustica 97 (2011) 24–33. [10] S. Föller, W. Polifke, D. Tonon, Aeroacoustic characterization of T-junctions based on large eddy simulation and system identification, Proceedings of the 16th AIAA/CEAS Aeroacoustics Conference, Stockholm, Sweden, 2010. [11] Bruggeman, J.C., Flow Induced Pulsations in Pipe Systems, PhD Thesis, Technical University of Eindhoven, 1987. [12] A. Michalke, On spatially growing disturbances in an inviscid shear layer, Journal of Fluid Mechanics 23 (3) (1965) 521–544. [13] M.S. Howe, The dissipation of sound at an edge, Journal of Sound and Vibration 70 (1980) 407–411. [14] J.C. Bruggeman, A. Hirschberg, M.E.H. van Dongen, A.P.J. Wijnands, Self-sustained aero-acoustic pulsations in gas transport systems: experimental study of the influence of closed side branches, Journal of Sound and Vibration 150 (3) (1991) 371–393. [15] G. Kooijman, A. Hirschberg, J. Golliard, Acoustical response of orifices under grazing flow: effect of boundary layer profile and edge geometry, Journal of Sound and Vibration 315 (2008) 849–874. [16] Belfroid, S.P.C., Schiferli, W., Peters, R.M.C.A.M., and Buffing., J., Flow induced pulsations caused by split flow in a T-branch connection, Proceedings of the Sixth FSI, AE & FIV þ N Symposium, 2006, Vancouver, Canada. [17] A. Holmberg, , M. Karlsson, M. Åbom, Aeroacoustics of rectangular T-junctions subject to combined grazing and bias flows – an experimental investigation, Journal of Sound and Vibration, 340 (31), 2015, 152–166, http://dx.doi.org/10.1016/j.jsv.2014.11.040. [18] P.K. Kundu, I.M. Cohen, Fluid Mechanics, Academic Press, USA, 2010. [19] S.B. Pope, Turbulent Flows, Cambridge University Press, England, 2000.
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
247
[20] R. Glav, M. Åbom, A general formalism for analyzing acoustic 2-port networks, Journal of Sound and Vibration 202 (5) (1997) 739–747. [21] S.W. Rienstra, A. Hirschberg, An Introduction to Acoustics, Eindhoven University of Technology, 2002, availabe at http://www.win.tue.nl/%7Esjoerdr/ papers/boek.pdf〉, corrections at 〈http://www.win.tue.nl/%7Esjoerdr/papers/corrections〉, last visited 17.06.13. [22] A. Kierkegaard, S. Boij, G. Efraimsson, Simulations of the scattering of sound waves at a sudden area expansion, Journal of Sound and Vibration 331 (2012) 1068–1083. [23] A. Kierkegaard, S. Allam, G. Efraimsson, M. Åbom, Simulations of whistling and the whistling potentiality of an in-duct orifice with linear aeroacoustics, Journal of Sound and Vibration 331 (2012) 1084–1096. [24] D. Kuzmin, O. Mierka, S. Turek, On the implementation of the k–ε turbulence model in incompressible flow solvers based on a finite element discretization, International Journal of Computing Science and Mathematics 1 (2) (2007) 193–206. [25] A.V. Johansson, P.H. Alfredsson, On the structure of turbulent channel flow, Journal of Fluid Mechanics 122 (1982) 295–314. [26] E. Salt, S. Mohamed, D. Arthurs, S. Ziada, Aeroacoustic sources generated by flow-sound interaction in a T-junction, Journal of Fluids and Structures 51 (2014) 116–131. [27] M.A. Leschziner., Modelling turbulent separated flow in the context of aerodynamic applications, Fluid Dynamic Research 28 (2006) 174–210.
Contents lists available at ScienceDirect
Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi
A frequency domain linearized Navier–Stokes method including acoustic damping by eddy viscosity using RANS Andreas Holmberg a,b,n, Axel Kierkegaard b, Chenyang Weng b,c a b c
KTH, CCGEx – the Competence Center for Gas Excange, Teknikringen 8, 100 44 Stockholm, Sweden KTH, Marcus Wallenberg Laboratory for Sound and Vibration Research, Teknikringen 8, 100 44 Stockholm, Sweden KTH, Linné Flow Centre, Teknikringen 8, 100 44 Stockholm, Sweden
a r t i c l e i n f o
abstract
Article history: Received 30 July 2013 Received in revised form 6 February 2015 Accepted 16 February 2015 Handling Editor: P. Joseph Available online 11 March 2015
In this paper, a method for including damping of acoustic energy in regions of strong turbulence is derived for a linearized Navier–Stokes method in the frequency domain. The proposed method is validated and analyzed in 2D only, although the formulation is fully presented in 3D. The result is applied in a study of the linear interaction between the acoustic and the hydrodynamic field in a 2D T-junction, subject to grazing flow at Mach 0.1. Part of the acoustic energy at the upstream edge of the junction is shed as harmonically oscillating disturbances, which are conveyed across the shear layer over the junction, where they interact with the acoustic field. As the acoustic waves travel in regions of strong shear, there is a need to include the interaction between the background turbulence and the acoustic field. For this purpose, the oscillation of the background turbulence Reynold's stress, due to the acoustic field, is modeled using an eddy Newtonian model assumption. The time averaged flow is first solved for using RANS along with a k-ε turbulence model. The spatially varying turbulent eddy viscosity is then added to the spatially invariant kinematic viscosity in the acoustic set of equations. The response of the 2D T-junction to an incident acoustic field is analyzed via a plane wave scattering matrix model, and the result is compared to experimental data for a T-junction of rectangular ducts. A strong improvement in the agreement between calculation and experimental data is found when the modification proposed in this paper is implemented. Discrepancies remaining are likely due to inaccuracies in the selected turbulence model, which is known to produce large errors e.g. for flows with significant rotation, which the grazing flow across the T-junction certainly is. A natural next step is therefore to test the proposed methodology together with more sophisticated turbulence models. & 2015 Elsevier Ltd. All rights reserved.
1. Introduction Duct junctions and sudden area changes are common elements in pipe networks. While sometimes added to scatter sound, e.g. expansion chambers in exhaust system configurations, especially T-junctions and wall corrugations are known to cause severe resonances, see e.g. [1–3]. The governing aeroacoustic phenomena are generally well known, but the prediction of frequency and magnitude of the acoustic field in the design stage is often beyond reach. Unfortunately, building and
n
Corresponding author at: KTH, CCGEx – the Competence Center for Gas Excange, Teknikringen 8, 100 44 Stockholm, Sweden. E-mail address: [email protected] (A. Holmberg).
http://dx.doi.org/10.1016/j.jsv.2015.02.030 0022-460X/& 2015 Elsevier Ltd. All rights reserved.
230
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
testing prototypes is a time-consuming and expensive process. There is therefore a demand for more efficient and accurate prediction tools. In this paper, we describe a numerically efficient method for calculating linear acoustic propagation, scattering and hydrodynamic interaction. The numerical method is a modification of a linearized Navier–Stokes equations method for the frequency domain, described by Kierkegaard et al. [4]. It is a hybrid approach, in which the flow is divided into a mean and a time dependent but coherent perturbation part. As the mean flow is calculated by RANS using a k–ε turbulence model and the perturbation equation set are linear, the method is very time efficient. The modification we propose begins with substituting the flow decomposition used in [4] by a triple decomposition, where the flow is divided into a time average, a phase average and a random part. This choice stems from the need to include the oscillation of the turbulence Reynold's stress, due to the acoustic field. This oscillation induces an additional acoustic dissipation process, which in general cannot be neglected [5]. The method of using a triple decomposition in order to include the turbulence damping in the linearized Navier–Stokes equations (LNSE) was first introduced in the Ph.D. thesis by Holmberg [6], which is the basis of this paper. Recently, Gikadi et al. [7] presented a paper similar to this, in which 3D LNSE calculations were performed, also using the triple decomposition to include the turbulence effect on the acoustic field. The background flow was in their work calculated using LES, while here RANS is used. We apply our method in the study of a T-junction of rectangular ducts, subject to grazing flow of Mach 0.1 across the opening, for which experimental results exist. The interaction between the acoustic field, the geometry and the shear layer across the junction is studied via a scattering matrix. From this matrix, the frequency ranges of amplification and attenuation can be predicted [8], as well as the frequencies of whistling [9]. Previously, the scattering matrix for several grazing flow Mach numbers was determined by Karlsson and Åbom [8] for cylindrical ducts. A numerical solution has also been obtained by Föller et al. [10], who superimposed plane waves on a large eddy simulation, for a grazing flow Mach number of 0.1. The only significant deviation between experimental and simulation results, was found in the scattering matrix element describing acoustic transmission from the upstream to the downstream main pipe. Föller et al. [10] argued that the much stronger transmission found in the experiments was unphysical. However, this large transmission is also found in a new set of experiments conducted for the rectangular duct junction studied here. We will show the Mach 0.1 grazing flow case, and compare to our numerical results. An interesting feature of the paper by Gikadi et al. [7] is that they also studied a T-junction – in their case a junction of circular pipes, hence the need for 3D simulations. The simulation results in [7] are compared to the experimental results of Karlsson and Åbom [8] – just as in this paper – for the case of Mach 0.1 grazing flow. In the conclusion section, their results will be discussed in light of the results presented in this paper. This paper is outlined as follows; in Section 2, the interaction between the hydrodynamic and acoustic fields in the Tjunction is described. The modification of the governing equations is derived in Section 3. The methods to analyze the calculated results are described in Section 4, while the numerical implementation of the equations and the analysis methods are discussed in Section 5. Results are shown and discussed in Section 6, while a summary of the paper and an outlook is presented in Section 7. 2. Theory of the T-junction When the acoustic particle velocity turns into the sidebranch at the upstream edge of the T-junction, part of the acoustic energy is shed as a small disturbance [11]. In the case of a mean flow in the system, a shear layer is formed across the Tjunction opening, and the disturbance is conveyed by the mean flow across this layer. Due to the instability of the shear layer, the disturbance grows as it travels across the junction [12]. The interaction between this disturbance, the mean flow and the acoustic field can be studied e.g. by the vortex-sound theory by Howe [13]. He expresses the transfer of power from the hydrodynamic field to the acoustic field via the integral Z Π ¼ ρ0 ω v u dV (1) where V is the volume of integration, ω is the distribution of vorticity and v and u are the incompressible and the acoustic velocity fields. As the integral is zero when u and v are parallel, the case of plane wave propagation in a straight duct yields no energy transfer. In the T-junction however, the angle between the particle velocity and the convection velocity is nonzero, and thus, energy transfer can occur, where the direction of transfer depends on the phase of the acoustic field. Since the disturbances are shed when the acoustic particle velocity turns into the sidebranch, the integral is positive during the first half of the acoustic period, thus resulting in a negative sound power contribution. During the second period, the acoustic particle velocity changes sign, resulting in a positive sound power contribution. A net increase is possible since the magnitude of the disturbance, shed by the acoustic wave grows as it is conveyed across the opening. Depending on the ratio between the acoustic and the hydrodynamic timescales, i.e. on the Storuhal number, there will be a net increase or decrease of the acoustic energy. For low acoustic velocity amplitudes, i.e. of the order of O( 3) and below when normalized by the mean flow velocity [14], the magnitude of the shed disturbance is invariant on (although correlated with) the incoming acoustic amplitude, which, according to Eq. (1) yields a linear relation between the incoming and amplified acoustic fields [14]. The rate of growth of the shed disturbance depends on the shear layer as well as the hydrodynamic wavelength of the disturbance. Michalke [12] investigated this for free inviscid parallel flows having a hyperbolic-tangent profile, using modal
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
231
stability analysis. His results are unfortunately not directly applicable to the T-junction, especially not for mixed grazing-bias flows, that is, flow in all three branches. Previous studies also clearly show that the power amplification depends strongly on e.g. the edges of the junction [15], the amount of grazing and bias flow [16,17] and the cross sections of the ducts [17]. It is not surprising therefore, that no reasonable versatile engineering models of the T-junction exist. 3. Derivation of the governing equations In this section, we present the derivation of a modified version of the frequency domain formulation of the linearized Navier–Stokes equations presented by Kierkegaard et al. [4]. As mentioned in the introduction, we propose this modification in order to include dissipation of the acoustic field in the presence of strong vorticity. As in [4], the starting point is the Navier–Stokes equations [18], here written as
ρ where
Dρ ∂u þρ i ¼ 0 Dt ∂xi
(2)
Dui ∂p ∂τij ¼ þ þ ρF i ∂xi ∂xj Dt
(3)
∂ui ∂uj 2 ∂uk þ δij ∂xj ∂xi 3 ∂xk
τij ¼ μ
(4)
Here, ρ, u and p, are the density, velocity and pressure fields. F is a volume force, μ the dynamic viscosity and δ the Kronecker delta function. By assuming an isentropic process, we can relate the density to the pressure, and decouple the equation system from the energy equation. Now, instead of dividing the flow field into a mean and a coherent perturbed field, as done in [4], we adopt the triple decomposition used by Reynolds and Hussain [5]. The velocity and pressure fields are here divided into three parts e i ðx; t Þ þ u0i ðx; t Þ ui ðx; t Þ ¼ ui ðxÞ þ u
(5)
0
e ðx; t Þ þ pðx; t Þ pðx; t Þ ¼ pðxÞ þ p
(6)
The last term vanishes when a phase average is calculated, and hence it is the turbulent part of the flow. The second term vanishes when a time average is calculated, and is therefore the coherent perturbed field, which includes the acoustic waves. The first term is simply the time averaged flow field. For the density we assume a simpler decomposition, namely
ρðx; t Þ ¼ ρ þ ρe ðx; t Þ
(7)
The turbulent part of the density is thus assumed negligible, and the time averaged density is assumed constant throughout the flow. In other words, the flow is approximated as a harmonic field superposed on to an incompressible flow field. This is likely to be accurate for low Mach numbers (below Mach 0.3 say). The mean static pressure is allowed to vary in
Fig. 1. A T-junction; definitions of coordinate systems, positive directions for the Mach numbers and numbering of side branches used in the paper.
Table 1 Physical size of the computational domain.
Duct length [mm] RANS Duct length [mm] LNSE Duct height/width [mm] Wave decomposition zone, distance from junction [mm] Buffer zone length [mm]
Branch I
Branch II
Branch III
1500 1500 25 300—1200 500
1500 1500 25 400–1200 500
1000 1500 25 400–1400 500
232
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
space due to velocity gradients, while the mean density is assumed constant throughout the flow. Now, applying Eqs. (5) and (7) to the continuity equation, we get e i þ u0i e ui þ u e ∂ ρþρ dρ þ ¼0 (8) ∂xi dt The solution we seek is the coherent perturbed field. Therefore we take the phase average of Eq. (8), linearize the result and introduce the assumption of monochromatic waves with ending up with e iωt time dependence. The result, transformed to the frequency domain, is e e ∂ρ ∂ui ∂u e ¼ ρ i ui þ iω ρ (9) ∂xi ∂xi ∂xi Before continuing with the momentum equation, we note that the turbulent part of Eq. (8) can be written as
ρ þ ρe
∂u0i e ∂ρ ¼ u0i ∂xi ∂xi
(10)
We now apply the flow decomposition, i.e. Eqs. (5)–(7) on the momentum equation, i.e. Eq. (3), and rearrange the result using Eq. (10). The result is ∂u þ u e i þu0i ei d u ∂u0i ∂u0i u0j e ∂ρ i e j þ u0j ej ρ þ ρ e e þ uj þ u ρ þ ρe þ ρþρ þu0i u0j þ uj þ u dt ∂xj ∂xj ∂xj ∂xj 0 1 0 e j þ u0j ∂ uj þ u e k þ u0k e i þ u0i e þp ∂ ρþρ ∂ @∂ u i þ u 2 ∂ uk þ u ¼ þ μ þ δij A þ ρ þ ρe F i (11) ∂xj 3 ∂xi ∂xj ∂xk ∂xi Now, nothing has been said about the volume force Fi yet. In this work we will implement it as a harmonically fluctuating e Fe i is zero, and we therefore volume force multiplied by the mean density. Hence, in the computations, the nonlinear term ρ e F i ¼ ρFe i . After introducing this implementation, we continue by subtracting the time average from the phase have ρ þ ρ
Fig. 2. Time average flow field: (a) velocity magnitude [m/s] and (b) eddy viscosity (νT) [m2/s].
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
233
average of Eq. (11). Linearizing the result thus yields e ∂e uj ∂e ρ uk ∂e ui ∂ 2 ∂e ρ ddteu i þ ue j ∂∂xeuji þ uj ∂∂xujj þ ρe uj ∂∂xeuji ¼ ∂x þ μ þ δ ∂xj 3 ∂xk ij ∂xj ∂xi i þ ρFei þ ρ
∂ u0i u0j 〈u0i u0j 〉 ∂xj
ρe
∂ 〈u0i u0j 〉
ρe
u0i u0j ∂ ∂xj ∂xj
! e ρ
∂〈u0i u0j 〉 ∂xj
(12)
We now have a closure problem due to the unknown terms ∂〈u0i u0j 〉 u0i u0j ∂e r ij x¼ ∂xj ∂xj
(13)
e e ρ ∂ρ ρe e ∂〈u0i u0j 〉 e ρ (14) r ij ¼ 〈u0i u0j 〉 þ u0i u0j ∂ ∂xj ∂xj ∂xj e ρ r ij . However, while this is true for small which must be modeled. To start with, we consider e r ij to be negligible compared to e e ρ e i =∂xj Þ þ perturbations, it should be kept in mind that e r ij is not necessarily smaller than the terms ð∂=∂xj Þμ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij Þ, which we have decided to keep. The reason why we keep the latter terms is that, depending on ð∂u how e r ij is modeled, they may not be negligible everywhere in the computational domain. Before continuing, it is interesting to note the physical interpretation of e r ij . Since e r ij is expressed as the difference between the time average and the phase average of the turbulent Reynold's stresses, it can be interpreted as the oscillation of the turbulent Reynolds's stress due to the acoustic field, as is also discussed in [5]. Now, in order to model e r ij , we start with the time average term. According to the Boussinesq hypothesis [19], the deviatoric Reynolds stress is proportional to the mean rate of strain 2 ∂ui ∂uj þ (15) ¼ 2νT Sij u0i u0j þ kδij ¼ νT 3 ∂xj ∂xi where νT is the eddy viscosity, and k is the turbulent kinetic energy. Now, for the phase average term, the idea of a Newtonian eddy model introduced in Reynolds and Hussain [5] could be used
Fig. 3. Perturbation density for sidebranch excitation: (a) 690 Hz and (b) 1370 Hz.
234
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
ei ej ek ∂ ui þ u ∂ uj þ u 2 ∂ uk þ u 2 þ δij þ kδij 3 3 ∂xj ∂xi ∂xk e e e ∂u ∂ u ∂ui ∂u i 2 ∂u k 2 j j ¼ νT þ þ δ þ þ kδij 3 ∂xj ∂xi ∂xj ∂xi 3 ∂xk ij 2 ¼ 2νT Sij þ e S ij þ kδij 3
〈u0i u0j 〉 ¼ νT
(16)
where the divergence of the time average velocity is set to zero, since it become zero if the time average of the continuity equation is linearized. Now, the above equation may be considered as improper, since the Boussinesq closure is an equilibrium expression that assumes the anisotropy tensor to depend only on the local instantaneous value of the mean strain rate tensor. The acoustic shear on the other hand, is a time dependent quantity, and thus the time response (phase) between the instantaneous value of the mean shear and of the acoustic shear should be taken into account. This is however not done here. Instead we assume that the effect of this is negligible. Thus, we get e j 2 ∂u e i ∂u ek ∂u e þ δij (17) r ij ¼ νT ∂xj ∂xi 3 ∂xk The eddy viscosity is readily obtained from the RANS calculations. Now remains the coupling between the pressure and the density. As previously mentioned, we here assume an isentropic process. The equation of state can thus be approximated by [21] Dp Dρ ¼ c2 Dt Dt
(18)
In all published work concerning the frequency domain linearized Navier–Stokes equations method [4, 22,23], the quiescent medium approximation of Eq. (18) is used, seemingly without any noticeable errors induced. We here chose to continue in the line of the previous work, due to the increased speed of the computations this approximation benefits from. Thus we write
ρe ¼ c2 ρe
e (x-component) for sidebranch excitation: (a) 690 Hz and (b) 1370 Hz. Fig. 4. Perturbation velocity u
(19)
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
The governing equations are finally e j 2 ∂u ∂u ei e e i ∂u ek e ∂u du ∂u ∂u ∂ρ ∂ e j i þ uj j þ ρ e uj i ¼ c2 þ þu ρ μ þ ρνT þ δij þ ρFei dt ∂xj ∂xj ∂xj ∂xi ∂xj ∂xj ∂xi 3 ∂xk ui
e e ∂ρ ∂ui ∂u e ¼ ρ i þ iω ρ ∂xi ∂xi ∂xi
235
(20)
(21)
Comparing these equations with those derived by Kierkegaard et al. [4], where the flow decomposition consisted of a division into a time average part and a coherent perturbed part, we can note that the only difference is the introduced eddy viscosity, which is added to the kinematic viscosity. 4. Analysis method In the case of a duct system where the gas flow only have linear effects on the acoustic field, N-port networks [20] can be used to calculate plane wave propagation and scattering. Each element in the network is then treated as an N-port – a linear time invariant system relating an input to an output. For T-junctions, orifice plates, shallow cavities and other elements in which strong shear layers are present in the regions of wave propagation, nonlinear coupling between the acoustic and the hydrodynamic fields may occur. Using an N-port network to study the installation effects of these elements is thus a questionable approach. However, recently N-ports determined in the low level – and thus linear – regime have been used to predict whistling frequencies of T-junctions [9] and orifice plates [23]. For planar waves, the N-port applied on the Tjunction is referred to as the 3-port scattering matrix, which is defined accordingly [8] 0 I 1 2 30 I 1 0 I 1s pþ pþ ρI τI;II τI;III p B II C 6 7B pII C B II C s B p þ C ¼ 4 τII;I B C ρ τ (22) II;III 5@ A þ II @ A @ p þ A or p þ ¼ Sp þ p þ pIII þ
τIII;I τIII;II
ρIII
pIII
pIII þ
e (y-component) for sidebranch excitation: (a) 690 Hz and (b) 1370 Hz. Fig. 5. Perturbation velocity v
236
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
where S is the scattering matrix, p- is the incoming acoustic waves, p þ is the outgoing waves, and psþ the sound that is generated independent on the incoming sound waves, see Fig. 1. The diagonal elements of the scattering matrix thus represent the reflection coefficients towards the T-junction, while the other elements correspond to transmission coefficients between the three branches. Another intrinsic property of the T-junction that can be studied using via the scattering matrix is the power balance between outgoing and incoming energy [8]. The scattering matrix can also be used to determine if a (linear) reflection connected to the T-junction is sufficient to generate a nonlinear response, resulting in a whistle, as was experimentally shown by Karlsson and Åbom [9]. Kierkegaard et al. [23] attempted to predict whistling using numerically obtained scattering matrices of orifice plates subject to bias flow. They found out however, that the predictions where quite sensitive to errors (3 percent) in the scattering matrices. Therefore a sensitivity analysis of the stability analysis should always be made. Due to its many uses it is clear that the scattering matrix is a powerful tool if calculated correctly. We will consider the error in the scattering matrix when the success of the suggested turbulence damping is discussed.
5. Implementation The calculation procedure consists of three main steps; first the time averaged part of the flow is solved for, using an incompressible, steady-state RANS method. Then the coherent perturbation equations are solved, and finally, the propagating waves are extracted and used to calculate the scattering matrix. In this work, we used Comsol 4.2 with a MatLab link. Due to computational cost, the computational realm is set to 2D, both for the RANS solution and the acoustic solution. The physical size of the T-junction is shown in Table 1. The meshes used in this work are described in detail in Appendix A.
5.1. RANS implementation The mean flow is calculated using incompressible RANS equations with k–ε turbulence modeling [19]. The reason for using the k–ε turbulence model is that it is one of the most common and simplest turbulence models found in commercial software. The turbulent viscosity is obtained from the turbulent kinetic energy k and the dissipation rate of turbulence 2 energy ε, using the model constant C μ ¼ 0:09, via the equation νT ¼ C μ ðk =εÞ.
Fig. 6. Magnitude of the scattering matrix elements at zero bias flow; black: experiments, red: calculations e r ij ¼ 0, blue: calculations e e i =∂xj Þþ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij . (For interpretation of the references to color in this figure legend, the reader is referred to the web r ij ¼ νT ð∂u version of this article.)
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
237
The time averaged solution used as input to the acoustic calculations in the next section is obtained without using wall functions. To achieve convergence however, wall functions [24] are used in two consecutive pre-steps; first, starting from a uniform inflow velocity profile, the flow is calculated in a duct of 30 m in order to obtain the outlet flow velocity, pressure, turbulence kinetic energy k, and the dissipation rate of turbulence energy ε. The outlet properties are then used as inlet condition in the main branch I, see Fig. 1. The second step involving wall-functions is a T-junction mean flow calculation, where the objective is to generate an initial condition of the flow, used to setup calculations without wall functions.
Fig. 7. Phase of the scattering matrix elements at zero bias flow; black: experiments, red: calculations e r ij ¼ 0, blue: calculations e i =∂xj Þ þ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij . (For interpretation of the references to color in this figure legend, the reader is referred to the web e r ij ¼ νT ð∂u version of this article.)
Fig. 8. Perturbation density for sidebranch excitation 1370 Hz (e r ij ¼ 0Þ.
238
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
5.2. NSE implementation The sound fields are excited by applying the source terms F in Eq. (28) to one duct branch at a time. This is done in order to determine the scattering matrix, which can only be done if there are at least three linearly independent sound fields. Hence, the perturbed field is solved for three times. The maximum for the continuous function for F is set 1500 mm from the junction in each branch, and the function is continuous and non-zero in a region of 0.05 m in both directions. The acoustic boundary layers are assumed to have negligible effects on the end result. Hence, we implement a slip condition for the perturbation velocity on all walls. In theory, the results obtained using the methodology suggested here, will not be influenced by the end reflections in the duct. However, for numerical stability reasons, the reflections should be kept to a minimum. We therefore apply non-reflecting boundary conditions at in- and outflow boundaries, which based upon the Helmholtz equations are formulated as e ¼ ikρ e; n^ U∇ρ
e ¼ iku e; n^ U ∇u
e ¼ ikv e n^ U∇v
(23a,b,c)
e is zero for vertical boundaries and u e is zero for horizontal boundaries. However, these boundary conditions break where v down in the presence of vorticity. For this reason, buffer zones with artificial viscosity are placed before each duct end. The mean flow is not calculated for the buffer zones instead the values of the mean flow quantities at the duct ends are extended to hold for the entire buffer zone as well. This is done in order to avoid acoustic scattering at the interface.
5.3. Wave decomposition To obtain the scattering matrix, the upstream and downstream propagating waves must be extracted from the solved perturbation fields. We here choose the acoustic density instead of acoustic pressure, but the choice is arbitrary, as the constant factor between the two quantities imposed by the isentropic relation disappears when calculating the scattering matrix. As we have implemented a wall slip condition for the perturbation velocity, there is no acoustic boundary layer, and assuming the acoustic damping due to bulk viscosity to be negligible we define the wavenumbers as being real-valued. We can then setup a nonlinear equation system with the real valued phase, wavenumber and amplitude of each wave as unknowns
ρe xj ¼ ρe þ eiðφ þ þ k þ xj Þ þ ρe eiðφ k xj Þ
(24)
The coordinate axis x is positive outwards and negative into the T-junction in each branch, and the scattering matrix is determined at the position xj ¼0. To solve the equation we apply a nonlinear curve-fitting algorithm, were for each xj, the input is the cross section mean value of the acoustic density. To suppress errors due to vorticity, the cross section mean is calculated at 30 points for each cross section, while the number of cross sections j is set to 150 per duct, in the wave decomposition zones shown in Table 1. In order to obtain a wave decomposition zone of 1000 mm for the sidebranch, it is made 500 longer for the acoustic domain than it is for the RANS calculations. The mean flow quantities taken at the cross section a distance of 1000 mm away from the junction, where the RANS domain ends, extended to hold for the extra duct length in the acoustic domain, just as was done for the buffer zones.
Fig. 9. Perturbation density for sidebranch excitation 1370 Hz (e r ij ¼ 0Þ, zoomed out view.
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
239
6. Results and discussion As already mentioned, the mean flow configuration we consider in this paper is that of a pure grazing flow of Mach 0.1, implying a Reynolds number of approximately 47600 based upon the duct height. To study the effects of the proposed modification of the numerical method, results are shown both with and without the modification, and are compared to experimental results. We also include a grid convergence study for both calculation methods and for the RANS calculations in Appendix A. The scattering matrix data is compared to the experimental results found in [17], for a T-junction consisting of rectangular ducts with a cross section of 120 25 mm2, subject to the same flow case as is studied here. While not being true 2D, the rectangular junction will to a large extent behave as a 2D edge in terms of coherent vortex shedding and acoustic propagation times. Concerning the cross-sectional dimensions of the experimental rig, it is argued in the experimental paper [17] that based on the findings by Alfredsson and Johansson [25], the width to height ratio is large enough for the rig to exhibit fully developed turbulent channel flows upstream of the T-junction. Another support for comparing 2D calculations with the rectangular rig experiments is found in the paper by Salt et al. [26], who performed PIV-measurements, acoustic measurements and FEM modelling on a T-junction with a width/height ratio of 117/70 in the sidebranch and 127/70 in the main branches. They conclude that neglecting the fluctuating component in the third dimension only have minor effects on the flow-acoustic coupling mechanism. As the width/height ratio is smaller in the paper by Salt et al. [26] than for the case compared to here [17], the experimental results are expected to essentially behave as if the rig was 2D. The magnitude of the velocity field and the turbulent viscosity, obtained from the final RANS calculation are shown in Fig. 2, where the distinct shear layer across the opening can be observed. There is also a clear recirculation in the sidebranch. It is interesting as well to note that even in the center of the duct, the turbulent viscosity is two orders of magnitude larger than the kinematic viscosity. A few samples of typical coherent perturbation fields are given in Figs. 3–5, where the perturbation density and velocity components are shown at 690 Hz (a Strouhal number of 0.50) and at 1370 Hz (a Strouhal number of 1.0), for sidebranch excitation (branch III). In these figures, we have used the closure assumption e e i =∂xj Þ þ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij . The resulting scattering matrices from this assumption, as well as data r ij ¼ νT ð∂u obtained using e r ij ¼ 0, are shown together with experimental data in Figs. 6–7. When approaching the Strouhal number of unity, a strong peak can be observed in the magnitude of the calculated e i =∂xj Þ þ ð∂u e j =∂xi Þ scattering matrix obtained setting e r ij ¼ 0. This peak disappear using the model e r ij ¼ νT ð∂u e k =∂xk Þ δij Þ. We can also note that at Strouhal numbers less than 0.5 say, the results using ð2=3Þð∂u e e i =∂xj Þ þ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij is better than when setting e r ij ¼ νT ð∂u r ij ¼ 0, while there is a midrange of Strouhal numbers where setting e r ij ¼ 0 results in almost perfect agreement (of the magnitude) to the experimental data. In order to further study the peak at the Strouhal number of approximately 1.0, we present two figures capturing the density perturbations at 1370 Hz for sidebranch excitation and the case where e r ij ¼ 0, see Figs. 8 and 9. In these figures strong instabilities are clearly observed which are convected with the time averaged flow, while at the same time forming a wavelike pattern. Their presence stems from an unbalance in perturbation energy, due to the lack of a proper dissipation term in the perturbation equations. Although this effect could probably be treated to some extent by artificially increasing the kinematic viscosity, it is more physically relevant to include the effect of the Reynolds stresses via the eddy viscosity. Mathematically the difference is that the eddy viscosity is spatially varying, which implies that the turbulence damping can be zero also when the shear of the perturbed quantities is not, and vice versa. 7. Conclusion Propagation, scattering and hydrodynamic interaction of acoustic waves have been calculated using a linearized Navier– Stokes method in the frequency domain. The decomposition of the flow into time average, phase average and random parts yields an additional term in the linearized equation set. This term was modeled by the Newtonian eddy model [5], i.e. proportional to the turbulent viscosity. The importance of retaining and modeling this term is studied for the case of a rectangular T-junction subject to grazing flow of Mach 0.1. It should be noted that the validation and analysis of the proposed method is in this paper limited in 2D, although the formulation is fully presented in 3D. Comparing to experimental data, the prediction of the acoustic-hydrodynamic interaction is several factors stronger, unless the modification of the equations suggested in this paper is implemented. As the turbulent viscosity is already solved for during the RANS calculations and simply added to the kinematic viscosity, there is no additional computation time involved. This modified version of the method presented by Kierkegaard et al. [4] is therefore very promising. A natural next step is to try other viscosity models, since the k-ε model is known to be inaccurate for e.g. flows with normal-stress anisotropy, high strain rates and curvatures [27]. Comparing with the results presented by Gikadi et al. [7], it is clear that their results support the observations made here; i.e. that the inclusion of the turbulence damping on the perturbation quantities yields stable and physical relevant results, while neglecting the turbulence effect may cause severe over predictions (mostly) of the acoustic amplification by the hydrodynamic shear flow. The results by Gikadi et al. [7] also shows frequencies where a small error is induced by including the turbulence damping, although at a very few number of frequency points compared to what is found here. It is not surprising that their results are more accurate however, since they solve the background flow using an LES in 3D, while we use 2D RANS with k-ε. The method to choose depends on the accuracy required. We note though that the method
240
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
presented here is easily extended in 3D, which would be faster than the LES method presented by Gikadi et al. [7]. It would therefore be a very interesting case to compute, preferably with several different turbulence models.
Appendix A In this appendix, the main results of conducted grid convergence studies are presented. For the RANS calculations, three different meshes were tested, see Table A1. In all three meshes, quadrilateral elements were used below yoh/25, where y is the wall normal direction and h¼0.025 m is the duct height. The height of the wall-adjacent element was kept constant as the mesh was refined. Everywhere else in the computational domain, triangular elements were used. The convergence was studied by plotting the x-component of the time averaged velocity at three different lines, all with the length of the duct height, see Fig. A1. The velocity profiles at the three lines are shown in Figs. A2–A4. In Fig. A5, the magnitude of the calculated scattering matrices are shown for each RANS mesh, where the acoustic closure e i =∂xj Þ þð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij is used. The finest acoustic mesh, see Table A2, was applied in hypothesis e r ij ¼ νT ð∂u this study. The velocity profiles obtained using the medium and fine RANS meshes are almost identical, while the course RANS mesh yields slightly different profiles. The influence of this on the acoustic result is in Fig. A5 shown to be tremendous. There is also a slight variation between the results obtained by the medium and the fine RANS mesh, which is seen in the highest and the lowest Strouhal number ranges. The fine RANS mesh solution was applied in the end, and was also used in the acoustic mesh convergence study. A closeup of the medium acoustic mesh, which was used for the final acoustic results in this paper, is shown in Fig. A6. The mesh consists of three main parts. At the downstream edge, a quadrilateral mesh was applied in a rectangle up to a wall normal direction of yoh/20, and with a streamwise length of h/2. Far from the T-junction, in the plane wave propagation regions of the ducts, a quadrilateral mesh was also applied. In between these two regions, a triangular mesh was used. The number of triangular elements across the duct opening, the number of wall normal quadrilateral elements, the maximum growth rate and the maximum element size was specified for three different meshes, see Table A2. In the plane wave propagation regions (far field), the wall normal size ratio between the largest element and the smallest, wall adjacent element was set to 3. A variant of the medium mesh was also tested, by setting this size ratio to 1, see Fig. A7. The result of the convergence study (for the size ratio set to 3) is shown Figs. A8 and A9, in terms of the magnitude of the scattering matrix elements at five frequencies, viz. 100, 400, 700, 1000 and 1400 Hz. The sensitivity to the wall normal size ratio is investigated by computing the scattering matrix in the full frequency range of interested. The result of this study is shown in Fig. A10 for the case without turbulent damping, and in Fig. A11 for the case with this damping implemented.
Table A1 RANS mesh properties. Mesh
Coarse
Medium
Fine
Triangular elements Quadrilateral elements Total number of elements Max. element size Min. element size Max. element growth Number of elements at yo h/25
445,676 329,840 775,516 h/25 h/2500 1.025 40
258,344 164,920 423,264 h/25 h/2500 1.05 20
156,301 491,60 205,461 h/15 h/2500 1.05 10
Fig. A1. Figure showing the lines at which the x-component of the time average velocity is plotted. .
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
241
Fig. A2. Time averaged velocity component in the x-direction at cut line A. Blue dashed line: coarse RANS mesh, red line: medium RANS mesh, green line: fine RANS mesh. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
Fig. A3. Time averaged velocity component in the x-direction at cut line B. Blue dashed line: coarse RANS mesh, red line: medium RANS mesh, green line: fine RANS mesh. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
Fig. A4. Time averaged velocity component in the x-direction at cut line C. Blue dashed line: coarse RANS mesh, red line: medium RANS mesh, Green line: fine RANS mesh. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
242
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
e i =∂xj Þþ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij , black: coarse RANS mesh, red: Fig. A5. Magnitude of calculated scattering matrix elements, using e r ij ¼ νT ð∂u medium RANS mesh, blue: fine RANS mesh. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
Table A2 Acoustic mesh properties. Mesh
Coarse
Medium
Fine
Triangular elements Quadrilateral elements Max. element size Max. growth rate Elements across opening No. of wall normal elements at downstream edge No. of wall normal elements in plane wave region Approximate computation time per frequency point on standard laptop PC
2042 6005 h/5 1.2 25 3 5 4
7668 23,900 h/10 1.1 50 5 10 5
29,689 95,600 h/20 1.05 100 10 20 45
Fig. A6. Close up of the medium mesh used for acoustic calculations. .
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
243
Fig. A7. Close up of the medium mesh without wall normal grid stretching in the far field. .
e i =∂xj Þ þ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij , black: coarse acoustic mesh, red: Fig. A8. Magnitude of calculated scattering matrix elements, using e r ij ¼ νT ð∂u medium acoustic mesh, blue: fine acoustic mesh. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
244
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
Fig. A9. Magnitude of calculated scattering matrix elements, using e r ij ¼ 0, black: coarse acoustic mesh, red: medium acoustic mesh, blue: fine acoustic mesh. For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
245
Fig. A10. Magnitude of calculated scattering matrix elements, using e r ij ¼ 0, black: experimental results, red: no wall normal grid stretching in the far field, see Fig. A7, blue: wall normal grid stretching in the far field applied, see Fig. A6. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
246
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
e i =∂xj Þ þ ð∂u e j =∂xi Þ ð2=3Þð∂u e k =∂xk Þ δij , black: experimental results, red: no Fig. A11. Magnitude of calculated scattering matrix elements, using e r ij ¼ νT ð∂u wall normal grid stretching in the far field, see Fig. A7, blue: wall normal grid stretching in the far field applied, see Fig. A6. Note that the blue curve lies almost exactly atop the red curve. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) .
References [1] S. Ziada, Flow-excited acoustic resonance in industry, Journal of Pressure Vessel Technology 132 (2010) 1–9. [2] P. Lafon, S. Caillaud, J.P. Devos, C. Lambert, Aeroacoustical coupling in a ducted shallow cavity and fluid/structure effects on a steam line, Journal of Fluids and Structures 18 (2003) 695–713. [3] S.P.C. Belfroid, D.P. Shatto, R.M.C.A.M. Peters, Flow induced pulsations caused by corrugated pipes, Proceedings of ASME Pressure Vessel and Piping Division Conference, 2007, San Antonio, Texas. [4] A. Kierkegaard, S. Boij, G. Efraimsson, A frequency domain linearized Navier–Stokes equations approach to acoustic propagation in flow ducts with sharp edges, Journal of Acoustic Society of America 127 (2) (2010) 710–719. [5] W.C. Reynolds, A.K.M.F. Hussain, The mechanics of an organized wave in turbulent shear flow. Part 3. Theoretical models and comparisons with experiments, Journal of Fluid Mechanics 54 (2) (1972) 263–288. [6] Holmberg, A., Aeroacoustic Characterization using Multiport Methods, PhD Thesis, Royal Institute of Technology, Stockholm, 2012. [7] J. Gikadi, S. Föller, T. Sattelmayer, Impact of turbulence on the prediction of linear aeroacoustic interactions: acoustic response of a turbulent shear layer, Journal of Sound and Vibration 333 (24) (2014) 6548–6559, http://dx.doi.org/10.1016/j.jsv.2014.06.033. [8] M. Karlsson, M. Åbom, Aeroacoustics of T-junctions – an experimental investigation, Journal of Sound and Vibration 329 (2010) 1793–1808. [9] M. Karlsson, M. Åbom, On the use of linear acoustic multiports to predict whistling in confined flows, Acta Acustica United with Acustica 97 (2011) 24–33. [10] S. Föller, W. Polifke, D. Tonon, Aeroacoustic characterization of T-junctions based on large eddy simulation and system identification, Proceedings of the 16th AIAA/CEAS Aeroacoustics Conference, Stockholm, Sweden, 2010. [11] Bruggeman, J.C., Flow Induced Pulsations in Pipe Systems, PhD Thesis, Technical University of Eindhoven, 1987. [12] A. Michalke, On spatially growing disturbances in an inviscid shear layer, Journal of Fluid Mechanics 23 (3) (1965) 521–544. [13] M.S. Howe, The dissipation of sound at an edge, Journal of Sound and Vibration 70 (1980) 407–411. [14] J.C. Bruggeman, A. Hirschberg, M.E.H. van Dongen, A.P.J. Wijnands, Self-sustained aero-acoustic pulsations in gas transport systems: experimental study of the influence of closed side branches, Journal of Sound and Vibration 150 (3) (1991) 371–393. [15] G. Kooijman, A. Hirschberg, J. Golliard, Acoustical response of orifices under grazing flow: effect of boundary layer profile and edge geometry, Journal of Sound and Vibration 315 (2008) 849–874. [16] Belfroid, S.P.C., Schiferli, W., Peters, R.M.C.A.M., and Buffing., J., Flow induced pulsations caused by split flow in a T-branch connection, Proceedings of the Sixth FSI, AE & FIV þ N Symposium, 2006, Vancouver, Canada. [17] A. Holmberg, , M. Karlsson, M. Åbom, Aeroacoustics of rectangular T-junctions subject to combined grazing and bias flows – an experimental investigation, Journal of Sound and Vibration, 340 (31), 2015, 152–166, http://dx.doi.org/10.1016/j.jsv.2014.11.040. [18] P.K. Kundu, I.M. Cohen, Fluid Mechanics, Academic Press, USA, 2010. [19] S.B. Pope, Turbulent Flows, Cambridge University Press, England, 2000.
A. Holmberg et al. / Journal of Sound and Vibration 346 (2015) 229–247
247
[20] R. Glav, M. Åbom, A general formalism for analyzing acoustic 2-port networks, Journal of Sound and Vibration 202 (5) (1997) 739–747. [21] S.W. Rienstra, A. Hirschberg, An Introduction to Acoustics, Eindhoven University of Technology, 2002, availabe at http://www.win.tue.nl/%7Esjoerdr/ papers/boek.pdf〉, corrections at 〈http://www.win.tue.nl/%7Esjoerdr/papers/corrections〉, last visited 17.06.13. [22] A. Kierkegaard, S. Boij, G. Efraimsson, Simulations of the scattering of sound waves at a sudden area expansion, Journal of Sound and Vibration 331 (2012) 1068–1083. [23] A. Kierkegaard, S. Allam, G. Efraimsson, M. Åbom, Simulations of whistling and the whistling potentiality of an in-duct orifice with linear aeroacoustics, Journal of Sound and Vibration 331 (2012) 1084–1096. [24] D. Kuzmin, O. Mierka, S. Turek, On the implementation of the k–ε turbulence model in incompressible flow solvers based on a finite element discretization, International Journal of Computing Science and Mathematics 1 (2) (2007) 193–206. [25] A.V. Johansson, P.H. Alfredsson, On the structure of turbulent channel flow, Journal of Fluid Mechanics 122 (1982) 295–314. [26] E. Salt, S. Mohamed, D. Arthurs, S. Ziada, Aeroacoustic sources generated by flow-sound interaction in a T-junction, Journal of Fluids and Structures 51 (2014) 116–131. [27] M.A. Leschziner., Modelling turbulent separated flow in the context of aerodynamic applications, Fluid Dynamic Research 28 (2006) 174–210.