impedance boundary in a time-domain CAA method
E. van der Pol
Master Thesis in Applied Mathematics
August 2012
method
Summary
In this thesis the modelling of a ground impedance boundary condition in a time-domain com- putational aeroacoustic method is investigated. A ground impedance boundary is a boundary that absorbs part of the incoming acoustic wave and reflects the rest. In the literature several models are proposed to model the quantity of absorption, the impedance, for a range of ground types. All of the proposed models are defined in the frequency domain, but the computational aeroacoustic method calculates in the time domain. Unfortunately, the transformation to the time domain is not straightforward, and therefore approximations are proposed. One ap- proximation is a broadband approach and is described by poles, another is a single-frequency approximation.
There are several ways to discretize the pole approximation. In the literature a recursive implementation is proposed. However, in this thesis a new implementation is proposed. The implementation is given as a system of ordinary differential equations. This system can be solved with the same fourth-order time-integration method, as is used in the compu- tational aeroacoustic method, which makes the entire system fourth-order accurate. Some two-dimensional test cases show the influence of the ground impedance boundary condition on a sound wave.
Master Thesis in Applied Mathematics Author: E. van der Pol
Supervisor(s): A.E.P. Veldman Second supervisor: A.C.D. van Enter
External supervisor: J.C. Kok and W. Rozema Date: August 2012
Institute of Mathematics and Computing Science P.O. Box 407
9700 AK Groningen The Netherlands
Executive summary
Modelling of a ground impedance boundary in a time-domain CAA method
Report no.
NLR-TR-2012-349 Author(s)
E. van der Pol Classi cation report Unclassi ed
Date August 2012 Knowledge area(s) Computational Physics &
Theoretical Aerodynamics Aeroacoustics & Experimental Aerodynamics
Descriptor(s)
Ground impedance boundary condition
Computational aeroacoustics Problem area
Aircrafts make noise during take- off and landing. This can be quite annoying to people living nearby airports, but also a danger to the en- vironment. Sound walls and other objects have been placed nearby noisy areas such as airports or high- ways to block part of the sound, but these objects only reflect the sound. This holds for most man- made surfaces that are basically solid and they are called acousti- cally hard. Naturally most grounds are not acoustically hard, but sur- faces covered with grass or other vegetation. These grounds absorb some of the incoming sound and therefore reflect only part of it.
NLR uses a computational aeroa- coustic (CAA) method to accurately predict the propagation of sound in a fluid. So far the method can only handle acoustically hard bound- aries. As said before, most surfaces are not acoustically hard, so an im- plementation of acoustically soft boundaries is wanted. Acoustically soft boundaries can be modelled as impedance boundaries where the impedance describes the quantity of the absorption.
Description of work
This investigation concers the modelling of impedance bound- aries and the implementation of the impedance boundaries in the computational aeroacoustic method of NLR. The impedance describes how an impedance boundary re- flects and absorbs sound waves for a certain ground type and fre- quency. Therefore the impedance is defined in the frequency domain.
For all acoustically soft surfaces an impedance can be prescribed, so models are necessary to deter- mine what the impedance is for different ground types. Since 1949 models can be found in the liter- ature that describe the impedance for different ground types and fre- quencies. Some of these models are based on analytical derivations and others are based on experimental values. For further investigation the experimentally-derived and physically-valid Miki model was chosen.
The computational method of NLR predicts propagation of sound in the time domain and the impedance is defined in the frequency domain.
Therefore the impedance needs to be transformed to the time do-
main with a Fourier transform.
Most of the models from the lit- erature are not easily or not at all transformable to the time domain.
Therefore approximations have to be made. There are different ways of approximation the impedance, including a single-frequency and a broadband approximation. In this case the broadband approximation that is described by poles is chosen for further investigation.
Once a choice for a model and approximation are made the impedance is discretized to im- plement into NLR’s computational method. In the literature some dis- cretizations can be found, but a new one based on ordinary differential equations is derived. This method is then applied to a one-dimensional finite-volume method in Matlab.
With this method the accuracy and stability of the system with an impedance boundary condition are
investigated.
Results and conclusions It is established that the new dis- cretization of the impedance boundary keeps the system fourth- order accurate and stable for small (enough) time steps. After expend- ing the discretization to higher order, the impedance boundary is applied into NLR’s computa- tional method. With this method two-dimensional test cases were simulated. The first test case simu- lates a point source and investigates the sound pressure level. Another test case initializes an acoustic pulse and simulates its propagation over a large domain.
Applicability
Now that the impedance boundary is implemented in NLR’s compu- tation method, calculations with acoustically soft surfaces can be carried out.
NLR-TR-2012-349
Modelling of a ground impedance boundary in a time-domain CAA method
E. van der Pol
No part of this report may be reproduced and/or disclosed, in any form or by any means, without the prior written permission of the owner.
Customer National Aerospace Laboratory NLR
Contract number ----
Owner National Aerospace Laboratory NLR
Division Aerospace Vehicles
Distribution Limited
Classification of title Unclassified August 2012 Approved by:
Author Reviewer Managing department
Contents
1 Introduction 5
2 Modelling of aeroacoustics 6
2.1 Linearized Euler Equations 6
2.2 Boundary conditions 7
2.2.1 Solid wall 8
2.2.2 Inflow and outflow boundaries 8
2.2.3 Partially absorbing surfaces 9
3 Impedance ground surfaces 9
3.1 Frequency vs. time domain 9
3.2 Material dependence 10
3.3 Hard-backed layered ground 11
3.4 Physical validity 11
3.5 Impedance models 12
3.5.1 Zwikker-Kosten 12
3.5.2 Delany-Bazley/Miki 15
3.5.3 Other models 18
3.6 Approximations 19
3.6.1 Poles 19
3.6.2 MSD+ 20
3.7 Conclusions 22
4 Numerical method 23
4.1 Computational problem 23
4.2 Time integration 25
4.3 Spatial schemes and artificial diffusion 26
4.3.1 Artificial diffusion 29
4.4 Discretized boundary conditions 29
4.4.1 Solid wall 30
4.4.2 Symmetry boundary 31
4.4.3 Impedance boundary with recursive convolution 31
4.4.4 Impedance boundary with a system of ODEs 33
4.5 One-dimensional test case 34
4.5.1 Accuracy 37
4.5.2 Stability 39
5 Two-dimensional problems 43
5.1 Single frequency 45
5.2 Gaussian pulse 48
6 Conclusions and recommendations 51
References 53
Appendix A ENSOLV routines 55
A.1 pol.h 55
A.2 bcpole.F 55
A.3 dim.h 61
A.4 bcond.F 61
A.5 rdbcd.F 61
A.6 updbcd.F 63
Appendix B Discretization methods for the MSD+ approximation 64
B.1 Integration discretization 64
B.2 Discretization with auxiliary parameters 65
B.3 System of ODEs 66
1
Introduction
An aircraft emits a lot of noise. This is not a problem for the environment when an aircraft is fly- ing at 30.000 feet, but during takeoff and landing an aircraft can make so much noise that it may be annoying to people and also an environmental hazard. All kinds of objects and effects can strengthen or damp the amount of noise that is heard. Think of a sound wall placed nearby busy highways to damp the noise to residential neighbourhoods, or when the wind is directed towards you, the noise seems to strengthen. When there is a solid ground like a runway, the noise reflects upon this runway. However looking around an airport like Schiphol, there is a lot of grass and vegetation. One of the properties of ground covered with grass or vegetation is that part of the noise is absorbed by the ground. This property is even stronger in a layer of snow covering the ground. Therefore it is very interesting to model the influence this property has on the propaga- tion of sound.
NLR uses a computational aeroacoustic method to model the propagation of sound, while taking into account ground interactions and atmospheric effects. At the moment the method models sound propagation with a perfectly reflecting ground boundary condition. However most of the time we do not deal with a solid ground, but with a vegetation covered ground such as grass, or in the winter sometimes with a ground covered with snow. These ground types can be modelled as an impedance boundary condition.
The property that indicates the amount of absorption and reflection of an incoming acoustic wave is called the impedance. In reality there is a large range of ground types which all have different impedances. The impedance depends on the material properties of the ground type and therefore models are used to determine the impedance for a specific ground type. In the literature a wide range of models can be found to describe the impedance of ground types. Some models were derived analytically and others were fitted to experimental data.
The impedance is not only dependent on the ground type, but also on the frequency of the incom- ing acoustic wave. The impedance is therefore defined in the frequency domain. However the goal of this report is to implement an impedance boundary into the CAA method used by NLR which calculates the propagation of sound in the time domain. Therefore the impedance needs to be transformed to the time domain with the Fourier transform. Most models cannot be (easily) transformed to the time domain and have to be approximated so that they can be transformed.
This can be done for a certain frequency (single band) or for a frequency band (broadband).
Several models and approximations are described and discussed in this report and a choice is made for a model and approximation. There are different methods to numerically discretize the
chosen approximation. These methods have a different set-up and their accuracy and stability are discussed with a one-dimensional simulation.
When the computational method for an impedance boundary is implemented, some test cases can be simulated. In the final chapter of this report some two-dimensional simulations can be found which were simulated with the aeroacoustic method of NLR. The first test case simulates a source that produces a continuous acoustic wave of a certain frequency. The second test case shows what an impedance boundary does to a Gaussian pulse, which contains a band of frequen- cies.
2
Modelling of aeroacoustics
Aeroacoustics is the study of the propagation of sound in a fluid. Sound can be heard due to pres- sure perturbations in fluids such as air. Therefore the study of the propagation of sound can be seen as a fluid dynamics study. The equations that describe the propagation of sound can be de- rived from the Euler equations for fluid dynamics problems with inviscid and compressible flow.
This will be shown in the next section.
2.1 Linearized Euler Equations
The Euler equations for inviscid and compressible flow in conservation form are given as
∂ρ
∂t + ∂(ρuj)
∂xj = 0
∂(ρui)
∂t + ∂(ρuiuj)
∂xj
+ ∂p
∂xi
= 0
∂(ρE)
∂t + ∂(ρHuj)
∂xj = 0. (1)
where p is the pressure, u the velocity vector and ρ the density. Also E is the total energy per unit mass and H = E + pρis the total enthalpy per unit mass (Ref. 7). For a calorically perfect gas, the pressure can be eliminated from the Euler equations by using the following relation
p = (γ − 1)(ρE −1
2ρuiui). (2)
Here γ is the ratio of specific heats and is written as γ = ccp
v, with cv and cp respectively the heat capacity at constant volume and the heat capacity at constant pressure. For air, γ is equal to 1.4.
The density, the velocity and the energy can be split in a time-independent mean flow and time- dependent perturbations of the mean flow as
ρ(x, t) = ρ0(x) + ρ0(x, t) u(x, t) = u0(x) + u0(x, t)
E(x, t) = E0(x) + E0(x, t). (3)
Here the subscript 0 denotes the mean flow and the primes denote the perturbations. The param- eter γ again appears, here in the relation for the speed of sound c0. This is given as c0 =qγp
0
ρ0 , where p0 is the pressure and ρ0 the density of the mean flow. The internal energy e for a calori- cally perfect gas is dependent on the temperature T by e = cvT .
As said before, aeroacoustics is the study of the time-dependent perturbations of the mean flow.
The total flow has to satisfy the Euler equations. When substituting the flow variables (3) into the Euler equations (1), and subtracting the steady Euler equations which have to hold for the mean flow variables, a system for the perturbations is left. Assuming that the perturbations are small compared to the mean flow variables, such that all quadratic or higher order terms of the pertur- bations can be neglected, then the perturbations have to satisfy the linearized Euler equations (LEE) in conservation form, given by
∂ρ0
∂t + ∂
∂xj
(u0jρ0+ ρ0u0j) = 0
∂(ρui)0
∂t + ∂
∂xj(u0j(ρui)0+ (ρ0u0i)u0j) + ∂p0
∂xi = 0
∂(ρE)0
∂t + ∂
∂xj(u0j(ρH)0+ (ρ0H0)u0j) = 0 (4)
with
(ρuj)0 = ρ0u0j+ u0jρ0 (ρE)0 = ρ0E0+ E0ρ0
p0 = (γ − 1)((ρE)0− u0i(ρui)0+ 1
2(u0iu0i)ρ0)
(ρH)0 = (ρE)0+ p0. (5)
2.2 Boundary conditions
Sound propagation is defined on a domain, thus boundaries for that domain have to be set. These boundaries can have different types such as a solid wall or an outflow boundary. A few of these boundaries will be described here for a 3-dimensional domain.
2.2.1 Solid wall
A solid wall is called acoustically hard, i.e. there can be no flow through the solid wall, so the normal momentum has to be zero. This is called a slip boundary condition. Say a vector n points out of the boundary then the condition can be written as
(ρu)0· n = 0. (6)
2.2.2 Inflow and outflow boundaries
Another possibility is an open boundary. At this boundary either the incident wave can enter or leave the domain. This open boundary can be modelled as an inflow or outflow boundary, where respectively an incident wave enters or leaves. The inflow and outflow boundary condition are modelled as characteristic boundary conditions based on the one-dimensional characteristic the- ory. This means that the solution is expressed in characteristic variables or Riemann invariants.
When n, s and t form an orthonormal system, with n pointing outside of the boundary, the Rie- mann invariants are given as
s0 = cv p0
p0− cp ρ0
ρ0 u0s= s · u0
u0t= t · u0 Ro0 = 1
2n · u0+ 1 2ρ0c0
p0 R0i = 1
2n · u0− 1 2ρ0c0
p0 with c0, cvand cpas before.
The perturbations of the mean flow can be expressed in terms of the Riemann invariants as ρ0 = 1
c20p0−ρ0
cps0
u0 = (R0o+ R0i)n + u0ss + u0tt p0 = ρ0c0(R0o− R0i).
For the inflow and outflow boundaries the incoming Riemann invariants have to be prescribed.
The inflow boundary has different incoming Riemann invariants than the outflow boundary. The incoming Riemann invariants at a far field inflow boundary are s0, u0s, u0tand R0iand are set to zero. At a far field outflow boundary the incoming Riemann invariant is R0iwhich is set to zero.
Because these characteristic boundary conditions are known to lead to significant reflections of outward-going waves in 2D and 3D, a buffer zone is added. This however will not be discussed here and more information can be found in Ref. 7.
2.2.3 Partially absorbing surfaces
In reality not all the boundaries are perfectly reflecting like a solid wall, or open like an outflow boundary. Most natural surfaces, such as grasslands, are porous surfaces and have the property to absorb a part of the incoming sound and reflect the rest. These surfaces are called impedance sur- faces or boundaries and can be seen as acoustically soft. Impedance is the measure of reflection and absorption by the boundary and is related to the pressure and velocity of the flow. In the next chapter it is explained how the impedance boundary condition is defined, how the impedance is modelled for different porous materials and surfaces, and how the impedance can be transferred to the time domain.
3
Impedance ground surfaces
The impedance is defined as the ratio of the Fourier transforms of the pressure and the normal velocity at the impedance boundary:
Z(ω) = P (ω)
Un(ω), (7)
with ω the angular frequency. The impedance is a complex number and can be split up in a real and imaginary part Z(ω) = R(ω) + iX(ω). The real part R(ω) is called the resistance and reflects the loss of acoustic energy from the system. The imaginary part X(ω), which is called the reactance, reflects how the kinetic energy is converted to potential energy.
3.1 Frequency vs. time domain
As can be seen in equation (7) the impedance Z is defined in the frequency domain and so are the pressure P and normal velocity Un. However, our goal is to implement an impedance surface in the CAA method which works in the time domain. Therefore the impedance needs to be trans- formed to the time domain. This will be done with the inverse Fourier transform. The Fourier transform F (ω) is defined as
F (ω) = Z ∞
−∞
f (t) e−iωtdt
so that the inverse Fourier transform f (t) is defined as f (t) = 1
2π Z ∞
−∞
F (ω) eiωtdω.
The impedance in the frequency domain can be written as P (ω) = Z(ω)Un(ω) and the con- volution integral is needed to transform the relation to the time-domain. The impedance in the time-domain is now given by the relation
p0(t) = Z ∞
−∞
z(τ ) (u0(t − τ ) · n)dτ
where u0(t) · n is the normal velocity in the time domain and will be written as u0n(t) from now on such that the previous equation becomes
p0(t) = Z ∞
−∞
z(τ ) u0n(t − τ )dτ. (8)
3.2 Material dependence
The impedance is dependent on which medium the air flows through. Different types of porous material have various acoustical properties that reflect how easy it is for air to flow through the material. The most important property is the flow resistivity σ which is the resistance that the rigid frame offers an air flow through the material. This can be described as the drop in pressure per unit length in the direction of the flow of air at unit speed (Ref. 1). The unit of flow resistiv- ity is Pa sm2. When the flow resistivity is high, like in asphalt, it is hard for air and thus for sound to travel through the ground. When the flow resistivity is small it is very easy for sound to travel through the ground, such as is the case with snow.
A difference is made between the flow resistivity σ of a material and the effective flow resistivity σe. The effective flow resistivity takes into account that the pores in the material sometimes have strange structures and that it is not possible for air to flow through all of the material. The struc- ture of the pore shapes can be described as the tortuosity q2which is the curviness of the pores (Ref. 1). When the pores in the porous material are straight in the flow direction, the tortuosity is equal to 1, when they have curves, lateral cavities or are in a different direction, the tortuosity is larger than 1. Another material property is the porosity Ω which is the volume fraction of air in the material (Ref. 14) and therefore lies between one and zero. Wilson et al. (Ref. 16) propose the following relation for the effective flow resistivity σe
σe = σΩ
q2 (9)
which is smaller than the flow resistivity σ for porous materials because of the values of the tor- tuosity and porosity.
The factorqΩ2 represents the air in the main pores per unit volume that gets moved by the pres- sure gradient. In lateral cavities next to main pores the air does not move because the pressure gradient is almost zero there, so not all the air in the main pores gets moved (Ref. 17). Some impedance models use just the effective flow resistivity σeand others use the flow resistivity σ in combination with other parameters such as the porosity Ω and tortuosity q2in the models.
When impedance is mentioned in the rest of this report, the characteristic impedance is meant.
This means it is assumed that the material properties do not depend on the frequency of the inci- dent acoustic waves.
3.3 Hard-backed layered ground
When calculating the impedance of a ground it is assumed that the ground is semi-infinite. Most of the ground types that are considered however have a layered structure, like grass-covered and snow-covered fields. Here the lower layer is acoustically hard, which means totally reflect- ing. For this case Rasmussen (Ref. 10) derived a formula for the impedance of a so called hard- backed layer
Z(d, ω) = Z(ω)ˆ
ρ0c0 coth(−ikd) (10)
where k is the complex wave number in the porous layer, d the thickness of the hard-backed layer and ˆZ(ω) the impedance of semi-infinite ground.
3.4 Physical validity
Impedance is a physical property of a porous ground surface, therefore it must be physically valid. For an impedance to be physically valid the next three physical conditions have to hold (Ref. 12).
Because the pressure p0(t) and the normal velocity u0n(t) are real, the impedance in the time- domain z(t) also has to be real. The consequence for the impedance in the frequency domain Z(ω) describes the first physical condition, the reality condition
Z(−ω) = Z∗(ω) for all ω ∈ R (11)
where Z∗ denotes the complex conjugate of Z. Note that this condition holds when R is even and X is odd in ω.
The second condition has to do with the energy in an acoustic flow. An impedance surface ab- sorbs part of the acoustic energy and therefore cannot add energy to the acoustic flow. As men- tioned before, the real part of the impedance represents a loss of acoustic energy if and only if it is positive. This means that for the impedance surface to not add energy, the real part of the impedance needs to be positive which leads to the passivity condition
R(ω) ≥ 0 for all ω ∈ R. (12)
The last condition for physical validity is the causality condition. This condition is based on the requirement that the pressure p0(t) cannot depend on the normal velocity u0n(t) of the future.
When looking at the relation between the pressure, velocity and impedance in the time-domain (8), note that for p0(t) to not depend on u0n(t) of the future, z(τ ) needs to be equal to zero for τ < 0. The result of this is that
Z(ω) is analytic in Im(ω) < 0.
Because it also holds that the velocity u0n(t) cannot depend on the pressure p0(t) of the future the condition has to be expanded such that the impedance is non-zero in Im(ω) < 0. Therefore the causality condition states that
Z(ω) is analytic and non-zero in Im(ω) < 0. (13)
A consequence of this condition is that the real and imaginary parts of the impedance are related through the Hilbert transform H[.] through X(ω) = H[R(ω)].
3.5 Impedance models
In the last few decades several models to describe an impedance were devised. Some models only use one parameter to describe the ground material like the Delany-Bazley and Miki models.
Others make use of more parameters, such as the phenomenological Zwikker-Kosten model, the two-parameter Attenborough model and Wilson’s relaxation model.
3.5.1 Zwikker-Kosten
A way to describe the impedance Z is by the phenomenological Zwikker-Kosten model (ZK model). This model was analytically derived by Zwikker and Kosten in 1949 by using linear acoustic equations (the momentum and continuity equation) for porous materials and a harmonic plane sound wave that satisfies the wave equation. The model is defined as
Z = Z∞
r1 + iωτ
iωτ (14)
with Z∞ = ρ0Ωc0q and τ = ρσΩ0q2. Here the three parameters of the model are the tortuosity q2, the flow resistivity σ and the porosity Ω.
Below the analytical derivation for the one dimensional case is given. The linear acoustic equa- tions for free air are
−∂u0
∂x = 1 ρ0c20
∂p0
∂t (15)
−∂p0
∂x = ρ0∂u0
∂t (16)
and the linear acoustic equations for porous materials are
−∂u0
∂x = Ω ρ0c20
∂p0
∂t (17)
−∂p0
∂x = ρ0
q2 Ω
∂u0
∂t + σu0. (18)
Here the first equation differs from the continuity equation for free air (15) through Ω. This is because the available volume of air is a factor Ω smaller than in free air. The second equation differs from the momentum equation for free air (16) in two ways. One of them uses Darcy’s law to add a term that deals with the use of porous materials through the flow resistivity σ (Ref.
9). The flow resistivity σ is defined as the measure of the drop in pressure per unit length in the direction of the flow of air at unit speed. In formula form this can be written as
σ = −∂p∂x0
u0 (19)
and can be found in the momentum equation for porous materials. The other difference comes from the lack of flow in some parts of the porous material. In section 3.2 this was described with the term qΩ2 which is also found in the momentum equation for porous materials.
From the linear acoustic equations for porous materials the impedance Z will be derived. Take a harmonic plane sound wave travelling in the positive x-direction where the pressure and velocity are given as
p0(x, t) = 1
2q(x, k) eiωt+1
2q∗(x, k) e−iωt (20)
with
q(x, k) = A e−ikx, (21)
and
u0(x, t) =1
2v(x, k) eiωt+1
2v∗(x, k) e−iωt (22)
with
v(x, k) = B e−ikx. (23)
Here A and B are constants and the superscript ∗ denotes the complex conjugate. Then the impedance is derived as Z = v(x,k)q(x,k) = AB. When substituting the first (or second) term for the pressure and velocity into the continuity equation (17), an equation for the complex wave number k expressed in the impedance Z is gained
−∂u0
∂x = Ω
ρ0c20
∂p0
∂t
⇒ −(−ik) 1
2B e−ikxeiωt
= Ω
ρ0c20(iω) 1
2A e−ikxeiωt
⇒ (ik)B = Ω ρ0c20(iω)A
⇒ k = ωΩ ρ0c20
A
B = ωΩ
ρ0c20Z (24)
When substituting the pressure and velocity into the momentum equation (18), an equation for the impedance Z = AB expressed in the complex wave number k is given as
−∂p0
∂x = q2ρ0
Ω
∂u0
∂t + σu0
⇒ −(−ik) 1
2A e−ikxeiωt
= q2ρ0
Ω (iω) 1
2B e−ikx eiωt
+ σ 1
2B e−ikxeiωt
⇒ (ik)A = iρ0q2ω
Ω + σ
B
⇒ Z = A
B = ρ0q2ω
Ω − iσ 1
k (25)
Now combining the complex wave number (24) and the impedance (25) gives Z = ρ0q2ω
Ω − iσ 1
k = ρ0q2ω
Ω − iσ ρ0c20 ωΩ
1 Z
⇒ Z2 = ρ20c20q2
Ω2 − iρ0c20σ ωΩ
⇒ Z = (ρ0c0) s
q2
Ω2 − i σ
ρ0ωΩ = ρ0c0q Ω
r 1 − i 1
ωτ = Z∞
r1 + iωτ iωτ with Z∞= ρ0Ωc0q and τ = ρσΩ0q2. This is exactly the Zwikker-Kosten model (14).
When using the second term of the velocity and pressure, the impedance becomes its complex conjugate and the reality condition (11) follows. The passivity condition (12) also holds because Z∞and τ are positive and so is the real part of the root. The root of the impedance can be writ- ten as
r1 + iωτ
iωτ =
r 1 − i 1
ωτ
= v u u t
q
1 +ω21τ2 + 1
2 − i
v u u t
q
1 +ω21τ2 − 1
2 ,
where the last expression is just the root of a complex number split into a real and imaginary part. From this it can be seen that the real part of the root is indeed bigger than zero and the pas- sivity condition holds. The causality condition (13) is also satisfied because
q1+iωτ
iωτ is non-zero for Im(ω) < 0 and also analytic, because the pole ω = 0 does not lie in the lower-half plane.
Therefore this description of the impedance is physically valid.
3.5.2 Delany-Bazley/Miki
One of the most popular and most-used models is the empirical one-parameter model that was derived from experimental data by M.E. Delany and E.N. Bazley in 1969 (Ref. 3). Their model for the impedance has the following form
Z = ρ0c0 1 + a ρ0ω 2πσe
b
+ ic ρ0ω 2πσe
d!
(26) where ρ0and c0are respectively the atmospheric pressure and speed of sound. In this report they are taken as ρ0 = 1.225mkg3 and c0 = 340ms. Delany and Bazley experimentally derived the following coefficients for their model (26)
a = 0.0595 b = −0.75 c = −0.0891 d = −0.73.
Here ρ2πσ0ω is dimensionless and thus the impedance has the same dimension as ρ0c0, which is
kg
m2s. This model only has one parameter, σe, while the modified ZK model had three, σ, q and Ω. Delany and Bazley based their model on experimental results of various fibrous absorbing materials like glass-fibre and mineral-wool materials and assumed that these materials had porosi- ties Ω and tortuosities q2close to 1. In that case σeis replaced by σ in (26).
Usually absorbing materials do not have porosities and tortuosities close to 1 and therefore an adapted parameter is used in the DB model, the effective flow resistivity σe. The effective flow resistivity allows the porosity and tortuosity of a material to differ from 1 by the following rela- tion for the effective flow resistivity as shown in Section 3.2: σe= σΩq2 .
Delany and Bazley did not only derive a model for the impedance, but also for the complex wave number k
k = ω c0
1 + q ρ0ω 2πσe
r
+ ip ρ0ω 2πσe
s
(27) where the coefficients are
q = 0.0989 p = 0.1972 r = −0.70 s = −0.59.
The complex wave number is needed when a hard-backed layered ground (section 3.3) is used.
In 1990 Yasushi Miki (Ref. 8) stated that not all the physical properties from the previous section hold for the DB model because of the choice of coefficients. Delany and Bazley assumed that the frequency only has positive real values. When that is assumed, the model for the impedance (26) satisfies the reality condition (11). The passivity condition (12) is almost satisfied, as can be seen in figure 1 for a flow resistivity σ of 200 kPa sm2 . The causality condition (13) is however not satisfied.
Fig. 1 Relative impedance ρZ
0c0 of the DB model for σ = 200 kPa sm2
Miki derived two relations between the coefficients which should hold for the model to be phys- ically valid. He derived these from the consequence of the causality conditions (13). The con- sequence says that X(ω) should be the Hilbert transform of R(ω). To require relations for the coefficients from X(ω) = H[R(ω)], where H[.] is the Hilbert transform, the following relation is used
H[|x|ν−1] = −cotνπ
2 sgn(x)|x|ν−1. (28)
When applying this relation to the DB model (26) and assuming that ω is positive, the following two relations for the coefficients are derived
d = b and c = −a cot(b + 1)π
2 . (29)
Here it was assumed that the reality condition (11) holds and the passivity condition (12) holds when a ≥ 0. With these relations in mind, Miki used the experimental data of Delany and Bazley to come to the following coefficients:
a = 0.0794 b = −0.632 c = −0.1218 d = −0.632.
Miki also derived relations for the coefficients of the complex wave number of the DB model.
Miki expresses the complex wave number by a single complex power-law relation
¯k(¯ω) = ω¯ ic0
(1 + iafb).
When substituting the complex frequency ¯ω = iω in the above relation
¯k(iω) = iω
ic0(1 + iaibωb)
= ω c0
1 + ia cos(bπ
2 )ωb+ i2a sin(bπ 2 )ωb
= ω c0
1 − a sin(bπ
2 )ωb+ ia cos(bπ 2 )ωb
and comparing this to the equation for the complex wave number of Delany and Bazley (27) k = ω
c0
1 + q
ρ0ω 2πσ
r
+ ip
ρ0ω 2πσ
s
= ω c0
1 + q ρ0 2πσ
r
ωr+ ip ρ0 2πσ
s
ωs
, (30)
the assumption of Miki gives the following relations for the coefficients s = r
and
q ρ0 2πσ
r
= −a sin(rπ 2 )
= −
a cos(rπ
2 ) sin(rπ2 ) cos(rπ2 ).
From equation (30), a cos(rπ2 ) should be equal to p 2πσρ0 s
such that the previous equation be- comes
q
ρ0
2πσ
r
= −p
ρ0
2πσ
s
tan(rπ 2 )
= −p ρ0 2πσ
r
tan(rπ 2 ), which gives
q = −p tanrπ 2 .
From these relations and again the experimental data of Delany and Bazley, Miki derived the following coefficients of the complex wave number
q = 0.1810 p = 0.1239 r = −0.618 s = −0.618.
3.5.3 Other models
Next to these two models there are many more to be found in the literature. Some models make use of more parameters to describe the porous materials and others make use of different param- eters. Below two more examples are presented of models for the impedance of porous materials, others can be found in additional literature like the book by Attenborough (Ref. 1).
The Attenborough two-parameter model (Ref. 1) makes use of the effective flow resistivity σe
and the effective rate of change of porosity with depth αe. He assumes that porous materials do not have a stable porosity but they change in depth. The impedance is given by Attenborough as
Z = 0.436(1 + i)
r2πσe
ω + 19.74i2παe
ω . (31)
Another model was derived by Wilson et al. (Ref. 16). This so-called relaxation model has 4 parameters and is defined as
Z = ρ0c0q Ω
1 + γ − 1
√1 − iωτe
1 − 1
√1 − iωτv
−1/2
. (32)
τv and τeare the aerodynamic and thermodynamic characteristic times and are respectively given by Wilson as
τv = 2ρ0q2 Ωσ τe= NP Rs2Bτv
with NP R the Prandtl number. The parameters q2, Ω, σ, and sBappear in the expression for the aerodynamic and thermodynamic characteristic times. sBis the only parameter that has not been seen before and denotes the pore shape factor.
This model mimics the DB model in the range 0.01 < ρσ0f
e < 1, but is also valid outside this range when using the following approximations for the characteristic times
τv ≈ 2.1ρ0 σ τe≈ 3.1ρ0
σ .
3.6 Approximations
To calculate the pressure from the relation P (ω) = Z(ω)Un(ω) the inverse Fourier transforms of the impedance
z(t) = 1 2π
Z ∞
−∞
Z(ω) eiωtdω and the normal velocity
u0n(t) = 1 2π
Z ∞
−∞
Un(ω) eiωtdω
are needed. The pressure can be calculated with the convolution product (8), which can also be written as
p0(t) = Z ∞
−∞
z(t − τ ) · u0n(τ )dτ. (33)
However it is not straightforward to calculate the inverse Fourier transform of the impedance
models of the previous section. It is therefore useful to approximate the expression of the impedance.
There are several ways to do this.
3.6.1 Poles
One of the ways was described by Reymen et al. (Ref. 11) for a broadband model and is the sum of S first-order systems with real poles and T second-order systems with complex conjugate poles. The impedance can now be written as
Z(ω) =
S
X
k=1
Zk(ω) +
T
X
l=1
Zˆl(ω), (34)
where
Zk(ω) = Ak
λk+ iω = Akλk
λ2k+ ω2 − i Akω
λ2k+ ω2 (35)
and
Zˆl(ω) = Aˆl
µl+ iω + Bˆl
µ∗l + iω = Dl+ iωCl
(αl+ iω)2+ β2l. (36)
Here λk ≥ 0 is called a real pole and µland µ∗l are a pair of complex conjugate poles αl± iβl, with αl≥ 0.
The inverse Fourier transforms of these systems are easy to calculate and are
zk(t) = Ake−λktH(t) (37)
and ˆ
zl(t) = e−αlt
Clcos(βlt) +Dl− αlCl βl
sin(βlt)
H(t) (38)
where H(t) is the Heaviside function which is zero for t < 0 and 1 for t ≥ 0. Then the inverse Fourier transform of the poles approximation becomes
z(t) =
S
X
k=1
zk(t) +
T
X
l=1
ˆ
zl(t). (39)
Some restrictions on the parameters are needed for this approximation to be physically valid. In the second-order system the coefficients Cland Dlcan be written as Cl = ˆAl+ ˆBland Dl = αl( ˆAl+ ˆBl) + iβl( ˆBl− ˆAl). For the reality condition (11) to hold Ak, Cland Dlneed to be real, which results in the condition ˆBl = ˆA∗l. For the causality condition (13) to hold, the poles have to satisfy λk ≥ 0 and αl ≥ 0 and the model is passive (12) when Ak ≥ 0 and for a good choice of Cland Dl. Ak ≥ 0 is a very strong restriction, because the model can also be passive when not all the coefficients Akare bigger than or equal to zero.
In all our models the imaginary part of the impedance is negative, therefore the imaginary part of the second-order system needs to be negative. This results in constraints for ˆAland ˆBland subsequently for Cland Dl.
The coefficients of this approximation can be found by a fitting method which fits this approxi- mation to one of the models or a data set. Two fitting methods were given by Cott´e et al. (Ref. 2) One is the Vector Fitting technique proposed by Gustavsen and Semlyen (Ref. 5) and the other makes use of the MATLAB function fmincon, an -constraint method. The first method can give complex conjugate poles, but when keeping the number of poles relatively small (below 10), only real poles appear. The second method assumes that only real poles are being used.
3.6.2 MSD+
The same Reymen et al. described a special case of the pole approximation, which is an approx- imation for a single frequency. This approximation is better known as the mass-damper-spring approximation
Z(ω) = a−1
(iω) + a0+ a1(iω).
In the article of Heutschi et al. (Ref. 6) this approximation was extended with an extra term a−2(iω)−2and in that way was made second order in ω. It turns out that this MSD+ approxi- mation is a better approximation then the first order approximation of Reyman et al. because now
the real part of the impedance is described with two parameters instead of a single constant. The MSD+ approximation is written as
Z(ω) = a−2(iω)−2+ a−1(iω)−1+ a0+ a1(iω). (40)
The inverse Fourier transform of the MSD+ approximation can be found by using the following relations
∂y(t)
∂t
←→ (iω)Y (ω)F
Z t
−∞
y(τ )dτ ←→ (iω)F −1Y (ω).
The inverse Fourier transform then becomes p0(t) = a−2
Z t
−∞
Z t0
−∞
u0n(τ ) dτ dt0 + a−1
Z t
−∞
u0n(τ ) dτ + a0u0n(t) + a1
∂u0n(t)
∂t . (41)
Also for this approximation it needs to be checked which restrictions are needed for all the phys- ical conditions to hold. The reality condition (11) follows directly from the approximation when all the parameters are real. The passivity condition (12) is easily maintained when setting a−2 ≤ 0 and a0 ≥ 0. The causality condition (13) requires a bit more effort. The expression for the impedance can be rewritten as a rational function
Z(ω) = a−2+ a−1(iω) + a0(iω)2+ a1(iω)3
(iω)2 .
The frequency can now be written as a real and imaginary part as ω = ωR + iωI. For the impedance to be causal it has to be analytic and non-zero in Im(ω) = ωI< 0. For the impedance to be analytic in ωI < 0 the poles have to lie in the upper-half plane. It can easily be seen from the previous equation that the poles of the impedance are located in 0, which is assumed to lie in the upper-half plane. Also needed for causality is the impedance to be non-zero in ωI < 0. This means that a−2 + a−1(iω) + a0(iω)2 + a1(iω)3 cannot be equal to zero. Already there are two constraints on the parameters: a−2 ≤ 0 and a0 ≥ 0. Now the other two parameters a−1and a1 need to be chosen in such a way that the zeroes lie in the upper-half plane.
Heutschi et al. suggest that finding the coefficients can be done by fitting the MSD+ approxima- tion to a model by minimizing the error of the squared differences of the model and the MSD+
approximation.
3.7 Conclusions
This chapter gave an overview of some different ground impedance models and ways to approxi- mate these models. For the next part of the report, the numerical discretization of the approxima- tions will be used and therefore a choice of model and approximation will be made.
First let us discuss the different models that were presented. These models were the Zwikker- Kosten model (14) and the Delany-Bazley model (26) with two choices for the coefficients. The ZK model is based on a mathematical derivation, but has more parameters. The DB model has a simple form, and with the Miki coefficients the model is also physically valid. In figure 2 the three different models for a range of frequencies and flow resistivity of σ = 200 kPa sm2 , porosity of Ω = 0.5 and the tortuosity of q2 = 1.3 were plotted. In the Miki and Delany-Bazley method the effective flow resistivity σeis used as parameter. The effective flow resistivity depends on the tortuosity and porosity and in this case is equal to σe≈ 76.92 kPa sm2
(a) Real part of the relative impedance for σ = 200kPa sm2 , Ω = 0.5 and q2= 1.3
(b) Imaginary part of the relative impedance for σ = 200kPa sm2, Ω = 0.5 and q2= 1.3
Fig. 2 Relative impedance for the DB, Miki and ZK models
As can be seen, the models do not differ that much. Therefore the choice of model is mainly one of convenience. For this reason the DB model with Miki coefficients is chosen, because it only has one parameter to work with and is a physically valid model.
Which approximation to use is another choice that has to be made. The choice is between the pole approximation (34), which is a broadband model, and the MSD+ approximation (40), which is a small-band model. Both models are easily transferred to the time-domain and can both be made physically valid by the right choice of coefficients. The discretization of both approxima- tions has not yet been discussed, but also plays a role in the decision which approximation to use.
In sections 4.4.3 and 4.4.4 two discretizations for the pole approximation are discussed and two
discretization for the MSD+ approximation can be found in appendix B. As can be seen in ap- pendix B, all the possible discretization methods for the MSD+ approximation have some fault.
From instability issues, needing to approximate the derivative of the normal velocity to division by zero. A discretization method of the pole approximation also needs approximations, but one method does seem to work well without any concerns. Therefore the pole approximation is pre- ferred.
Also important is the broad usability of the approximations. Because the pole approximation can be fitted for a large band of frequencies, our preference goes to this approximation. However only real poles for the approximation are considered such that the pole approximation is written as
Z(ω) =
S
X
k=1
Ak
λk+ iω =
S
X
k=1
Akλk
λ2k+ ω2 − i Akω
λ2k+ ω2. (42)
4
Numerical method
To acquire solutions for an acoustical problem with an impedance boundary, the problem has to be transformed to a computational problem. In the first section the two-dimensional discretized form of the linearized Euler equation is described which can be adapted for one-dimensional and higher-dimensional problems. In the next two sections the numerical methods for the time integration and spatial discretization as used by NLR are discussed. In section 4.4 two ways to discretize an impedance boundary condition are discussed. The two ways are then put to the test in a one-dimensional test case in terms of the accuracy and stability of the discretized impedance boundary conditions.
4.1 Computational problem
Assuming a two-dimensional domain D ∈ R2, the linearized Euler equations (4) prescribe the acoustic flow on this domain. Before the acoustical problem can be solved numerically it has to be defined as a computational problem.
As seen in section 2.1, the flow can be split into the mean flow and the perturbations of the mean flow (3). The computational problem considers a prescribed mean flow with velocity u0 = (u0, v0). Also prescribed are the density ρ0of the mean flow and the speed of sound c0. Because it is assumed that the domain only contains air, which is assumed to be a calorically perfect gas, the pressure p0of the mean flow can be calculated from
c0 =r γp0
ρ0 ⇒ p0 = c20ρ0 γ .
Here γ is the heat capacity ratio which is 1.4 for air. The energy per unit mass of the mean flow for a calorically perfect gas E0is calculated from the relation (2) as
E0 = p0
(γ − 1)ρ0 +1
2(u20+ v02).
The mean flow variables ρ0, p0, E0, c0 and u0 are prescribed, but not necessarily constant through- out the domain. Another mean flow variable that can be seen in the linearized Euler equations (4) is the total enthalphy per unit mass H0, which can be calculated from the energy of the mean flow E0by
H0= E0+ p0 ρ0.
Now having all the mean flow variables prescribed, the linearized Euler equations (4) in conser- vative form of the 2D-domain D can be rephrased as
∂U
∂t = ∇ · F (U ) (43)
with U the flow state vector
U =
ρ0 (ρux)0 (ρuy)0 (ρE)0
.
Here ∇ denotes the two-dimensional gradient and F (U ) is defined as
F (U ) =
u0ρ0+ ρ0u0x v0ρ0+ ρ0u0y u0(ρux)0+ (ρ0u0)u0x+ p0 v0(ρux)0+ (ρ0u0)u0y
u0(ρuy)0+ (ρ0v0)u0x v0(ρuy)0+ (ρ0v0)u0y+ p0 u0(ρH)0+ (ρ0H0)u0x v0(ρH)0+ (ρ0H0)u0y
. (44)
When integrating the conservative linearized Euler equations (43) over a fixed control volume Ω ⊂ D and applying Gauss’ divergence theorem the integral form of the equations is obtained:
Z
Ω
∂U
∂t dV + Z
δΩ
F (U ) · n dA = 0, (45)
where δΩ is the boundary of the fixed volume Ω and n is the unit normal pointing outside of the volume. The entire domain D can be discretized into grid cells Ωi,j. When applying the integral form of the equations to a grid cell, the discretized form of the integral form of the linearized Euler equations (45) is written as:
Vi,j
∂U
∂t
i,j
+ Bi,j = 0 (46)
where Vi,jis the volume of the grid cell Ωi,j: Vi,j =
Z
Ωi,j
dV
and Bi,j is the discretized form of Z
δΩi,j
F (U ) · n dA
which is called the flux balance. In section 4.2 time-integration methods are discussed and in section (4.3) the derivation of Bi,jis discussed.
4.2 Time integration
In this section the time integration of the discretized linearized Euler equations (46) is discussed.
Reorganizing the discretized equations (46) gives
∂U
∂t
i,j
= −Bi,j
Vi,j = f (t, U ).
There are various methods to solve this differential equation and a well-known and used range of solutions methods are the Runge-Kutta methods. The general structure of a Runge-Kutta method is
Un+1= Un+ ∆t
s
X
i=1
bif (tn+ ci∆t, Yi) (47)
where
Yi = Un+ ∆t
s
X
j=1
aijf (tn+ cj∆t, Yj) for i = 1, ..., s. (48) A special case of the Runge-Kutta method is used by NLR to discretize the time-integration.
This is the standard fourth-order low-storage Runge-Kutta method and is given as Un+1= Y4
and
Yi = Un+ γi∆t f (tn+ γi−1∆t, Yi−1) for i = 1, ..., 4 (49) with coefficients given as
γ0 = 0, γ1= 1
4, γ2 = 1
3, γ3 = 1
2, γ4 = 1. (50)
Each step of the fourth-order low-storage Runge-Kutta is linear with a different time step de- pending on the coefficients, but together they result in a fourth-order method for linear equations.
The method is called low-storage because only the previous stage has to be stored to calculate the next one. The stability region of the standard fourth-order low-storage Runge-Kutta is the same as that of the classical fourth-order Runge-Kutta method.
4.3 Spatial schemes and artificial diffusion
The second part of the discretized equations (46) that needs to be discussed is the spatial scheme.
Consider a one-dimensional grid which is uniformly divided into grid cells of length ∆x. On this grid a function f (x) is defined. A finite-difference approximation of the first derivative of the function f (x) at point x is given as
∂f
∂x(x) ≈ 1
∆x
M
X
j=−N
ajf (x + j∆x) (51)
where the coefficients aj determine the kind of approximation and the order of accuracy. Exam- ples of finite difference approximations are the forward difference approximation of order 1
∂f
∂x(x) ≈ 1
∆x(f (x + ∆x) − f (x)) and central difference approximation of order 2
∂f
∂x(x) ≈ 1
∆x
1
2f (x + ∆x) − 1
2f (x − ∆x)
.
The coefficients of these and other examples and their order of accuracy are found in table 1.
Accuracy a−3 a−2 a−1 a0 a1 a2 a3
1 -1 1
2 -12 0 12
4 121 -23 0 23 -121
6 -601 203 -43 0 34 -203 601
Table 1 The coefficients of some finite difference approximations for a first-order derivative and their order of accuracy
The finite-difference approximation that is used by NLR is the central 7-points dispersion-relation- preserving (DRP) scheme of Tam and Webb (Refs. 15, 7 and further in this section). The coeffi- cients for this approximation are
a0= 0
a1= −a−1= 0.770882380518 a2= −a−2= −0.166705904415 a3= −a−3= 0.0208431427703.
The coefficients are chosen in such a way that the approximation is fourth-order accurate and low dispersive. The low dispersiveness of the method means that the effective numerical wave num- ber approximates the real wave number of the solution for waves longer than 8 mesh spacings (Ref. 13).
