Aero and vibroacoustical prediction of the noise generated by turbulent boundary layers

(1)

KTH Royal Institute of Technology

Department of Vehicle Engineering

Aero and vibroacoustical prediction of

the noise generated by turbulent

boundary layers

Alberto ALONSO PINAR

supervised by

Romain Leneveu and

Dr. Mats ˚_Abom

February 2019

Submitted in part fulfilment of the requirements for the degree of Master in

Aerospace Engineering of KTH Royal Institute of Technology

(2)

(3)

Abstract

Flow noise is a major noise source in the transportation industry from the automotive up to the aircraft industry. The characterization of the aerodynamic excitations and of the structure transmission is of primary interest in order to improve the passenger and crew comfort. One important component of the flow noise excitation comes from the wall pressure fluctuations induced by the development of turbulent boundary layers around the vehicle.

This master thesis is focused on the validation of a method to calculate structural vibrations based on aerodynamic data by the comparison with experimental data. The analysis will be done with a flat plate configuration in a wind tunnel. The strategy yields good results and it has been validated: it could be used for more complex problems.

(4)

(5)

iv CONTENTS 2.5 Auto spectra . . . 15 2.5.1 Amiet . . . 15 2.5.2 Efimtsov . . . 15 2.5.3 Chase - Howe . . . 16 2.5.4 Goody . . . 16 2.6 Cross spectra . . . 17 2.6.1 Corcos . . . 18 2.7 Vibrations . . . 18 3 Flat plate 21 3.1 Strategy . . . 21 3.2 Configuration . . . 22 3.3 CFD Simulation . . . 23 3.3.1 Mesh . . . 24 3.3.2 Computation . . . 28 3.3.3 Post-Processing . . . 29 3.4 FEM calculation . . . 32 3.4.1 Setup . . . 32 3.4.2 Calculation . . . 34 3.5 Results . . . 37

3.5.1 CFD and experimental validation . . . 37

(7)

3.5.3 Microphone array results . . . 48 3.5.4 FEM results . . . 52 4 Conclusion 55 4.1 Applications . . . 56 Bibliography 56 v

(8)

(9)

List of Tables

3.1 Comparison of correlations formulas . . . 27

3.2 Modal base comparison . . . 33

3.3 Point 19 at U0 = 30m.s−1 . . . 40 3.4 Point 22 at U0 = 30m.s−1 . . . 40 3.5 Point 19 at U0 = 50m.s−1 . . . 40 3.6 Point 22 at U0 = 50m.s−1 . . . 40 3.7 Point 19 at U0 = 60m.s−1 . . . 41 3.8 Point 22 at U0 = 60m.s−1 . . . 41 vii

(10)

(11)

List of Figures

1.1 Vibratec’s logo . . . 1

1.2 Schematic describing the strategy and the physical phenomenons behind it. . . . 5

2.1 Schematic describing the boundary layer. Source: wikipedia . . . 8

2.2 Schematic describing the flow over a flat plate. Source: wikipedia . . . 9

2.3 u+ vs y+ graph showing the different sub layers. Source: wikipedia . . . 11

2.4 Vibrating flat plate excited by a TBL. U∞ is the outer fluid velocity and δ is the

boundary layer thickness. Extracted from [1]. . . 14

2.5 Wavenumber frequency spectra at a given frequency. See [2] for more details. . . 17

2.6 Comparison of cross spectra values with experimental and numerical results. The

experimental data is in traingles, dots and cross. More details on [3]. . . 19

3.1 Experimental setup at LMFA. Upper left is the wind tunnel in the anechoic

chamber, upper right is the instrumented flat plate and lower is a schematic of

the wind tunnel with the main dimensions. The test section if where the flat

plate is placed. . . 23

3.2 Illustration of a stl geometry (source: cfd.direct ). . . 24

3.3 Creation of the background mesh (source: cfd.direct ). . . 25

(12)

x LIST OF FIGURES

3.4 Refinement of the mesh (source: cfd.direct ). . . 25

3.5 Elimination of inside cells (source: cfd.direct ). . . 25

3.6 Surface snapping and layer addition (source: cfd.direct ). . . 25

3.7 Illustration of boundary layer solving strategies (source: Leap Australia). . . . 26

3.8 The geometry that was used for this study. One should note the inlet offset. . . 27

3.9 Zoom on the near wall zone of the mesh. . . 28

3.10 Residuals for U0 = 30m.s−1 and y+ = 30 . . . 30

3.11 FEM mesh model of the flat plate . . . 32

3.12 Main modes of the flat plate with their corresponding frequency at their right. . 34

3.13 Cross correlation matrix for U0 = 30m.s−1 and f = 355Hz. The axis represent the FEM node ID. . . 35

3.14 Edat file showing the implementation of the boundary layer parameters. . . 36

3.15 Edat file showing the wall pressure spectra construction for Actran. . . 36

3.16 Point 19 at U0 = 30m.s−1 . . . 38 3.17 Point 22 at U0 = 30m.s−1. . . 38 3.18 Point 19 at U0 = 50m.s−1 . . . 39 3.19 Point 22 at U0 = 50m.s−1. . . 39 3.20 Point 19 at U0 = 60m.s−1 . . . 41 3.21 Point 19 at U0 = 60m.s−1 . . . 42 3.22 Pressure field at U0 = 30m.s−1 . . . 43 3.23 Pressure field at U0 = 50m.s−1 . . . 43

(13)

LIST OF FIGURES xi

3.24 Pressure field at U0 = 60m.s−1 . . . 44

3.25 Axial velocity field at U0 = 30m.s−1 . . . 44

3.28 Friction velocity along x. . . 45

3.29 Boundary layer thickness along x. . . 46

3.30 Momentum thickness along x. . . 46

3.31 Displacement thickness along x. . . 47

3.32 Cross spectra matrix, 100Hz. . . 49

3.36 Wavenumber-frequency diagram for U0 = 10m.s−1 at left and U0 = 20m.s−1 at right . . . 50

3.37 Wavenumber-frequency diagram for U0 = 30m.s−1 at left and U0 = 40m.s−1 at right . . . 50

3.38 Wavenumber-frequency diagram for U0 = 30m.s−1: heterogeneous TBL at left and homogeneous TBL at right . . . 51

3.39 Wavenumber-frequency diagram for U0 = 60m.s−1: heterogeneous TBL at left and homogeneous TBL at right . . . 52

3.40 Corcos α vs frequency . . . 52

(14)

3.42 Acceleration at the center of the plate for different realizations. The curves are

smoothed and presented in a 1/3 octave decomposition. The vertical step is

equal to 1dB. . . 53

3.43 Acceleration at the center of the plate for different outer velocities. The vertical

step is equal to 1dB. . . 53

(15)

Chapter 1 Introduction

1.1 VibraTec

VibraTec is the main company of the group VibraGroup, the french leader in the acoustics and

structural dynamics industry. Two other companies are included in the group : Vibrateam

which provides engineering services to other companies and MicrodB which develops and

man-ufactures antennas for acoustic imaging and software for source localization.

Figure 1.1: Vibratec’s logo

The headquarter is in Ecully (a commune in the metropolis of Lyon) and has two international

branches: one in Malaysia and the other in Germany. Its sales revenue is around 10 million

euros and 20% of it is invested in research and development. The group has around a hundred

(16)

2 Chapter 1. Introduction

employees and around 60 of them are permanently based in the Ecully’s site which holds a

laboratory and a calculation office.

Vibratec offers advice for very technical and complex problems to a variety of companies in

different sectors. Their services usually include preliminary drafts giving recommendations in

order to optimize a product from an acoustical or vibratory point of view or advice on more

advanced projects where malfunctions have appeared. Finally, the company is highly engaged

in applied research in order to better understand the problems and to propose new services:

it usually has two or three PhD students working at the company. These projects are usually

presented at international conferences.

Vibratec is organised in four different business units: Railway, Automotive, Aerospace and

OilGas. The production is ensured by a transverse operational unit.

Each of these business units (BU) has a number of engineers and technicians specialized in

experimental tests or simulations. The ability to couple measurements with simulation

tech-niques is one of the main skills of the company making it especially strong when dealing with

results for the calibration or validation of a numerical model.

During my master thesis project I was part of the Aerospace BU. My supervisor, Romain

Leneveu leads the Aerospace BU and coordinates, among others, the CANOBLE project for

Vibratec. I also had the chance to interact with Martin Rissmann, who made some CFD and

vibroacoustics simulations before I started the project. Finally, I also interacted with LMFA

(Laboratoire de M´ecanique des Fluides et Acoustique, Ecole Centrale de Lyon) researchers like

(17)

1.2. Context and application 3

1.2 Context and application

This master thesis project is part of the European research program CleanSky aiming at

im-proving the performance of aircrafts and to level up aeronautics research in Europe. Indeed,

the air traffic is constantly growing (the number of passengers should double in 20 years) so big

investments in research and development are well justified. The CANOBLE (CAbin Noise from

Boundary Layer Excitation) is one of the CleanSky projects and is focused on the interior noise

in an aircraft cabin due to the structural excitation by the turbulent boundary layer developing

along the fuselage. Indeed, the interior noise can heavily influence the passenger comfort and

the structural fatigue of the aircraft, this is particularly true in long range flights.

The interior noise is created by three main sources:

1. jet-propulsion noise

2. turbulent boundary layer excitation

3. internal systems such as air conditioning

DLR (German Aerospace center) researchers have measured the different noise contributions

doing in-flight measurements ([4]) and have concluded that internal systems have a minor

impact on the interior noise far away from it. The TBL noise is the major contributor to the

interior noise particularly in the back of the plane where the boundary layer is thicker and

more developed. Finally, propulsion noise is very important in the lift-off and during the cruise

flight where it has a similar level to the turbulent boundary layer noise. It is then of major

importance to correctly characterize the flow around the aircraft and the transmission of the

excitation in order to better control the interior noise.

The CANOBLE main goal is to better understand and characterize the aerodynamical noise:

numerical simulations will be developed in order to quantify the noise radiated by the structure

and by the flow and a full-scale cabin model of a commercial jet will be used to conduct

(18)

4 Chapter 1. Introduction

1.3 Strategy

The project contains an important simulation part which goal is to quantify the interior noise

generated by the structural subject to pressure fluctuations due to the turbulent nature of the

boundary layer.

The goal of this master thesis project is to deploy calculation and simulation techniques in a

multiphysics context for this project that can also be used in a more industrial way for future

projects by generating high fidelity CFD simulations, coupling between different softwares

(be-tween OpenFoam and Actran for example) and being time efficient. One of the main aspects

of the thesis was to develop RANS simulations and vibroacoustical calculations for a flat plate

study case that will be used to verify and validate the whole simulation process. Indeed, RANS

simulations are still the only industrial simulation technique that can offer good results while

demanding less time and calculation resources than LES or DNS for highly complex

geome-tries. These simulations offer the possibility to recreate the unsteady flow parameters that

would serve as the temporal excitation for the structure. However, due to their complexity and

cost (both in time and computational resources) they are still prohibitive in an industrial case

so RANS simulations will be used and coupled with semi-empirical wall pressure spectra based

on aerodynamic parameters (such as the boundary layer thickness for example).

To sum up, the simulation strategy is to calculate the flow field using OpenFoam, then to

ex-tract the interesting parameters using Python scripts. These parameters can be directly given

to Actran which has its own implementation of the wall pressure spectra or can be stored in

order to give Actran the wall pressure spectra calculated using Python scripts. Finally, Actran

calculates the structural response to the pressure excitation and using a weak coupling strategy

the interior noise can be calculated. See figure 1.2 for a more graphical explanation of the

strategy flow.

During the thesis, not only the simulation results were of interest for Vibratec but also the

knowledge management of all the CFD scripts, Actran scripts and Python scripts that were

(19)

1.4. Objectives 5

written in order to allow the use of all the scripts for future projects. Finally, two abstracts

were presented at two conferences: one for the 25th AIAA/CEAS Aeroacoustics conference (see

[5]) to be held in May at TU Delft and the other, in collaboration with the LMFA, for the 26th

International Conference on Sound and Vibration (ICSV 19, (see [6])) to be held in Montr´eal

during July 2019.

Figure 1.2: Schematic describing the strategy and the physical phenomenons behind it.

1.4 Objectives

The main objectives of this master thesis are:

• Implementation of a TBL source generator taking into account the steady aerodynamic quantities of a boundary layer from experimental data or CFD calculations

• Creation of an analysis tool of boundary layer quantities and wall excitation quantities • Validation of the strategy on a flat plate

(20)

Chapter 2 Background Theory

This chapter will concentrate on giving a short description of turbulent boundary layers. The

main characteristics and theoretical aspects necessary to the comprehension of the thesis will

be explained.

2.1 Modelling a flow

The equations that best describe the flow of a fluid are named the Navier-Stokes equations.

There are three partial differential equations: the mass conservation (equation 2.1, involving

the density ρ, and the velocity vector V), the momentum conservation (equation 2.2, involving

among other quantities the volumic forces f, the shear stress τ and the pressure p) and the

energy conservation (equation 2.3, involving among other quantities the internal energy e, the

heat flux due to convection q and the heat flux due to radiation qr).

Dρ

Dt + ρ · div(V) = 0 (2.1)

ρDV

Dt = ρf + div(τ ) − grad(p) (2.2)

(21)

2.2. Turbulent boundary layer 7

ρD Dt(e +

V2

2 ) = ρf · V − div(pV) + div(τ · V) + div(q) + div(qr) (2.3)

The momentum conservation expresses the fact that the change of momentum is the result of

the forces applied to, or by, the fluid. The forces can be divided into three different terms: the

pressure forces, the body forces (such as gravity) and finally the viscosity forces. The last type

of forces is particularly interesting when modelling a fluid flow.

2.2 Turbulent boundary layer

2.2.1 Boundary Layer

An external flow around a certain object will exhibit many different aerodynamic phenomena.

The behaviour of the flow field will depend on a great number of parameters such as the

object geometry (its size and orientation), the fluid properties (the viscosity or the density)

and of course the freestream conditions (velocity and pressure). One of the most important

parameters that will give a rough description of the flow field is the Reynolds number defined

by Re = ρV L_µ where L is the characteristic length of the flow, ρ is the density of the fluid, V is the characteristic velocity and µ is the dynamic viscosity. It is the ratio of inertial to

viscous forces and will then predict if the flow is laminar (when the viscosity is dominant over

the inertial forces) or turbulent (when the inertial forces dominate over the viscous forces). An

aircraft cruise situation is an example of a high Reynolds flow. The flow can usually be divided

into two main regions:

1. A viscous boundary layer close to the surfaces of the aircraft

2. An almost inviscid flow outside this boundary layer

The boundary layer is the result of the no-slip condition that states that the velocity of a

(22)

8 Chapter 2. Background Theory

velocity of the fluid next to it is then zero. As a result of this condition and of friction forces

(due to viscosity) a velocity gradient region is created close to the walls known as boundary

layer. The force exerted by the closest fluid particles to the wall is known as the wall shear

stress and is defined by:

τw = µ ∂u ∂y y=0 (2.4)

Where y is the normal distance to the wall and u the flow velocity. For a Newtonian fluid, the

shear tensor is proportional to the velocity gradient: τ = µ∇u.

One of the main characteristics of the boundary layer is its thickness, usually denoted δ in

the literature, that can be assumed to be a function of the flow principal direction. As the

transition from boundary layer to the outer flow is continuous, a precise thickness value of the

boundary layer does not make any sense. The thickness is often and arbitrarily given as the

point where the velocity reaches 99% of the outer velocity (see figure 2.1).

Figure 2.1: Schematic describing the boundary layer. Source: wikipedia

For simplicity reasons, one can consider a flat plate subject to a flow as in figure 2.2 . In this

case, the boundary layer will start at x = 0 as a region where fluid particles move smoothly

one along the others. This is the laminar boundary layer. Due to the increase in thickness and

the increase of the axial Reynolds number the flow becomes unstable at the position xcr (see

figure 2.2) which corresponds to a Reynolds number of around Re = 5 · 105_{. At some point,}

(23)

2.2. Turbulent boundary layer 9

turbulent boundary layer at around a Reynolds number of Re = 3 · 106. This transition is random, unsteady an 3D-dependent making it very difficult to characterize (see [7]). Finally,

in the turbulent boundary layer the fluid moves in more complicated paths. The turbulent

boundary layer is of main interest for vibration analysis due to the pressure fluctuations as it

will be discussed later.

Figure 2.2: Schematic describing the flow over a flat plate. Source: wikipedia

Characterizing a boundary layer is then quite difficult. Two main thickness are usually

calcu-lated to describe a boundary layer:

1. The displacement thickness δ1 =

R∞

0 (1 − u(x)

U0 ) dx is the distance by which a surface would

have to be moved in the direction perpendicular to its normal vector away from the

reference plane in an inviscid fluid flow to give the same flow rate as occurs between the

surface and the reference plane in a real fluid. (see [8])

2. The momentum thickness δ2 =

R∞

0 u(x)

U0 (1 −

u(x)

U0 ) dx is the distance by which a surface

would have to be moved parallel to itself towards the reference plane in an inviscid fluid

stream of velocity to give the same total momentum as exists between the surface and

(24)

2.2.2 Turbulent boundary layer description

In a turbulent boundary layer, the fluid motion is chaotic creating as a result pressure

fluctu-ations both in space and time domain. It can be divided into different regions (see 2.3 for a

graphical visualization).

The closest region to the wall is called the viscous sublayer. After this, there is the buffer layer

and finally the log layer. The main parameters that influence the behaviour of the flow in a

turbulent boundary layer are: the fluid density, the kinematic viscosity, the distance to the wall

and the wall shear stress.

A normalized (or non-dimensional) velocity and distance to the wall can be defined as:

u+= u uτ

(2.5)

y+= yuτ

ν (2.6)

Here, uτ is the friction velocity defined by uτ = pτw/ρ. Each sublayer is defined by its

behaviour.

In the viscous sublayer, the viscosity dominates over the inertia forces due to the fact that the

velocity is zero at the wall: as a result, µ∂_∂22u_y = 0 so u = Cy. This is the region up to around

y+ _{= 5. The velocity profile is then linear and is given by: u}+_{= y}+_.

The buffer layer is where the inertia forces start to dominate over the viscosity. This is the

region from y+= 5 up to y+ = 20 − 30. In this layer no law holds entirely. It is then of primary importance to correctly resolve the boundary layer in order to have an accurate description of

the flow near the walls.

Finally, the log layer goes from y+ _{= 30 up to y}+ _{= 300. Indeed, in this region the velocity}

profile is given by: u+ ₌ 1 κlog(y

+_{) + B. Where κ = 0.41 is the Von Karman constant and}

B = 5.1 is a general constant.

Analytical formulations that blend the viscous sublayer linear behaviour with the log behaviour

(25)

2.3. Turbulence modelling 11

Musker formulation and is given by:

y+ = u++ e−κB eκu+ − 1 −(κu +₎2 2 − (κu+₎3 6 (2.7)

Figure 2.3: u+ _{vs y}+ _{graph showing the different sub layers. Source: wikipedia}

These kind of laws are well verified in zero pressure gradient flows without separation. However,

for complex, unsteady, 3D flows it is common to find zones with pressure gradients that can

lead to recirculation regions in the boundary layer.

2.3 Turbulence modelling

In order to keep all calculations affordable for industrial applications, the CFD calculations will

use the Reynolds Average Navier Stokes (RANS) calculations. The idea behind these equations

is to assume that an instantaneous velocity can be decomposed into a fluctuating term and a

mean or time-averaged term. This is expressed as: u(x, t) = u(x) + u0(x, t).

(26)

12 Chapter 2. Background Theory cartesian coordinates: ∂ui ∂t + uj ∂ui ∂xj = fi− 1 ρ ∂p ∂xi + ν ∂ 2_u i ∂xi∂xj − ∂ 2_u iuj ∂xj (2.8)

The main problem of this expression is that the term ρuiuj, known as Reynolds stress and

named τi,j , can’t be explicitly calculated: the problem is not closed. In order to calculate

it, one option is to model it with a transport equation. It can be shown that the transport

equation includes terms with high order terms and depends on the pressure. The Boussinesq

approximation proposes a way to model the Reynolds stress:

τij = νt ∂ui ∂xj + ∂uj ∂xi −2 3ρkδij (2.9)

As a result, the Reynolds stress is defined by a turbulent viscosity νt, the trace-less mean strain

rate sensor Si,j and the turbulent kinetic energy k. The idea is then to calculate the turbulent

viscosity νt using transport equations. For external aerodynamics, the most used models are

the k − , k − ω and Spalart Allmaras.

The model that was used in this thesis is the k − ω SST which stands for Shear Stress Transport

and is a blend of the k − and k − ω models: near the wall, the model uses the k − ω which

is highly adapted for detached flows and far away from the walls, the model used it the k − .

The implementation of this model is available in OpenFoam can be used for both completely

resolved boundary layers (mesh with y+_{< 1) and with wall functions (y}+ _{> 30 ).}

2.4 Modelling the excitation

The interior noise of an aircraft is mainly due to the aerodynamical noise created by the

turbulent boundary layer developing on the aircraft. The TBL excitation can be divided into

two main sources:

(27)

2.4. Modelling the excitation 13

2. The hydrodynamic contribution related to the flow and so to the convective velocity. It

can be divided into two main sources:

• A temporal source represented by the auto spectra of the pressure fluctuation. It is linked to the turbulence time scale.

• An spatial source represented by the cross spectra of the pressure fluctuation. It is linked to the turbulence length scale.

Due to the randomness of the wall pressure fluctuations, one of the best ways to describe them

is to use weakly stationary distributed random processes: the spatial and temporal properties

of the wall pressure fluctuations are completely given by the space-time correlations. The given

formulas are valid for a boundary layer developing over a flat plate in a low mach configuration:

the boundary layer thickness increases so slowly that it can be considered homogeneous in all

directions and stationary in time. A summary of the different phenomena and the coordinates

definition can be seen in figure 2.4.

The wall pressure fluctuation can be decomposed in two different parts:

• The aerodynamic component characterized by the convective wavenumber: kc= ω/Uc

• The acoustic component characterized by the acoustic wavenumber : kac = ω/c

Uc is the convective velocity and c the sound speed. The acoustic component of the pressure is

related to the waves generated by the turbulence in the boundary layer that propagate in all

space directions and to the contributions of the plate vibration.

For a flow at low Mach number and over a flat plate along the x direction, the pressure

correlation between two points situated at (x, y) and (x + ξ1, y + ξ2) at the time delay τ can be

expressed as :

R(ξ1, ξ2, τ ) = E[p(x, y, t)p(x + ξ1, y + ξ2, t + τ )] (2.10)

Where E is the ensemble average and p(x, y, t) is the fluctuating pressure at the point (x, y)

(28)

Figure 2.4: Vibrating flat plate excited by a TBL. U∞ is the outer fluid velocity and δ is the

boundary layer thickness. Extracted from [1].

1. Autocorrelation in the time domain: R(0, 0, τ )

2. Space correlation : R(ξ1, 0, τ ) in the streamwise direction and R(0, ξ2, τ ) crosswise

Taking the Fourier transformation of these two quantities, one obtains:

1. The auto spectra (known as the single-point wall pressure frequency spectra): φ(ω) =

1 2π

R∞

−∞R(0, 0, τ )e −iωt_dτ

2. The cross spectra (the space/frequency spectra): Sn(ξ1, ξ2, ω) = _2π1

R∞

−∞R(ξ1, ξ2, τ )e −iωt_dτ

The cross correlation power spectral density (PSD) matrix is then the result of the

multiplica-tion of the reference spectrum given here by Goody’s model and the spatial correlamultiplica-tion matrix

calculated here using Corcos’ model: S(ξ1, ξ2, ω) = φ(ω) · Sn(ξ1, ξ2, ω).

The measurement of these quantities is difficult to make for complex geometries and is usually

limited to other functions such as the auto correlation function for example (see [9]).

Fur-thermore, DNS and LES calculations that would allow unsteady flow calculations are still too

costly for industrial applications as it was stated before. As a consequence, semi-empirical

models capable of predicting these spectra have been developed for the last sixty years. These

semi-empirical models allow a fast and reliable way of calculating the auto and cross spectra

(29)

2.5. Auto spectra 15

2.5 Auto spectra

The TBL depends on a variety of lengths, velocities and pressure scales. Usually, the viscous

sublayer gives the first set of length, velocity and pressure scales and the outer layer the second

set. At lower frequencies, the big turbulent structures corresponding to the outer scales of

the boundary layer are the ones that heavily influence the fluctuations of pressure ([10]). For

higher frequencies, it is the small structures that are in the inner part of the boundary layer that

influences the most the behaviour of the pressure field. As a result, for lower frequencies, the

pressure spectra can be scaled using the outer scales such as the boundary layer edge velocity,

the boundary layer thickness and the displacement thickness. For high frequencies, the pressure

spectra can be scaled using the inner layer scales such as the wall shear stress and so the friction

velocity.

The main WPS models that will be reviewed here are: Amiet, Efimtsov, Chase-Howe, Goody

and Smol’yakov.

2.5.1 Amiet

The single point WPS model given by Amiet [11] is:

φ(ω) ρδU3 e = 10−5· F (ω) (2.11) With F (ω) = (1 + ωδ_U∗ e + 0.217( ωδ∗ Ue ) 2_{+ 0.00562(}ωδ∗ Ue ) 4₎−1_{, and 0.1 <} ωδ∗ Ue < 20. This model is a

good starting point as it represents a good mean value but the data is considerably spread out

specially at high frequencies.

2.5.2 Efimtsov

Efimtsov [12] states that the WPS should depend on the mach number, the Reynolds number

and the strouhal number. He considered a number of measurements ranging from M = 0.41 to

(30)

be fully developed and laying on a zero pressure gradient. Efimtsov’s single point wall pressure

spectrum is given by φ(ω) = 0.01τ 2 wδ Ue(1.0 + 0.02(ωδ/Uτ)2/3) (2.12)

This model gives a good mean value at lower frequencies but is not capable of giving the

satisfactory shape at higher frequencies.

2.5.3 Chase - Howe

Howe [13] uses a number of different variables and the Chase theoretical model in order to

propose a new WPS model for low mach configurations: the use of inner and outer scales

allows a better agreement of the numerical model to the experimental data. The model is

proportional to ω2 _{at low frequencies and to ω}−1 _{at higher frequencies but does not include the}

ω−5 slope at even higher frequencies. It states that:

φ(ω) = 2(δ ∗_/U ∞)3(τwω)2 [(ωδ∗_/U ∞)2+ 0.0144]3/2 (2.13)

2.5.4 Goody

The Goody’s model was published in 2004 and is based on the Chase-Howe model ([14]). Goody

added another 19 spectral measurements that covered the Reynolds range 1.4 · 103 _{to 2.34 · 10}4_.

He added the ω−5 decay at high frequencies and modified the exponents and constant to better agree with the experimental data. It is given by:

φ(ω) = 3.0(δ/U∞) 3_(ωτ w)2 [(ωδ/U∞)0.75+ 0.5]3.7+ [(1.1R−0.57_T )(ωδ/U∞)]7 (2.14) Where RT = U 2 τδ U∞ν.

This model has proven to be the best description of the wall pressure spectra and the one that

(31)

2.6. Cross spectra 17

2.6 Cross spectra

The cross spectra gives a quantification of the spatial dependence of the pressure. Bull stated in

1967 [9] that the wall pressure field has two main contributions from two different wave-number

components. The first one is a high wave-number component associated with the turbulence and

it’s laterally coherent over distances proportional to their wavelength. The other component

is associated with the large scale eddies in the boundary layer. Later in 1996, he defined four

different regions: by order and for a fixed frequency, the acoustic ring (where k < ω/c), the

subconvective region (where k ≈ ω/c), the convective region (where ω/c < k < ω/Uc) and

finally the viscous region (where k >> ω/Uc) . The acoustic ring or region corresponds to low

wave-numbers. The subconvective region has a major importance in underwater applications

but no such much for aicraft applications. The convective region corresponds to the region

where the turbulence is convected at the convective speed of the flow Ucand where most of the

energy is concentrated. Finally, the viscous region corresponds to the small scales of turbulence.

These regions can be seen in figure 2.5.

(32)

2.6.1 Corcos

The spatial correlation proposed by Corcos [15] is based on two different coherence lengths: an

streamwise length and a crosswise length. The proposed cross power spectral density is based

on the assumption that the streamwise variation along x is independent from the crosswise

variation along y. As a result, the CPSD proposed by Corcos is:

Sn(ξ1, ξ2, ω) = e−|ξ1|/L1e−|ξ2|/L2e−iωξ1/Uc (2.15)

Where L1 = _αωUc is the coherence length in the streamwise direction and L2 = U_βωc is the coherence

length in the crosswise direction. Experimental observations lead to the mean values of α = 0.1

and β = 0.77. The formulation of the Corcos model in the wavenumber/frequency domain is:

Pn(k1, k2, ω) = αβu2 c π2_ω2_(α2 _{+ (1 −} k1uc ω )2(β2+ ( k2uc ω )2) (2.16)

Corcos model has proven to give good results in the convective region. However, it also

over-estimates the level for small wavenumbers. As it can be seen in the figure 2.6, the experimental

data lies between the Chase and Smol’yakov models and the values are smaller than for the

Corcos model. However, for simplicity reasons it was decided to implement this model in Actran.

2.7 Vibrations

The aero and vibroacoustic modelling is based on a Finite Element approach. The pressure

loading is calculated using the semi-empirical models and is applied on the external skin of

the vibroacoustic model to get the vibration and acoustic results. In this paper, only a weak

coupling between the flow and the structure is considered meaning that the vibration of the

(33)

2.7. Vibrations 19

Figure 2.6: Comparison of cross spectra values with experimental and numerical results. The experimental data is in traingles, dots and cross. More details on [3].

Element Modelling can be expressed as:

X(ω) = A(ω)y(ω) (2.17)

Where A(ω) is the stiffness matrix composed of different blocks: the structure block, the

acoustic block and the coupling block between the acoustical and structural load, y(ω) is the

response vector composed of the structural displacement and the acoustic pressure field and

finally, X(ω) is the load vector composed of a structural and acoustic load.

For a TBL excitation, the pressure fluctuations are random so the system’s response has also a

random character. The problem can then be expressed as a Power Spectral Density formulation:

Sy(ω) = A−1(−ω)Sx(ω)(A−1(ω))T

The pressure is converted into an effective load using a coupling matrix C such that X(ω) =

(34)

Sx(ω) = CTSp(ω)C. So, the system to solve is written as: Sy(ω) = H(−ω)CSp(ω)CTHT(ω).

In this paper, a sampling approach in the space/frequency domain is used to compute the

response of the structure. Firstly, Sp(ω) is decomposed using the Cholesky decomposition:

Sp(ω) = L(ω)LH(ω). Secondly, N samples of Sp(ω) are approximated using N phase angles

with a uniform distribution on [0, 2π] such that the pressure field for the kth realization is

(35)

Chapter 3 Flat plate

The vibroacoustic behaviour of an aluminium flat plate is analysed here. The first challenge

consists in performing a CFD calculation capable of correctly reproducing the flow over the

flat plate. The spatial and temporal representations of the wall pressure spectra are then

cal-culated using semi-empirical models. The next challenge is the reproduction of the mechanical

behaviour of the plate with the pressure field as the exciting force. Finally, the last challenge

is the calculation of the radiated noise. This final step was not investigated thoroughly as it

does not represent a major challenge for VibraTec.

This flat plate study was done in order to verify and validate the methodology that will then

be used for more general cases such as the Canoble project. All the absolute scales on the plots

are hidden or normalized due to confidentiality issues.

3.1 Strategy

The main steps of the method developed are summarised here:

1. A CFD simulation is done using the OpenFOAM software

2. The main aerodynamic parameters of the flow are extracted using Python scripts

(36)

22 Chapter 3. Flat plate

3. The auto and cross spectra of the pressure fluctuations are calculated using Python scripts

4. The displacement of the flat plate is calculated using a FEM solver (MSC Actran) and

the pressure as the exciting force

5. The acoustic field and the noise radiated by the plate is calculated based on the

displace-ment of the flat plate

Only the main steps were named here, several small and dependent steps are also done.

3.2 Configuration

The flat plate is an aluminium plate of length l = 650mm, width w = 300mm and thickness

t = 1mm. It is decoupled via an interface of 20mm. Therefore, it can be considered that the

boundary conditions are free-free for the plate: it is simply supported. A modal analysis of

the plate was done in order to quantify the difference between the numerical model and the

experimental one: the results are discussed later.

To impose a turbulent boundary layer excitation on the plate, the plate is placed on a wind

tunnel at the LMFA (see [17]). The instrumented flat plate is situated at a distance of 1047mm

from the end of the convergent so as to have a fully developed turbulent boundary layer.

The experimental tests were done in the subsonic wind tunnel of the Centre Acoustique at Ecole

Centrale de Lyon in France [17], see figure 3.1 for a visualization of the setup. A subsonic stream

powered by a double stage centrifugal fan (Howden, 800kW , 9 tons) makes it possible to reach

a Mach number M = 0.5 with a mass flow rate of 19.4kg/s in a 30cm× 40cm cross-section and M = 0.8 with a mass flow rate of 10.4kg/s. The incoming air first passes through a settling

chamber with honeycomb and several wire meshes designed to reduce free stream turbulence.

This results in an air flow at ambient temperature with a low background noise and low residual

turbulence intensity, less than 1%. The velocity profiles have been measured with a Dantec

55P01 hot-wire operating in constant voltage mode using a Streamline anemometer (see [18]

(37)

3.3. CFD Simulation 23

point of interest: x = 1272mm (named point 19 afterwards) and 1372mm (named point 22

afterwards).

Figure 3.1: Experimental setup at LMFA. Upper left is the wind tunnel in the anechoic chamber, upper right is the instrumented flat plate and lower is a schematic of the wind tunnel with the main dimensions. The test section if where the flat plate is placed.

3.3 CFD Simulation

The chosen software for the CFD simulation is OpenFoam due to two main reasons: i) it is free

so there is no license limitation and ii) it is opensource, so it is possible to modify the solver and

add, delete or modify coefficients or terms in the equations and it is relatively easy to extract

data.

The goal of this simulation is to reproduce the wind tunnel configuration and to obtain the

heterogenous boundary layer developing along the walls. Once the experimental configuration

has been correctly reproduced it will be possible to extract the main aerodynamic parameters

(38)

3.3.1 Mesh

One of the main issues with CFD calculations is the fact that there is a great number of

param-eters that will affect the accuracy of the solution. Starting from the mesh, it is crucial to have a

good quality mesh as it will affect the rate of convergence (or lack of it for some configurations),

the accuracy of the solution and the calculation time. The starting point is the geometry of

the problem. In this case, the geometry is composed of the walls of the wind tunnel. An “.stl ”

file containing the geometry is imported to OpenFOAM and needs to be correctly reproduced

in order to have a meaningful solution (see figure 3.2).

Figure 3.2: Illustration of a stl geometry (source: cfd.direct ).

The OpenFoam mesher is SnappyHexMesh, a robust but rather complex mesher, which creates

the mesh in a five-step process. The first step will create the background mesh of the problem

(see figure 3.3). Then, the mesher will refine some parts specified by the user (see figure 3.4).

Next, the cells that are not used (inside the obstacle for an external aerodynamic problem or

outside the pipe for an internal problem) are deleted (see figure 3.5). Then, the mesh is snapped

to the surface in order to reproduce the geometry. This step has been proven to be crucial for

the creation of a good quality mesh as it heavily influences the creation of inflation layers. The

final step is the creation of layers near the walls (see figure 3.6).

The mesh are always structured ones because OpenFOAM is only capable of working with this

type of mesh. Structured mesh are usually very efficient for simple geometries but calculations

for more complicated ones are rather inefficient.

The problem we deal with is essentially 2D because the flow will only exhibit 3D effects near

(39)

Figure 3.3: Creation of the background mesh (source: cfd.direct ).

Figure 3.4: Refinement of the mesh (source: cfd.direct ).

Figure 3.5: Elimination of inside cells (source: cfd.direct ).

Figure 3.6: Surface snapping and layer addition (source: cfd.direct ).

tunnel the crosswise flow can be assumed to be homogeneous. This was confirmed by the data

provided by the LMFA: the velocity profiles in the crosswise direction gave similar results. The

CFD calculation is a 2D simulation which is done in order to simplify the CFD simulation by

having less cells : the mesh generation and the calculation are faster.

(40)

condition to the patch whose normal direction is aligned to the geometric direction for which

no solution is required..

Another parameter of influence is the size of the cells. To have a good quality mesh, the cells

should have an aspect ratio of 1, this means that the ratio of the length and width and thickness

should be close to 1. This will be crucial for the snapping phase and will allow the creation of

good quality layers. Python scripts were developed in order to give the layer parameters based

on the freestream cell size.

The first layer thickness of the mesh is of great importance as it will determine the accuracy

of the simulated boundary layer. Depending on the problem, one can choose to fully resolve

the boundary layer down to the wall or to use wall functions that will model the near-wall

behaviour of the flow. The main parameter to judge the applicability of a wall function is the

dimensionless wall distance y+_{of the mesh. When using wall functions, the first grid cell should}

be somewhere between 30 < y+_{< 300. When no wall functions are used y}+ _{should be close to}

1, for the k − ω model.

Figure 3.7: Illustration of boundary layer solving strategies (source: Leap Australia).

Of course, the dimensionless wall distance is a function of a priori unknown aerodynamic

quantities because some quantities (as the friction velocity uτ) need to be calculated in order

to estimate the dimensionless wall distance. It is possible to estimate this quantity based on

experimental correlations. The method is the following:

1. Calculate the Reynolds number based on the freestream velocity and the length scale of

(41)

problem)

2. Estimate the skin friction coefficient Cf based on an experimental correlation.

3. Compute the wall shear stress defined by τw = Cf1₂ρUf ree2 and the friction velocity

4. Compute the wall distance

A comparison of different correlations was done in order to asses the influence of the formula.

The comparison can be found on table 3.1.

Correlation Desired y+ Cf ds (mm) Actual y+

Raymer 30 0.0031 0.176 32 KTH course 30 0.0030 0.178 32.3 Von Karman (see [19]) 30 0.0035 0.165 29.98 Schlichting (see [20]) 30 0.0027 0.191 34.7

Table 3.1: Comparison of correlations formulas

The y+ _{estimation is close and coherent for all correlations. As a result, there is no need to}

choose a correlation instead of another one for industrial cases.

Figure 3.8: The geometry that was used for this study. One should note the inlet offset.

The mesh that was used for this study is presented in figure 3.8. The cells can not directly be

seen as their size relative to the global geometry is very small. It should however be noted that

there is an offset in the geometry: when comparing figure 3.8 with 3.1, the used geometry has

an inlet offset of 2.5D = 1.4m, where D = 0.56m is the inlet diameter. This was done in order

to better initialize the calculation by growing a boundary layer along the inlet tunnel and thus

making the numerical simulation closer to the experimental configuration. This proved to be

very efficient in order to have comparable results with the experimental data. In figure 3.9, one

can see the refinement of the mesh close to the wall. In this case, the expansion rate of the

(42)

Figure 3.9: Zoom on the near wall zone of the mesh.

3.3.2 Computation

Once a good quality mesh has been created it is then necessary to choose the solver, the

nu-merical schemes, the initial conditions and the boundary conditions. The boundary conditions

for this case are as follows:

• the walls of the wind tunnel are modelled as walls with no-slip conditions meaning that the flow will be attached to the walls.

• the inlet has fixed velocity U = (Ux, Uy, Uz) = (Ux, 0, 0) as the flow is supposed to depend

only on the axial direction and a “zero gradient” pressure boundary condition meaning

that the gradient of the pressure is null in the perpendicular direction to the inlet patch

• the outlet will have a zero gradient velocity meaning that the velocity gradient is null in the perpendicular direction to the outlet patch and a fixed relative static pressure of

Poutlet = 0.

(43)

The only quantities that remain to be calculated are the turbulent ones: the turbulent kinetic

energy and the specific turbulent dissipation rate for a k − ω turbulence model. The equations

used for calculating these quantities are: k = 3₂(U I)2 _{and ω = C}−0.25 µ

√ k

l . Where I is the

turbulence intensity, µt the turbulent viscosity and l the turbulent characteristic length. The

turbulent kinetic energy and the specific turbulent dissipation rate were initialized using an

inlet turbulence intensity of I = 1% as in the experimental setup.

The solver potentialFoam is used to initialize the flow field. It is a solver that works under the

assumption of a potential flow (meaning that the velocity vector is the gradient of a scalar field

named the velocity potential: V = grad(φ) ). It is very useful to provide a conservative initial

field which is very useful to avoid instabilities. Running this solver for 5 iterations has proven

to provide good initial fields which rapidly converged to a steady solution.

The numerical scheme used for all quantities was a second order upwind scheme. Using an

upwind scheme instead of a centered one is a good choice for stabilizing such a flow field

dominated by the convection. It is a second order one in order to have an accurate solution

and avoid the dissipation caused by first order schemes. No divergence problem was observed

as the mesh were of good quality and the boundary conditions are well defined.

The solver used for the RANS simulation is the simpleFoam solver. This solver is based on

the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm and solves the

continuity equation and the momentum equation. At least 4000 and up to 10000 iterations are

then needed in order to have a converged solution (see for example figure 3.10).

3.3.3 Post-Processing

Once the calculation is finished, OpenFOAM saves the different fields in a raw format. It is then

needed to use scripts in order to extract the quantities of interest. There is essentially a single

quantity that needs to be extracted: the velocity vector along the vertical axis y and the wall

shear stress resulting from the “no-slip” condition. With these informations it is then possible

to calculate the boundary layer thickness, the displacement and momentum thicknesses and

(44)

Figure 3.10: Residuals for U0 = 30m.s−1 and y+ = 30

1. Create a file containing the points of interest where these two quantities will be extracted

in the OpenFOAM syntax.

2. Read the raw data from OpenFOAM and organize it in order to make it easier to work

with and to make it easier to interact with other softwares.

Boundary layer thickness

The definition of the boundary layer thickness is confusing as the velocity profile is continuous

and it’s difficult to establish the thickness of it. A good indicator is the use of δ99 defined by

δ99 = y(u = 0.99Ue). The problem with this convention is the definition of the outer velocity

Ue. After discussion with the LMFA it was decided to choose Ue= max(U ). The procedure to

calculate the boundary layer thickness is then:

1. Interpolate the velocity as a function of the vertical distance to the wall

2. Calculate this function on a vector containing the same spatial limits (ymin = 0 and

(45)

3. Search for the first point where u > 0.99 · Ue

This procedure gave good results for this configuration.

The procedure for calculating the displacement and momentum thickness is similar to the

calculation of the boundary layer thickness:

1. Interpolate the velocity as a function of the vertical distance to the wall

2. Calculate this function on a vector containing the same spatial limits (ymin = 0 and

ymax = 100mm) but with a greater number of points

3. Integrate the calculated velocity using the definitions of the displacement and momentum

thickness.

Finally, the friction velocity depends on the wall shear stress that was directly extracted from

OpenFoam. However, the value was not coherent with the experimental values provided by the

LMFA. It was decided to calculate the wall shear stress by fitting the velocity profile to the

Musker profile ([21]).

The Musker profile uses the eddy viscosity model and is provided by the solution to the following

equation:

du+

dy+ =

(y+₎2_{/κ + 1/s}

(y+₎3_{+ (y}+₎3_{/κ + 1/s} (3.1)

Musker showed that s = 0.001093 by fitting the solution u+(y+) of the equation to the log-linear region.

The velocity profile is firstly normalized into wall units (u+ and y+) by selecting the value of

wall shear stress (and so friction velocity) that best fits to the Musker profile. Practically, this

process was automated using Python scripts by searching for the minimum difference between

the normalized data and the model equation using the Nelder-Mead simplex algorithm which

proved to be very efficient.

(46)

sampling may not capture this point. The constants for the velocity profile were the

commonly-used boundary layer constants (κ = 0.41, B = 5.1) following the advice of [21]. This method

proved to give good results coherent with the experimental data.

3.4 FEM calculation

A FEM model was created in Actran to simulate the behaviour of the flat plate subject to

the flow. The model can be seen in figure 3.11: the black points represent the position of the

accelerometers, it is also possible to see the flat plate in the middle (with a more refined mesh)

and the interface on the outer (the mesh is coarser).

Figure 3.11: FEM mesh model of the flat plate

3.4.1 Setup

The flat plate is modelled as a thin shell of thickness 1mm and made of aluminium (Young

modulus E = (70 + 1.4i)GP a, Poisson ratio ν = 0.3 and density ρ = 2700kg.m−3). The air is supposed to be incompressible with a sound speed of c = 340m.s−1 and a density of 1.225kg.m−3. The flat plate is simply supported on a 20mm border by an interface which is

(47)

3.4. FEM calculation 33

modelled as a thin shell of thickness 1mm. The accelerometers are modelled as point masses

of mass 0.003kg. As the problem can be considered an acoustic radiation in a half-space, the

radiation surface is modelled by a Rayleigh surface. The acoustic radiation problem can be

formulated as: ∇p(r) + k2_{p(r) = 0 where r is a point, k the wavenumber and p the acoustic}

pressure. The boundary conditions are: _∂n∂p = ρω2_u

n and limr→∞r(_∂n∂p + ikρ) = 0 with r

a spherical radial distance. Then, the acoustic pressure in the half-space Σ is given by the

following boundary integral:

p(r) = −ρω2 Z Σ un(r0) e−ikR 2πRdΣ(r) (3.2)

The element size of the flat plate is ∆x = ∆y = 5mm and the support elements are of

size ∆x = 10mm and ∆y = 5mm. The mesh is composed of 7056 elements. These values

were chosen so as to have a valid mesh up to 1250Hz. However, the adaptivity method (see

Actran user guide) was used in order to have 6 elements per wavelength for every frequency

up to 5623Hz. It is considered 180 logarithmically distributed frequencies between 355Hz and

5623Hz. As it will be seen, this method yields good results up to 2000Hz.

Table 3.2: Modal base comparison

Experimental fexp (Hz) Numerical fnum (Hz) ∆ (%)

59 77.3 23.7 95 93.3 1.82 132 122.3 7.93 168 164.1 2.38 189 204.5 7.58 216 215.2 0.37 226 220.5 2.49 265 247.8 6.94

The flat plate is modelled as simply supported (Free-Free conditions). The validity of the

model and its suppositions can be verified by comparing the modal base of the plate: the table

3.2 shows the comparison between the numerical eigenfrequencies and the experimental ones.

The relative difference ∆ is equal to ∆ = (fexp − fnum)/fnum. The values are relatively close

(48)

match the experimental with the numerical modal base was made due to the great number of

modes for this plate, as it will be seen this creates some differences between the experimental

measurements and numerical data at lower frequencies. A few examples of the plate modes can

be seen in figure 3.12.

Figure 3.12: Main modes of the flat plate with their corresponding frequency at their right.

3.4.2 Calculation

With this setup, a direct calculation was done on MSC Actran. A direct calculation is usually

not very efficient as the equation 2.17 has to be solved for every desired frequency, meaning

that the system of equations has to be inversed. However, due to the low number of frequencies

we calculate, a modal approach was not interesting here due to the fact that the plate has a

high number of modes.

The interest of the method here presented is the capacity of customizing the generation of the

reference spectrum and the spatial correlation PSD matrix: it would be possible to use any

available semi-empirical model and tune its coefficients. Indeed, Python scripts were developed

in order to calculate those quantities and use them as an input for MSC Actran. It is then

possible to calculate the spectrum and the spatial correlation using 3D fields that are extracted

from OpenFOAM. The figure 3.13 shows an example of the spatial correlation PSD matrix

(49)

3.4. FEM calculation 35

can be seen that the values of the spatial correlation PSD matrix are maximal in the diagonal

and decrease as the distance to the diagonal increases which should be related to the coherence

of the pressure in the space domain.

Figure 3.13: Cross correlation matrix for U0 = 30m.s−1 and f = 355Hz. The axis represent

the FEM node ID.

First, Actran calls a python script in order to retrieve the auto and cross spectra. These spectra

are calculated using the CFD data that is firstly interpolated over the CFD domain. Secondly,

the interpolated function (boundary layer thickness, friction velocity etc) is applied to the FEM

domain at every node of the mesh. As a result, one has the aerodynamic parameters at every

point of the mesh and it is then possible to calculate the spectra.

This process is shown in figures 3.14 and 3.15 with an example of a ”.edat ” file for Actran. The

boundary layer parameters such as the friction velocity or the outer velocity for example are

calculated over the FEM mesh using an interpolation of the fields coded in Python (3.14). The

wall pressure models are given to Actran using a Python script. One should note that in figure

3.15, the model that was used is not mentioned: all the calculations are done in the Python

(50)

Figure 3.14: Edat file showing the implementation of the boundary layer parameters.

Figure 3.15: Edat file showing the wall pressure spectra construction for Actran.

The interpolation used is done using the Python function LinearNDinterpolator of the package

scipy.interpolate that performs a linear interpolation over the dataset. A linear interpolation is

(51)

3.5. Results 37

3.5 Results

In this part, the results from both the numerical analysis (CFD and FEM) and the experimental

campaign (hot-wire and accelerometers) are presented in order to conclude of the validity and

accuracy of the method.

3.5.1 CFD and experimental validation

CFD calculations were made with an inlet velocity of Ux = 30m.s−1, Ux = 50m.s−1 and

Ux = 60m.s−1 and for each velocity three meshes were used with y+ = 1, y+ = 30 and

y+ _{= 100. Quantities of interest were extracted along the axial direction every 1mm and at}

two particular points of interest where the hot-wire anemometry was performed (previously

defined as point 19 and point 22).

It is presented the comparison between the velocity profiles from the experimental analysis and

the CFD calculations for two different points situated at an axial distance of 100mm. Figures

3.16 and 3.17 show the normalized velocity profiles for a velocity of U0 = 30m.s−1. The

numer-ical profile that best agrees with the experimental data profile and the theoretnumer-ical boundary

layer laws is the one for y+ _{= 1. Figures 3.18 and 3.19 show the normalized velocity profiles for}

a velocity of U0 = 50m.s−1. The numerical profile that best agrees with the experimental data

profile and the theoretical boundary layer laws is once again the one for y+ _{= 1. However good}

the match between the normalized profiles seems, one should take a closer look at the calculated

values that will be used in order to calculate the wall pressure spectra. The tables 3.3, 3.4, 3.5

and 3.6 show the inner and outer parameters of the boundary layer for the experimental and

numerical case at the two reference points and for two different velocities: Ux = 30m.s−1 and

Ux = 50m.s−1 respectively.

In this case, the best overall prediction is not always realized with a y+ _{= 1 mesh. Indeed, as it}

can be seen in the figures 3.16 to 3.19, the boundary layer thickness is very close for the three

meshes (around y+ _{= 2000) and the log law is well reproduced. As a result, the displacement}

(52)

are small compared to the whole boundary layer so its influence in the integration of the profile

is minor. Finally, as the friction velocity is calculated through an optimization process in order

to fit the Musker model the values are also close for the three meshes. The CFD simulation

has been validated and can be used for other inlet speeds and conditions.

Figure 3.16: Point 19 at U0 = 30m.s−1

(53)

3.5. Results 39

Figure 3.18: Point 19 at U0 = 50m.s−1

(54)

40 Chapter 3. Flat plate Case Ue Uτ δ δ1 δ2 Experimental 31.8 1.18 29.54 4.04 3.14 CFD y+_{= 1} _30.78 _1.13 _26.66 _4.04 _3.00 CFD y+_{= 30} _33.02 _1.14 _29.03 _4.36 _3.15 CFD y+= 100 32.98 1.12 28.92 4.41 3.08 Table 3.3: Point 19 at U0 = 30m.s−1 Case Ue Uτ δ δ1 δ2 Experimental 31.87 1.19 34.71 4.33 3.34 CFD y+_{= 1} _30.71 _1.12 _27.99 _4.24 _3.16 CFD y+= 30 33.95 1.10 30.02 4.56 3.31 CFD y+= 100 32.91 1.11 29.77 4.61 3.24 Table 3.4: Point 22 at U0 = 30m.s−1 Case Ue Uτ δ δ1 δ2 Experimental 52.1 1.89 34.4 3.37 2.58 CFD y+_{= 1} _52.15 _1.86 _26.94 _3.53 _2.74 CFD y+= 30 52.71 1.76 26.51 3.95 2.94 CFD y+_{= 100} _52.65 _1.70 _26.7 _4.12 _2.85 Table 3.5: Point 19 at U0 = 50m.s−1 Case Ue Uτ δ δ1 δ2 Experimental 51.8 1.89 30.8 3.6 2.8 CFD y+_{= 1} _51.85 _1.86 _28.67 _3.81 _2.94 CFD y+= 30 52.59 1.75 27.58 4.15 3.09 CFD y+_{= 100} _32.52 _1.69 _27.92 _4.31 _2.99 Table 3.6: Point 22 at U0 = 50m.s−1

3.5.2 CFD results

The CFD calculation yields a good overall result. It is then possible to use the results first to

calculate the flow for other inlet velocities and secondly to use them as an input for the

semi-empirical models. It was decided to simulate the flow for an inlet velocity of Ux = 60m.s−1 so

a mach number of around M = 0.2.

The results are presented in the tables 3.7 and 3.8 and the profiles can be seen in figure 3.7

(55)

3.5. Results 41 Case Ue Uτ δ δ1 δ2 CFD y+= 1 63.53 2.21 25.16 3.56 2.72 CFD y+_{= 30} _63.63 _2.11 _26.09 _3.82 _2.86 CFD y+_{= 100} _63.41 _2.03 _27.30 _4.13 _2.91 Table 3.7: Point 19 at U0 = 60m.s−1 Case Ue Uτ δ δ1 δ2 CFD y+_{= 1} _63.8 _2.19 _26.03 _3.74 _2.86 CFD y+_{= 30} _63.48 _2.09 _27.10 _4.01 _3.00 CFD y+= 100 63.41 2.03 27.30 4.13 2.91 Table 3.8: Point 22 at U0 = 60m.s−1 Figure 3.20: Point 19 at U0 = 60m.s−1

(56)

(57)

3.5. Results 43

Once again, the calculated values are very close for the three different meshes. As a conclusion,

a mesh of y+ _{= 30 should be used for industrial cases in order to reduce the number of cells}

and so reducing the time of the solution. For simple and incompressible flows, the accuracy of

the prediction does not seem to be affected by the refinement of the mesh. Indeed, the wall

functions used for y+ _{= 30 and y}+ _{= 100 give reasonable results. It should be noted that for}

more complex flows exhibiting recirculation zones and pressure gradients, a y+ = 1 should be used in order to have a good boundary layer profile.

It is also presented here the pressure (figures 3.22, 3.23 and 3.24) and velocities (figures 3.25,

3.26 and 3.27) fields in the whole wind tunnel. As it can be seen the boundary layers develop

linearly from the end of convergent up to the outlet and the pressure field in the test section is

almost constant due to the diverging upper wall.

Figure 3.22: Pressure field at U0 = 30m.s−1

(58)

Figure 3.24: Pressure field at U0 = 60m.s−1

Figure 3.25: Axial velocity field at U0 = 30m.s−1

Figure 3.26: Axial velocity field at U0 = 50m.s−1

(59)

3.5. Results 45

It is shown here the results of friction velocity (figure 3.28), boundary layer thickness (figure

3.29), displacement thickness (figure 3.30) and momentum thickness (figure 3.31) along the

wind tunnel for three different velocities. It should be noted that the evolution of all quantities

is linearly dependent on the axial direction and has the same slope for all three velocities.

(60)

Figure 3.29: Boundary layer thickness along x.

(61)

3.5. Results 47

(62)

The interest of this calculation is to compute 2D space varying fields that will be used to

calculate the response of the flat plate and quantify the difference between an homogeneous

boundary layer load and an heterogeneous boundary layer load. The shape factor was calculated

in the test section and has a mean value of H = 1.34, so the boundary layer is indeed turbulent

([20]).

3.5.3 Microphone array results

A flat array of microphones was installed in the lower wall of the wind tunnel so as to measure

pressure using 40 different and remote microphone probes: n = 24 microphones were placed in

the streamwise direction and m = 16 in the cross wise direction, the spacing between them is

irregular. The pressure signals are measured simultaneously with a sampling frequency of 50.0

kHz. The experimental setup was completely organized by the LMFA (see [22], [23] and [24]

for more information).

For stationary and stochastic signals, the cross spectral density is defined by:

Rp,p(x, r, ω) = lim T →∞

2π

T E[p(x, ω)p(x + r, ω)] (3.3)

Where r is the separation between two points and E is the ensemble average. Practically,

the frequency cross spectra is calculated using the Welch periodogram method. A rectangular

window is used and applied to every block: the average is done on 12501 blocks with 50%

overlapping. The cross spectra matrix can be seen in figure 3.32 for a frequency of f = 100Hz,

figure 3.33 for a frequency of f = 200Hz, figure 3.34 for a frequency of f = 500Hz and

figure 3.35 for a frequency of f = 1000Hz. The axis represent the number of the microphone

in the array: the first 24 are the streamwise microphones and the last 16 are the crosswise

microphones. Only the streamwise microphones will be used here. The values are given in dB,

the colormap can not be shown for confidentiality issues.

Once the frequency cross spectra matrix has been calculated it is possible to calculate the

(63)

3.5. Results 49

Figure 3.32: Cross spectra matrix, 100Hz. Figure 3.33: Cross spectra matrix, 200Hz.

(64)

Figure 3.36: Wavenumber-frequency diagram for U0 = 10m.s−1 at left and U0 = 20m.s−1 at

right

Figure 3.37: Wavenumber-frequency diagram for U0 = 30m.s−1 at left and U0 = 40m.s−1 at

right written: Φ(k, ω) = 1 4π2 Z Z R Rp,p(x, r, ω)e−ik·rdr (3.4)

Practically, the CPSD matrix was calculated by defining a vector of wavenumbers regularly

spaced with a step lower than the antenna resolution and calculating the discrete fourier

trans-formation along the streamwise direction. It was chosen ∆k = 5 rad.m−1. The discrete formula that was used here is:

Φ(k1, y, ω) = 1 2π n X i

Rpi,pj(xi, yj, ω)e−i·k1·ξi (3.5)

Aero and vibroacoustical prediction of the noise generated by turbulent boundary layers

KTH Royal Institute of Technology

Department of Vehicle Engineering