Localization of aeroacoustic sound sources in viscous flows by a time reversal method

Journal of Sound and Vibration 332 (2013) 3655–3669

Contents lists available at SciVerse ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Localization of aeroacoustic sound sources in viscous flows by a time reversal method Philippe Druault a,b,n, Re´gis Marchiano a,b, Pierre Sagaut a,b a b

University Pierre and Marie Curie- Paris 6, UMR 7190, Institut Jean Le Rond d’Alembert, 75252 Paris cedex 5, France CNRS, UMR 7190, Institut Jean Le Rond d’Alembert, 75252 Paris cedex 5, France

a r t i c l e i n f o

abstract

Article history: Received 7 November 2011 Received in revised form 19 November 2012 Accepted 4 February 2013 Handling Editor: P. Joseph Available online 13 March 2013

A time reversal procedure used previously for aeroacoustic source localization is extended to the viscous flows case. Successive flow configurations including steady flows and unsteady ones have been performed to test and validate the proposed procedure. It is demonstrated that it provides a useful methodology for the detection and the localization of aeroacoustic sources, not only in presence of a dissipative medium but also in unsteady flow configurations. When dealing with dissipative media, the application of time-reversal procedure allows for the recovery of both the shape and the location of the source of sound even if the source amplitude is damped due to the viscous energy loss. For unsteady flows (in presence or not of a dissipative medium), the aeroacoustic source detection method remains always effective even if some theoretical assumptions are not fully satisfied. & 2013 Elsevier Ltd. All rights reserved.

1. Introduction The localization of aeroacoustic sources from the knowledge of the far-field pressure remains one of the outstanding problems in the aeroacoustic community. One of the main difficulties is that different kinds of sound sources may induce the same far acoustic pressure field. The non-uniqueness of the solution in that kind of inverse problem explains on the one hand the difficulty of that problem and on the other hand the variety of methods proposed to solve it such as acoustic analogy approach, causality approach or direct resolution of the inverse problem. The first group of method relies on the famous Lighthill’s analogy [1] or improved versions [2–4]. These approaches are based on the rearrangement of the flow equations leading to the appearance of source terms referred to as the Lighthill tensor. The application of such analogies requires to choose the source term model along with the propagation operator. Therefore, an a priori knowledge about the location and the nature of the acoustic sources is assumed. For instance, in absence of external forces and for isentropic flows, the source term corresponds to the double divergence of the Lighthill tensor and is therefore of quadrupole type source. Then, even if previous approaches for aeroacoustic source investigation relied mostly on acoustic analogy approaches, they may suffer from limitations and ambiguities. For instance, previous turbulent jet noise investigations have not directly demonstrated that the noise sources actually consist in quadrupoletype sources. In this sense, based on acoustic analogy, equivalent mathematical source may then not correspond to real physical aeroacoustic sources. The second group of aeroacoustic source localization procedures is based on the concept of causality. It consists in extracting the flow region which exhibits the highest correlation with the far acoustic pressure signal. In such a way, it is

n

Corresponding author. Tel.: þ33 1 44 27 54 72. E-mail address: [email protected] (P. Druault).

0022-460X/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jsv.2013.02.006

3656

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

assumed that the source is defined in a statistical sense. These approaches may suffer from limitations due to the difficulty in well extracting the aerodynamic events or aerodynamic regions which are best correlated to the far pressure field [5,6]. The other procedures are based on the resolution of the inverse problem. For such investigations, several methods can be applied. For instance, the adjoint method [7,8] is a promising method but requires important storage capacities. Another class of methods consists in solving the inverse problem using the time reversal method [9]. Time reversal method is a powerful method to solve inverse problems. The key idea is to use the reversibility of equations to back propagate information in time and retrieve original state for a given configuration. Hence, time reversal method has been widely used in acoustics for problems involving imaging, source localization, wavefront synthesis or communications [10]. The basic assumptions require that equations be symmetric for the time variable (invariance with respect to time) and that information can be collected on a surface enclosing the domain of interest. Those assumptions correspond to the classical time reversal cavity [11]. If the information is collected only on a portion of the surface enclosing the domain of interest, the method is called Time Reversal Mirror. The method can be applied and leads to good results slightly degraded compared to time reversal cavity [11]. This is a first example of the robustness of the time reversal method even though it is applied in configurations which do not completely fulfill the basic assumptions. Time reversal method relies mainly on two basic assumptions conservation of energy and reciprocity of the medium [12–14]. Conservation of energy ensures that all the energy can be reversed to the original state. So, it excludes any dissipation process. Reciprocity means that a signal emitted at point A and received at point B is the same as a signal emitted at point B and received at point A. In terms of operators, those two properties are equivalent to the fact that the operator modeling the problem is self-adjoint. For instance, that property is true for the wave equation which is the operator describing the propagation of wave in homogeneous lossless media [15]. It is also true for wave equation in heterogeneous lossless media [15]. In the framework of aeroacoustics, the presence of flows breaks the property of reciprocity and flows are generally viscous. Indeed, even for a simple configuration such as the propagation of a pulse between two points (namely A and B) in presence of a uniform flow along direction AB, the signal emitted at point A and received at point B, is not the same as a signal emitted at point B and received at point A. So, the viscosity of the medium and the presence of flows do not allow to use time reversal method stricto sensu in aeroacoustics. To apply it in aeroacoustics, one has to adapt the method to recover reciprocity, and to exploit the intrinsic qualities of robustness of time reversal method. In a previous paper, Deneuve et al. [9] showed that it was possible to use time reversal method in order to localize the region where sound was generated. Time reversal method was applied to configurations for which the source was created by injection of mass, vibration of surfaces and noise generated by flow itself. In that previous paper, authors deal with Euler equations without dissipation in the medium. The proposed methodology overcomes the problem of reciprocity by a numerical approach consisting in computing the solution thanks to a medium where flows are reversed (in direction). Recently, the time reversal technique was also performed by the Lattice Boltzmann method [16]. Nevertheless, the question of the applicability of the method in presence of dissipative media has not been yet treated in the context of aeroacoustics. The goal of this paper is to investigate carefully these two problems (reciprocity and dissipation) and to propose solutions to overcome them and to be capable to use the extended time reversal method in flows without solid boundaries. In the present paper, the previous methodology is analyzed more deeply. In particular, the theoretical framework of time reversal in aeroacoustics for viscous medium is derived. Depending on the modeling equations it is shown that three cases are possible: (i) configurations for which time reversal method can be fully applied, (ii) configurations for which time reversal method can be applied but yields approximate solutions only and (iii) configurations for which time reversal method cannot be applied theoretically. Numerical simulations are proposed to support these theoretical results. Moreover, in the third case, it is numerically observed that for some configurations even though the time reversal method cannot be applied theoretically, it gives reliable results. The first part of the paper is devoted to the theoretical analysis of time reversal in aeroacoustics. The three above cases are discussed and precise rules are derived to use time reversal method in aeroacoustics. Description of the numerical algorithms and the methodology are discussed in the second part. In the third part, numerical simulations involving the three kinds of configurations are presented and discussed. 2. Time reversal in aeroacoustics: theoretical analysis As mentioned in the Introduction, two assumptions are necessary to use time reversal method for a given set of equations. First of all, the reciprocity relation is essential. Reciprocity is characterized by the possibility to exchange emitter and receiver. Obviously, the presence of flows breaks that property [17]. Let us imagine an experience of source localization where the initial configuration is a Gaussian spot of mass injected in a medium with a uniform flow. During propagation, the Gaussian pulse propagates and is distorted according to mean flow effect. Using the direct time reversal method, the wave will not back propagate to the original source point. Because of the flow, a shift will appear between the founded position and the original one. That shift depends on the magnitude of the flow. That last property can be used to determine characteristics of the flow as the mean flow or vorticity [17,18]. For configurations with steady flows, it is possible to retrieve reciprocity by inverting the direction of the flows. That relation known as reverse flow reciprocity was introduced by Howe [4] and used recently by Deneuve et al. [9]. Even if it seems very delicate in an experiment, from

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

3657

a numerical point of view, it consists only in inverting time and flow direction ðvðx,tÞ-vðx,tÞÞ between the direct simulation and the time reversed one. Thus, applying reverse flow reciprocity allows to recover an essential property required by time reversal method. Note that it is not always possible to apply reverse flow reciprocity, in particular, it holds only for steady flows. Consequently, time reversal method will be valid only for this kind of flows. That limitation seems to prohibit the application to configurations where sound is generated either by turbulent flows or by the interaction of the flow with solid boundaries. The latter configuration is beyond the scope of this paper but the former one is examined in the last section of this paper. The second assumption is the conservation of energy. Let us imagine waves propagating in a dissipative medium during the direct way evolution. In a time reversal experiment (numerical or experimental one), waves are back-propagated in the medium and undergo the dissipative effects twice. Note that dissipation does not break reciprocity. Indeed, in the last example, emitter and receiver can be exchanged. The effects due to dissipation are a loss of amplitude and a degradation of the original wave field. In acoustics, it has been shown that the quality of focusing using time reversal method is degraded [15]. The degradation consists in a reduction of the amplitude of the main lobe compared to a medium without dissipation but also in an increase of the levels of side lobes. Nevertheless, the maximum of amplitude is well-localized, even though energy conservation is broken. That is an illustration of the robustness of the time reversal method. The robustness is due to the fact that the time reversal method acts as a spatio-temporal matched filter [19]. That property means that even if energy is lost, using time reversal method ensures to back-propagate the maximum of amplitude at the good point x0 and at the correct time Tt 0 where x0 is the original point of emission, t0 is the original time of emission, and T is the time for which the time reversal method is applied. Note that the result is only valid for point x0 . There is no guarantee on the control for the other points of the wave field, that is why side lobes may increase [19]. Hence, the matched filter property is very important and justifies the use of time reversal method in absorbing media. It shows that localization will be correct (in both time and space) even though the absolute level and eventually the side lobes will be degraded. From the above considerations, three cases must be distinguished: 1. configurations with steady flows and no dissipation; 2. configurations with steady flows and dissipation; 3. configurations with unsteady flows. In the first case, the time reversal method can be fully applied. Assumption of energy conservation is correct, and the assumption of reciprocity can be made by using the reverse flow reciprocity property. Consequently, in such a case, the time reversal method will retrieve the source with the correct localization and correct level. In the second case, the spatio-temporal matched filter property is essential. Indeed, that property ensures that sources localization and time of emission of sound sources will be retrieved even if a slight degradation of the original state will occur. Consequently, the time reversal method can be applied to localize position and time of emission of sound sources to all configurations with steady flows even if dissipation exists. Again, it is necessary to ensure the reciprocity by reversing the flows. In the third case, it is not possible to fully justify theoretically the use of the time reversal method. Indeed, reciprocity cannot be retrieved by applying the flow reverse technique. That holds for viscous or inviscid flows. Nevertheless, as will be shown in the section devoted to numerical results, in some cases it is reasonable to use that method even though it is not theoretically justified.

3. Methodology for sound sources localization 3.1. General methodology The methodology is built in two steps: (i) numerical computation of the direct problem, (ii) numerical computation of the inverse problem to localize sources by a time reversal method and a sensitivity analysis. Then, two numerical simulations have to be performed. The first stage deals with the direct computation of the aeroacoustic problem. The numerical method detailed below is used to compute the flow variables at each point of a numerical domain for each time step of the simulation. During this stage, the flow variables are stored only at the boundaries of that domain for each time step. Besides, at the final time of the computation, flow variables are stored in the whole computational domain. Storage capacities are then limited compared to other methods like adjoint approach. The second stage consists in solving the inverse problem thanks to the time reversal procedure and to perform a sensitivity analysis with the complex differentiation technique [20]. The stored data in the first stage are time and flow reversed. Then, they are used as inflow boundary conditions for the second numerical simulation. In parallel, a sensitivity analysis about acoustical disturbances measured at the boundaries of the computational domain is made. The sensitivity is not computed with numerical derivation but with a smart complex differentiation technique. This technique permits to compute the sensitivity to acoustical disturbances just by monitoring the imaginary part of the flow variables. The analysis is performed simultaneously with the second numerical simulation.

3658

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

3.2. Numerical solver for the Navier–Stokes equations To demonstrate the capabilities of the extended time reversal method, it is chosen to use it in conjunction with a twodimensional solver of the Navier–Stokes equations. This choice allows to address problems of sound propagation as well as sound generation in many kinds of flows. The study is here restricted to the Newtonian fluids with the Stokes’ assumption ! ! ! and the viscous tensor is assumed to be m½r v þðr v ÞT 2=3r: v 1 (where m is the dynamic viscosity and 1 is the unity tensor). It is well known that this thermo-viscous absorption is weak for acoustic waves at moderate frequencies. Most of the time, other dissipation processes are stronger for acoustical waves such as molecular relaxation. Nevertheless, we choose to use only this term to introduce dissipation. Even though this attenuation is not due to thermo-viscous attenuation in real configurations, we use it to mimic eventual other sources of dissipation by overestimating the value of the parameter m in order to have only one parameter controlling the dissipation. A numerical code solving the two-dimensional Navier–Stokes equations was previously developed for aeroacoustic analysis [9,21]. After recalling the pseudo-characteristic formulation [22] of the Navier–Stokes equations in the dimensionless form, the numerical method is presented with a special emphasis on boundary conditions for viscous inflow. To investigate the effect of the viscosity onto the aeroacoustic source localization method, it is necessary to use the dimensionless version of the pseudo-characteristic formulation of the Navier–Stokes equations. To this end, the following dimensionless variables for the velocity, length, density and temperature are defined, respectively, as follows: un ¼ u=U 0 , vn ¼ v=U 0 , xn ¼ x=L0 ðyn ¼ y=L0 Þ, rn ¼ r=r0 , mn ¼ m=m0 and T n ¼ T=T 0 . U0 represents a reference velocity like for instance the convection velocity when dealing with plane mixing layer flow configuration (see Section 4.2). Then using previous flow dimensionless variables mentioned with a 0 index, additional characteristic variables (time, pressure, speed of sound, entropy, viscosity and thermal conductivity) are then deduced: t 0 ¼ L0 =U 0 , p0 ¼ r0 U 20 , c20 ¼ gU 20 (where g is the heat g capacity ratio), s0 =Cv ¼ p0 r0 , m0 ¼ r0 U 0 L0 and k0 ¼ m0 C p =Pr. In this work, m0 corresponds to the physical dynamic viscosity of air at about T ¼ 25 1C and p0 ¼ 101 325 Pa. It is worth noting that the dimensionless viscosity is the inverse of the Reynolds number Ren ¼ r0 U 0 L0 =m. The two-dimensional version of the non-dimensionalized equations of the pseudo-characteristic formulation becomes    pn s0 qsn qpn rn cn  n þ ns ns X þ X n þY n þ þ Y n þ n ¼ n þ X þY qt 2 C v qt   qun 1 nþ 1 ðX X n Þ þY nu þ X nm n ¼ 2 qt r   qvn 1 nþ 1 nv n ðY ¼  X þ Y Þ þ Y nm 2 qt n r qsn R ¼ X ns Y ns þ n Snm p qt n in which X n 7 ¼ ðun 7 gcn Þ Y n 7 ¼ ðvn 7 gcn Þ

Y

with U 2 mn F ¼ 0 s0 T 0 3 n

1 qpn qun 7 rn cn qxn qxn



1 qpn qvn 7 r qyn qyn

 

n cn

(2)

(3)

qvn , qxn

X ns ¼ un

qsn qxn

(4)

Y nu ¼ vn

qun , qyn

Y ns ¼ vn

qsn qyn

(5)

¼

S



X nv ¼ un

X nm ¼

nm

(1)

nm

mn 3

mn 3

4

q2 un q2 un q2 vn þ 3 n2 þ n n n 2 qx qy qx qy

q2 vn q2 vn q2 un 4 n2 þ 3 n2 þ n n qy qx qy qx

! (6) !

! 1 n q2 T n q2 T n ¼ k þ n2 þ Fn s0 qxn2 qy

 n 2  n 2  n 2 ! qu qv qv U 2 mn 4 þ4 þ3 þ 0 qxn qyn qxn s0 T 0 3

!  n 2 qu qvn qun qun qvn 3 þ6 n n 4 n n qyn qx qy qx qy

(7)

(8)

(9)

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

3659

One of the main interests of the pseudo-characteristic formulation is that it enables an easy implementation of efficient boundary conditions separating acoustic and hydrodynamic disturbances. Two kinds of boundary conditions are involved in the present study: inflow and outflow conditions. The implementation of outflow condition is straightforward by using the characteristic formulation. Indeed, the outflow conditions are equivalent to set the inflow fluxes to zero. For the left and right boundaries, X n þ ¼ 0 and X n ¼ 0 are prescribed, respectively, and for the bottom and top boundaries: Y n þ ¼ 0 and Y n ¼ 0. The inflow condition specification is often used when flow variables and their associated fluxes are known along a specific boundary. In the homoentropic case, the inflow condition implementation is given in the previous paper [9]. We then propose to detail these inflow conditions when dealing with Navier–Stokes equations. As a function of the direction treatment, two different inflow condition implementations are considered. Left and right inflow boundary condition. Taking left and right boundaries of the computational domain, we assume that the time evolution of the flow variables (un, vn, pn, sn, Tn) is available. Then, the following temporal and spatial derivatives: qpn =qyn , qun =qyn , qvn =qyn , qsn =qyn , q2 un =qyn2 , q2 vn =qyn2 , q2 T n =qyn2 , qpn =qt n , qun =qt n , qvn =qt n , qsn =qt n can be computed. Along the y-boundary, all the fluxes are then prescribed as follows:     1 qpn qvn 7 Y n 7 ¼ vn 7 g c n rn cn qyn qyn

X nv ¼

Y nu ¼ vn

qun qyn

Y ns ¼ vn

qsn qyn

 qvn 1  n 1 Y Y n þ  n þ n Y nm 2 qt r

X ns ¼ Y ns 

Xn 7 ¼

qsn R þ Snm qt n pn

 1  n 1 qpn qun g1 1 Y þ Y n þ 7Y nu  n n n 7 n þ n n Snm 7 n X nm 2 r c qt qt rc r

(10)

The last three equations cannot be exactly calculated from available data stored at the left or right boundary due to the presence of viscous terms. Indeed the computation of Snm , X nm and Y nm terms is not possible (see Eqs. (6)–(8)) due to the lack of data information related to the computation of x-fluxes. Then it is only possible to provide an estimation of viscosity fluxes: Snm , X nm and Y nm using the following expressions: X nm ¼

Y

Snm ¼

kn q2 T n s0

qy2

nm

! þ

mn

¼

3

3

mn 3

U 20 mn s0 T 0 3

q2 un qyn2

!

q2 vn 4 n2 qy

4

!

 n 2  n 2 ! qv qu þ3 qy qy

(11)

Bottom and top inflow boundary condition. The treatment of inflow condition along the y-direction is similar to the previous one detailed above, apart from the order of computation of the fluxes, knowing that in this case the x-derivatives of the flow variables are now available. Numerical algorithms. The pseudo-characteristic formulation of the Navier–Stokes equations provides a decomposition of the pressure, velocity and entropy fluxes. This enables a very simple and natural use of upwind schemes. To enforce both numerical stability and accuracy for wave propagation problems, a high order upwind Dispersion Relation Preserving (DRP) scheme is used. In the interior nodes the fourth-order accurate upwind biased DRP scheme is used and this scheme is modified near the computational domain boundaries [21]. Time integration is performed using a third-order TVD Runge–Kutta scheme [23].

3660

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

3.3. Time reversal method Time reversal method consists in using the invariance of equations with respect to time and to reverse the direction of the flow: 8 t-t > > > > > r > < ðx,tÞ-rðx,tÞ vðx,tÞ-vðx,tÞ (12) > > > pðx,tÞ-pðx,tÞ > > > : sðx,tÞ-sðx,tÞ With such a change, equations remain formally unchanged if the dissipation term is omitted. If it is taken into account, the equations are not invariant under the proposed transformation and theoretically time reversal method cannot be applied. Nevertheless, due to its property of matched filter, time reversal method still holds and provides good result as will be shown in the section devoted to results. From a numerical point of view, that operation consists in flipping the stored data at the boundaries of the direct computation and then to prescribe the time reversed data as inflow boundary conditions for the second numerical simulation. 3.4. Complex differentiation method Complex differentiation first introduced by Lyness and Moler [24] is used to compute the sensitivity derivatives of a function [25]. For instance, a function f(z) of a complex variable z ¼ x þ iE with E 5 1 can be expressed thanks to a Taylor series expansion around x and the derivative is obtained with the imaginary part of the expansion, provided an approximation of E2 of the gradient. A complex differentiation method is used in parallel of the time reversal method. It is capable to tag information in the stored variables and to follow it during the reversed simulation. So, thanks to differentiation complex it is possible to follow the propagation of the acoustical waves emitted by the boundaries of the computational domain and to distinguish in the pressure field the aerodynamic and the acoustical parts, respectively. From a numerical point of view, that technique is implemented by extending physical flow variables in the complex plane in the numerical solver [9]. Then, a small imaginary part is added to acoustical disturbances measured and time reversed on the boundaries of the computational domain. The pseudo-characteristic equations are then numerically solved for both the real and the imaginary parts, simultaneously. The monitoring of the imaginary part gives information about the acoustical pressure field. 3.5. Numerical and viscous dissipation investigation In order to avoid any misinterpretation of the results, it is necessary to distinguish between dissipation due to numerical artifact and dissipation due to physical process. Here, the numerical dissipation due to viscous effects is briefly analyzed. To this end, the propagation of a plane wave in a two-dimensional lossless medium ðm ¼ 0Þ and also in a homogeneous but dissipative viscous medium ðmn a0Þ is considered. Different test-cases are performed using multiple values of the viscosity parameter mn . Special attention is paid to the kind of inflow condition. As shown in the section devoted to the boundary condition, the inflow condition has to take into account the fluid viscosity. Then, one also studies the wave propagation in viscous medium but using standard inflow condition with no dissipation process (inviscid inflow condition). That is done by examining the effect of inflow specification onto the physical dissipation. Different Cartesian grids with a varying number of grid points per wavelength are successively used to examine the wave propagation dissipation. The time-step is specified such that the CFL number remains constant. In each case, the simulation is performed during a time corresponding to a distance of propagation of more than 30 wavelengths of the initial acoustic wave. Then, the amplitude of the pressure associated with the acoustic wave after 30 wavelengths is compared to the initial reference one. As an illustration, for a fixed mesh grid which corresponds to approximately 100 points per wavelength, the resulted logarithm representation of the physical dissipation is represented as a function of mn (Fig. 1). In this figure, the dashdotted line indicates the level of numerical dissipation when dealing with Euler equations ðmn ¼ 0Þ. The dashed line indicates the results deduced from numerical simulations using the viscous inflow condition and the black line corresponds to the results for the standard inviscid inflow condition (fluxes given in Eqs. (10) with mn ¼ 0). Note that, for high values of mn ð 4 5:5  103 Þ, one has to impose the viscous inflow condition ðmn a0Þ otherwise no physical acceptable solution is obtained. It is then emphasized that there is no effect of the nature of the inflow condition on the viscous dissipation until mn r 5:5, corresponds to a physical viscosity of 10  4 Pa s  1. This underlines that when dealing with air flow, using the standard inflow condition provides similar results as the ones deduced from a numerical simulation based on the viscous inflow condition. Moreover, for air flow simulation ðmn ¼ 1Þ, the level of viscous dissipation is only slightly superior to the one of purely numerical dissipation ðmn ¼ 0Þ. Those results show that physical dissipation for standard configurations

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

3661

Fig. 1. Logarithm representation of the physical dissipation of the acoustic wave as a function of mn . Dashed line corresponds to numerical results using viscous inflow condition. Black line indicates results based on the standard inflow condition (for which mn ¼ 0). Dash-dot line indicates the level of numerical dissipation.

is dominant versus the numerical one. Consequently, the present solver can be used to study configurations with viscous flows. 4. Numerical localization of sound sources in viscous flows: results In the theoretical analysis, it was shown that three configurations have to be distinguished in order to apply time reversal method in aeroacoustics. In a previous paper [9], configurations without dissipation were considered with an emphasis on different physical mechanisms responsible for sound generation. In the present part, one focuses on the two configurations with viscous effects for steady and unsteady flows. 4.1. Time reversal method for steady flow with viscosity To investigate the possibility to use the time reversal method in order to localize sound sources in viscous steady flows, two configurations are simulated. The first one consists in localizing the position where mass injection is made in a medium at rest (no flow) for different values of viscosity. The second configuration deals again with the localization of the position where mass is injected but in a steady flow (a steady vortex) for different values of viscosity. 4.1.1. Mass injection in a medium without flow but with viscosity Mass injection is known to generate pressure waves in the general case [26]. In that case, an injection of mass is made by imposing the following initial impulsive wave at t n ¼ 0:

rn ðxn ,yn ,tn ¼ 0Þ ¼ 1 þa expðxn2 þyn2 Þ 4

(13)

with a ¼ 10 . The computational domain size ð5,5; 5,5Þ corresponds to a uniform mesh of (300  300) points. The numerical simulation (first direct problem) is performed for different values of mn from mn ¼ 0 to mn ¼ 5:5  104 . Fig. 2 shows two instantaneous fluctuating pressure fields extracted from the direct simulations with mn ¼ 0 (inviscid gas) and mn ¼ 5:5  103 . That value of mn is extremely important and does not correspond to any physical value. Again, this choice was made in order to outline the applicability of the time reversal method even for extreme configurations. During each direct simulation, instantaneous variables (un, vn, pn, sn and Tn) are stored at four boundaries and at the end of simulation, once the initial impulsive pulse is entirely evacuated from the computational domain, all variables are stored in the whole computational domain. That database is used as initial conditions for the time-reversed simulation. In that configuration there is no flow. Consequently, the reverse flow reciprocity method does not need to be applied in order to recover reciprocity. Hence, only the robustness of time reversal versus dissipation of energy is tested. Test cases associated with mn ¼ 0 and with mn ¼ 5:5  103 are presented in Fig. 2 for times corresponding to the direct simulation. The same colormap is used on each figure. From a global point of view, it is observed that the initial pulse shape is effectively preserved. Also the location of the initial acoustic source is well retrieved in each test case through time reversal. Time reversed fluctuating pressure distribution obtained at the end of the computation is plotted along yn ¼ 0 (see Fig. 3). On this figure, the pressure field deduced from the time-reversed Euler equations is superimposed onto the ones computed from the Navier–Stokes equations using different values of mn . Moreover, the initial reference acoustic pressure

3662

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

Fig. 2. (a) Direct sense mn ¼ 0 and (b) Direct sense mn ¼ 5.5 103: direct simulation. Fluctuating pressure field at two instants. Top: initial pulse shape ðt n ¼ 0Þ. Bottom: instantaneous fluctuating pressure field at a further instant. (c) Reversed Simul. mn ¼ 0 and (d) Reversed Simul. mn ¼ 5.5 103: time reversed simulation. Fluctuating pressure field at similar time-reversed instants. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

3663

Fig. 3. Zoomed view of instantaneous pressure distribution along yn ¼ 0 obtained at the end of the time-reversed simulation.

Fig. 4. Logarithm representation of the relative pressure amplitude error computed from the difference between the time-reversed pressure maximum (using different mn values) and time-reversed pressure maximum deduced from Euler simulation ðmn ¼ 0Þ.

field generated in the first direct simulation at t ¼0 is also plotted for comparison. As a first result, we do not observe any difference between viscous results with mn inferior to 55 and results dealing with Euler equations ðmn ¼ 0Þ. This confirms that for small mn values, the viscous dissipation is of an order of 104 (see Fig. 1) which cannot be visible on such a representation. In these cases, according to the spatial discretization ðDxn ¼ 0:033Þ, the slight amplitude decay obtained in the reversed simulation has its main origin in the numerical dissipation, of an order of 102 . Conversely, when mn is superior to 552, the pulse pressure amplitude decreases as the fluid viscosity increases. Fig. 4 presents the logarithm representation of the normalized maximum amplitude error as a function of mn . For a particular test case with mn , the maximum amplitude error is the difference between the current maximum amplitude and the one obtained with mn ¼ 0. The error value is normalized with reference pressure value p0. In this configuration, the error evolves linearly with the normalized viscosity mn . Indeed, the high dissipative mechanism due to viscosity occurs for high mn values. The current case shows that the time reversal method is able to localize sources both in space and in time even if levels of amplitude are not well retrieved because of dissipation effects. This result is similar to previous studies [15,27]. It results

3664

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

Fig. 5. (a) Direct sense mn ¼ 0 and (b) Direct sense mn ¼ 552: instantaneous fluctuating pressure fields deduced from the interaction of an initial pressure pulse and a steady vortex. Top to bottom: it ¼0, it1 and it2. (c) Reversed Simul. mn ¼ 0 and (d) Reversed Simul. mn ¼ 552: instantaneous fluctuating pressure fields deduced from the reversed flow simulation. Top to bottom: N it it 2 , N it it 1 and Nit.

from the spatio-temporal matched filter property [19]. The next paragraph shows that this result is still valid in the presence of flows. 4.1.2. Mass injection in a medium with flow and viscosity The second test case deals with the interaction between a mass injection and a viscous vortex. As a first step, a steady vortex is generated using Taylor’s model with an exponentially decaying velocity profile leading to a zero circulation vortex. This vortex located in the center of the computational domain has a core radius r0 of 5 and a Mach number, M of 0.1. Lengths are normalized

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

3665

with the core radius of the vortex. The streamwise and transverse velocity components are specified as follows:   1r n2 un ¼ Myn exp 2

vn ¼ Mxn exp



1r n2 2



with r n2 ¼ xn2 þyn2 . The computational domain size is ð40,40; 40,40Þ corresponding to a uniform mesh of (800  800) points. Once this vortex has reached a steady state, an injection of mass is made at t n ¼ 0 and located at ðxn ¼ 1:5,yn ¼ 0Þ (see the top graph of the first two columns presented in Fig. 5). The shape of the initial pulse is the same as the one given with Eq. (13). The mass injection induces a wave that propagates and interacts with the steady vortex. Note that the amplitude of the associated acoustic pulse is very small compared to the one arising from the vortex. Again time reversal method is applied for different value of the viscosity parameter. Fig. 5 displays three instantaneous fluctuating pressure fields corresponding to time iteration: 0, it1 and it2, respectively. In this figure, both test cases are represented using inviscid fluid ðmn ¼ 0Þ and viscous fluid ðmn ¼ 552Þ, respectively. The scattering effect due to the presence of the steady vortex is clearly visible on these figures and especially for iteration it2. Note that the effect of viscous dissipation is not visible on these representations. During the flow simulations, at each time step, flow variables are saved at the four boundaries. These simulations are performed for Nit time iterations. This time duration is sufficient for the initial pulse pressure to exit the computational domain. At that last time step, the whole set of flow variables are stored in the whole computational domain. After timereversing stored data, they are used as inflow and initial conditions, respectively, for the second time-reversed simulations. Fig. 5 represents three instantaneous fluctuating pressure fields obtained at Nit it2 , N it it1 and Nit time iterations. These

Fig. 6. Zoomed view of the reversed pulse. Instantaneous pressure distribution along y¼ 0 deduced from reversed simulations using different mn values.

Fig. 7. Pulse amplitude error as a function of mn values. Logarithm representation.

3666

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

time iterations are the same time-reversed instants as the ones presented during the direct simulation (first two columns). When dealing with a fluid having an important dissipation ðmn ¼ 552Þ, the shape and the location of the initial pulse are well recovered. Again, the main difference concerns the decay of the amplitude of the initial pulse. Fig. 6 displays a zoomed view of the fluctuating pressure distribution along yn ¼ 0 and around the initial pulse location, xn 2 ½1:8,1:3. That figure shows that the maximum of amplitude, which corresponds to the localization of the source, is well located (exactly in yn ¼ 0) and comes exactly at the correct time. The only difference concerns the decay of the amplitude of the reversed acoustic pressure field according to the dissipation effects. Fig. 7 displays the evolution of the amplitude error of the acoustic pulse as a function of mn . The computation of this error is similar to the one detailed in the previous section. Due to the presence of high intensity vortex, the amplitude error is higher than in the previous test case (see Fig. 4) and does not vary linearly as a function of mn . Indeed, in the previous case all the portions of the wavefront experimented the same attenuation as their acoustical path was the same. In the current case, Fig. 5 shows that the wavefront is strongly distorted and so the paths followed by the different portions of the wavefront do not have the same length. Consequently, their amplitude more or less decreases during the propagation in the direct and the reversed propagation. The result is the nonlinear variation of the error with parameter m. Nevertheless, as in the previous case, only levels are affected. Position and time of emission are well retrieved.

4.2. Noise generated by flow instability: configuration with unsteady viscous flows In the first part, devoted to theoretical analysis, it was shown that time reversal method cannot be applied to configurations with unsteady flows with or without viscosity. Such situations may occur for source generation by the flow itself (in jet noise for instance) or for sound generated by interaction of the flow with solid boundaries. As previously mentioned the latter case is beyond the scope of this paper even if it concerns interesting applications. Concerning the generation of sound by the flow itself the method has been successfully applied to localize sources of sound in a plane mixing layer [9]. In the present paper, the sound generated by a plane mixing layer is investigated again in order (i) to explain why it works even if theoretical assumptions are not respected (ii) to extend previous results to cases with viscosity. The numerical simulation of a two-dimensional forced plane mixing layer is then considered. Such a flow configuration allows the implementation and the validation of the current methodology because aeroacoustic sources are well located in regions where pairings of vortices occur [28] . The computational domain has a length equal to 220do0 in the stream-wise direction (x-direction) and a width equal to 252do0 in the transverse direction (y-axis), where do0 is the inlet vorticity thickness used to normalize length scales. A Cartesian grid of 1800  1800 points is used. A sponge zone of 225 grid points along the x-direction is used to prevent spurious reflections.

Fig. 8. Instantaneous vorticity field. Top to bottom: (a) mn ¼ 0; (b) mn ¼ 5:5; (c) mn ¼ 55 and (d) mn ¼ 552.

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

3667

Inflow disturbances are generated at the left boundary ðxn ¼ 0Þ according to the configuration studied by Colonius et al. [29]. These perturbations are composed of a combination of linear eigenfunctions obtained from the linear stability analysis of the plane mixing layer. These eigenfunctions can be obtained from the Rayleigh equations [29] or from the linear viscous theory [30]. In the present work, inflow conditions correspond to the fundamental frequency f0 and the first three subharmonics of the eigenfunctions computed as in [29]. Even if there exist some small differences in the resulting eigenfunctions [30], we expect that for the current study the inflow used for the spatial forcing has very weak effects on the aeroacoustic source localization results. Configurations for various values of the viscosity parameter (mn : 0; 0.55; 5.5; 55; 552) are computed. Results deduced from computation with viscous fluid, mn ¼ 0:55 are quite similar to the ones with mn ¼ 0. For each test case, 20 000 time steps are computed once the flow is fully developed (the initial transitory flow is finished), Fig. 8 displays the instantaneous vorticity field of selected test cases at the same time. The x- and y- domains are limited to better appreciate the effects of fluid viscosity on the first developing flow scales. The first pairing process occurs at a similar location, for computations using mn ¼ 0, mn ¼ 5:55 and mn ¼ 55, but the spatial development of instantaneous larger scale structures is quite different in each case. Furthermore, when using a high viscosity value, mn ¼ 552, the flow field is very different from the other ones even for smaller flow scales. This is directly related to the high value of the viscosity which leads to some important energy dissipation. During the flow simulation, instantaneous flow variables are stored at the top and bottom boundaries as well as at the end of the physical computational domain. Fig. 9 provides spectra of the acoustical pressure signals received at the top boundary ðxn ¼ 100Þ for different values of the parameter of viscosity. To enhance the computation of the spectra, instantaneous velocity signals are cut up into several blocks with 75 percent overlapping and an Hann window is applied to each block. These spectra mainly exhibit frequency peaks at subharmonics of the fundamental frequency f0. The time reversal method is then used to localize the origin of the acoustic waves at frequency f 0 =2 detected at the top and bottom boundaries. Flow variables stored at the outflow boundaries during the first direct computation are time and flow reversed. Then, the pressure field filtered at the frequency f 0 =2 is added in the imaginary part of the complex pressure field to tag the chosen acoustic wave associated

Fig. 9. Pressure spectra computed at the top boundary computational domain ðX n ¼ 100Þ. (a) mn ¼ 0; (b) mn ¼ 5:5; (c) mn ¼ 55 and (d) mn ¼ 552.

3668

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

Fig. 10. Instantaneous real vorticity field superimposed onto the fluctuating pressure imaginary part filtered to f 0 =2 frequency. (a) mn ¼ 0; (b) mn ¼ 5:5; (c) mn ¼ 55 and (d) mn ¼ 552.

with f 0 =2 frequency. The imaginary part of the other flow variables is put to zero value. These complex reversed flow variables are finally used as inflow conditions for the second time-reversed simulations. Fig. 10 shows at the same time the instantaneous imaginary part of the pressure field which is superimposed onto the corresponding real vorticity field. Note that the real vorticity field is obtained from the preliminary direct simulation.The first test case using mn ¼ 0 allows to recover previous results [9]: the imaginary part of the pressure field converges towards the origin of the acoustic wave of f 0 =2. This origin coincides with the first pairing observed in the instantaneous vorticity field known to be the region where sound is generated [28]. That result shows that the time reversal method works to localize sources generated in a mixing layer. But because of the unsteadiness of the flow, it has been theoretically predicted that time reversal method cannot be used. Nevertheless, in that configuration, the region where flow is unsteady is confined in the zone of high vorticity. The extent of this zone is measured by the vorticity thickness (here do  1). In comparison, the wavelength associated to acoustical waves (characteristic length for acoustic process) is higher. For the highest frequency considered (and so the smallest acoustical wavelength) which is f 0 =2, the wavelength is about equal to 25. That difference between the characteristic scales of aerodynamic and acoustical processes explains why time reversal method seems to be working when the interest is focused on acoustical variables. Time reversed acoustical waves are only slightly altered by the no-reciprocity which is essentially localized in a very small region of space compared to the acoustical wavelength. Numerical simulations based on viscous fluid with mn ¼ 5:5 and with mn ¼ 55 permit to also recover the origin of the acoustic wave even if the level of the imaginary part is smaller than the one computed during the reference simulation ðmn ¼ 0Þ. It is then demonstrated that the time-reversal procedure coupled to the complex differentiation technique is effective in localizing the aeroacoustic sources in unsteady viscous flows while acoustical lengths are high compared to aerodynamic ones. An interesting feature of this methodology is observed for the simulation based on high viscous fluid, mn ¼ 552. Indeed, the instantaneous imaginary part of the pressure field shows that acoustic wave detected at the computational boundary domains has two origins: the first pairing process and the inflow conditions. This is due to the inflow condition imposed during the first direct simulation contains also a wave of f 0 =2 frequency. Then knowing that the amplitude of the aeroacoustic source wave is very small due to high energy dissipation, its amplitude may be of an order of the one imposed at inlet. The implementation of the sensitivity analysis allows the recovery of both acoustic sources: an artificial one (due to forced inflow condition imposed at the first direct simulation) and a more realistic one arising from vortex pairing process. It is noticeable that the pairing area is shifted compared to the other case (see Fig. 8). This shift is recovered by the method since the waves – except those ones which retropropagate to the origin – go back to the same location. Previous theoretical time-reversal method analysis is confirmed. Indeed, the time-reversed simulation using specific inflow conditions (with large associated wavelengths) cannot restitute the aerodynamic small flow scales. However, solving the time-reversed equations associated with the imaginary part of the variables permits to follow a particular acoustic wave. More precisely, the back-propagation of such a wave is identified and its aeroacoustic origin is located thanks to the acoustic wavelength resolution. Furthermore, adding viscous effects in the reversed simulation acts as an energy filtering process without any modification onto the aeroacoustic source localization method.

P. Druault et al. / Journal of Sound and Vibration 332 (2013) 3655–3669

3669

5. Conclusion Successive flow configurations including steady (with viscous fluid) flows and unsteady ones (with dissipation) have been performed to test and validate the time-reversal procedure method for aeroacoustic source localization. According to the steady flow configuration results, it is demonstrated that the time reversal procedure allows to detect the origin of the aeroacoustic source but the source amplitude is reduced due to the viscous energy loss. It is then confirmed that the application of the methodology permits to recover the shape and the location of the source of sound. By coupling this time-reversal procedure with the complex differentiation technique, it is also demonstrated that it is possible to detect and localize a particular acoustic wave generated in a turbulent flow. The presence of fluid viscosity only acts as an energy filtering without any modification of the aeroacoustic source localization methodology. Furthermore, it is emphasized that even if the presence of an unsteady flow (including viscous fluid) theoretically breaks the reversibility of the fluid motion equations, such an approach is robust and effective for the detection of particular aeroacoustic sources. This methodology proposes a new way for aeroacoustic source identification method without any a priori of the source. Furthermore, further analysis of the imaginary part associated with the back-propagation of the sound wave confirms the nature (quadripolar-type) of the aeroacoustic source in a plane mixing layer flow [31]. The present study is focused on sound generation of free shear flows. Further investigations are required to extend the present method to the case of turbulent wall-bounded flows, which raise new theoretical and practical problems dealing with breakdown of reciprocity. References [1] M. Lighthill, On sound generated aerodynamically: I. General theory, Proceedings of the Royal Society of London 211 (1952) 564–587. [2] J.E. Ffwocs Williams, D. Hawkings, Sound generation by turbulence an surfaces in arbitrary motion, Proceedings of the Royal Society of London A 264 (1969) 321–342. [3] N. Curle, The influence of solid boundaries upon aerodynamic sound, Proceedings of the Royal Society 231A (1955) 505–514. [4] M. Howe, Acoustics of Fluid–Structures Interactions, Cambridge University Press, 1998. [5] P. Druault, M. Yu, P. Sagaut, Quadratic stochastic estimation of far field acoustic pressure with coherent structure events in a 2D compressible plane mixing layer, International Journal for Numerical Methods in Fluids 62 (2010) 906–926. [6] P. Druault, X. Gloerfelt, T. Mervant, Investigation of flow structures involved in sound generation by two- and three-dimensional cavity flows, Computers and Fluids 48 (1) (2011) 54–67. [7] A. Jameson, L. Martinelli, N. Pierce, Optimum aerodynamic design using the Navier–Stokes equations, Journal of Theoretical and Computational Fluid Dynamics 10 (1998) 213–237. [8] M. Giles, N. Pierce, Analytic adjoint solutions for the quasi-one-dimensional euler equations, Journal of Fluid Mechanics 426 (2001) 327–345. [9] A. Deneuve, P. Druault, R. Marchiano, P. Sagaut, A coupled time-reversal complex differentiation method for aeroacoustic sensitivity analysis: towards a source detection procedure, Journal of Fluid Mechanics 642 (2010) 181–212. [10] M. Fink, D. Cassereau, A. Derode, C. Prada, P. Roux, M. Tanter, J.-L. Thomas, F. Wu, Time-reversed acoustics, Progress Report in Physics 63 (2000) 1933–1995. [11] D. Cassereau, M. Fink, Time-reversal of ultrasonic fields. iii theory of the closed time-reversal cavity, IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control 39 (5) (1992) 579–592. [12] A. Tarantola, Theoretical background for the inversion of seismic waveforms, including elasticity and attenuation, Pure and Applied Geophysics 128 (1988) 365–399. [13] J.-L. Thomas, Etude des mirroirs a retournement temporel. Applications a la focalisation des ondes ultrasonores et a la lithotritie. (Study about Time Reversal Mirrors. Applications to Ultrasonic Waves Focusing and Lithotritie), PhD Thesis, Univerisite´ Pierre et Marie Curie - Paris 06, 1994. [14] R. Carminati, J. Sa enz, J.-J. Greffet, M. Nieto-Vesperinas, Reciprocity, unitarity, and time-reversal symmetry of the S matrix of fields containing evanescent components, Physical Review A 62 (2000). 012712-1–7. [15] J.L. Thomas, Ultrasonic beam focusing through tissue inhomogeneities with a time reversal mirror: application to transskull therapy, IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control 43 (6) (1996) 1122–1129. [16] E. Vergnault, O. Malaspinas, P. Sagaut, A time-reversal lattice Boltzmann method, Journal of Computational Physics 230 (2011) 8155–8167. [17] P. Roux, M. Fink, Experimental evidence in acoustics of the violation of time-reversal invariance induced by vorticity, Europhysics Letters 32 (1995) 25–29. [18] P. Roux, J. de Rosny, M. Tanter, M. Fink, The Aharonov–Bohm effect revisited by an acoustic time reversal mirror, Physical Review Letters 79 (1997) 3170–3173. [19] M. Tanter, J.-L. Thomas, M. Fink, Time reversal and the inverse filter, Journal of the Acoustical Society of America 108 (2000) 223–234. [20] S. Lu, P. Sagaut, Direct sensitivity analysis for smooth unsteady compressible flows using complex differentiation, International Journal for Numerical Methods in Fluids 53 (2007) 1863–1886. [21] S. Lu, P. Sagaut, Pseudo-characteristic formulation and dynamic boundary conditions for computational aeroacoustics, International Journal for Numerical Methods in Fluids 53 (2007) 201–227. [22] J. Sesterhenn, A characteristic-type formulation of the Navier-Stokes equations for high order upwind schemes, Computers & Fluids 30 (2001) 37–67. [23] C. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes, Journal of Computational Physics 83 (1989) 32–78. [24] J. Lyness, C. Moler, Numerical differentiation of analytic functions, SIAM Journal on Numerical Analysis 4 (1967) 202–210. [25] W. Squire, G. Trapp, Using complex variables to estimate derivatives of real functions, SIAM Reviews 10 (1) (1998) 110–112. [26] L. Kovasznay, Turbulence in supersonic flow, Journal of the Aeronautical Sciences 20 (10) (1953) 657–682. [27] J. Garnier, A. Nachbin, Eddy viscosity for time reversing waves in a dissipative environment, Physical Review Letters 93 (15) (2004) 1–4. [28] C. Bogey, C. Bailly, D. Juve, Numerical simulation of sound generated by vortex pairing in a mixing layer, AIAA Journal 40 (2) (2000) 235–243. [29] T. Colonius, S. Lele, P. Moin, Sound generation in a mixing layer, Journal of Fluid Mechanics 330 (1997) 375–409. [30] A. Babucke, M. Kloker, U. Rist, DNS of a plane mixing layer for the investigation of sound generation mechanisms, Computer and Fluids 37 (4) (2008) 360–368. [31] P. Druault, R. Marchiano, P. Sagaut, Time reversal method coupled to complex differentiation technique for the aeroacoustic source detection in viscous flow, 18th AIAA/CEAS Aeroacoustic Conference, AIAA 2012-2285, 2012.