Determination of the Acoustic Power Generation Potentiality of One-Dimensional Acoustic Elements 

Technische Universit¨at M¨unchen Prof. Dr.-Ing. Thomas Sattelmayer Prof. Wolfgang Polifke Ph.D. (CCNY)

– Bachelor Thesis –

Determination of the Acoustic Power Generation Potentiality of One-Dimensional

Acoustic Elements

Author:

Alexander Avdonin Registration Number:

3601787 Supervisor:

Prof. Wolfgang Polifke Dipl.-Ing. Tobias Holzinger

4. Juli 2011

Lehrstuhl f¨ur Thermodynamik, Technische Universit¨at M¨unchen

Best¨ atigung der eigenst¨ andigen Arbeit

Ich erkl¨are hiermit, dass ich diese Arbeit ohne fremde Hilfe angefertigt und nur die in dem Literaturverzeichnis angef¨uhrten Quellen und Hilfsmittel benutzt habe.

Garching, den 4.07.2011

Name: Alexander Avdonin

Matrikelnummer 3601787

Abstract

In this bachelor thesis it is shown, how to compute the minimal and maximal potentiality of the generated acoustic power from the scattering or transfer matrix of a linear one-dimensional acoustic element in presence of mean flow. The function for the computation of the ratio of the scattered acoustic power to the incident acoustic power is derived from the scattering or transfer matrix, the amplitudes of the incident Riemann invariants and their phase difference.

By a small modification the derived function can calculate the acoustic power ratio from the Riemann invariants at the inlet of an element. The functions are implemented in a computation tool for MATLAB, which is used to investigate some acoustic systems. The obtained results are compared to the analytical solutions.

Contents

Contents VI

List of figures VII

List of tables VIII

Notation IX

1. Introduction 1

2. Characterisation of Acoustic Elements 2

2.1. Thermodynamical Variables in Linear Acoustics . . . 2

2.2. Matrix Representation of Acoustic Systems . . . 4

3. Acoustic Power Dissipation or Generation 6 3.1. Acoustic Power Dissipation in One-Dimensional Ducts in Presence of Uniform Flow . . . 6

3.2. The Scattering Matrix for the Acoustic Power Amplitudes and the Potentiality 9 3.3. The Ratio of the Scattered Acoustic Power to the Incident Acoustic Power . . 10

4. Description of the MATLAB tool 13 4.1. potentiality computation.m . . . 14

4.2. potentiality plot.m . . . 16

4.3. find data index.m . . . 17

4.4. pgen plot f1g2.m . . . 17

4.5. pgen fmincon f1g2.m . . . 19

4.6. pgen plot f1g1.m . . . 20

5. Illustrative Examples 21 5.1. Uniform Tube with Non-Viscous Adiabatic Mean Flow . . . 21

5.2. Uniform Tube with Mean Flow and Thermo-Viscous Loss . . . 22

5.3. Sudden Area Change with Incompressible Non-Viscous Flow . . . 27

5.4. Investigation of a Complex Acoustic System . . . 33

6. Conclusion 41

A. Appendix 42 A.1. Computation of the Transfer Matrix for the Riemann Invariants . . . 42 A.2. Computation of the Scattering Matrix for the Riemann Invariants . . . 42 A.3. Computation of the Scattering Matrix for the Acoustic Power Amplitudes . . 43 A.4. How to use the MATLAB Computation tool . . . 44 A.5. Computation Parameters of a Complex Acoustic System . . . 47

Bibliography 49

List of Figures

2.1. The Riemann invariants, the incident and scattered quantities of an acoustic

element. . . 3

2.2. The upstream and downstream quantities of an acoustic element. . . 5

3.1. The incident and scattered acoustic power amplitudes of an acoustic element. . 9

4.1. The structure arrayinput struct. . . 15

4.2. The structure arrayoutput struct . . . 15

4.3. Potentiality plot created by potentiality plot.m. . . 16

4.4. A polar plot created by pgen plot f1g2.m. . . 18

4.5. A Cartesian plot created by pgen plot f1g2.m. . . 18

5.1. The potentialities of the generated acoustic power of a uniform tube with mean flow and thermo-viscous loss, M = 0.05 . . . 25

5.2. The potentialities of the generated acoustic power of a uniform tube with mean flow and thermo-viscous loss, M = 0.1 . . . 25

5.3. η_si of a tube with mean flow and thermo-viscous loss. . . 26

5.4. The potentialities versus the cross-sectional area ratio. . . 30

5.5. The dominating potentialities versus the Mach number at the inlet with two different cross-sectional area ratios. . . 31

5.6. η_si of an expanding element. . . 31

5.7. η_si of a contracting element. . . 32

5.8. The virtual construction of a composed acoustic system. . . 33

5.9. The potentialities of a complex system in comparison to a flame element. . . . 35

5.10. The potentialities of the complex system in comparison to the product of the potentialities of each separate element of this system. . . 36

5.11. The potentialities of the complex system with different lengths of the tube at the middle of the compound. . . 37

5.12. The potentialities of the complex system with different α₁₂. . . 38

5.13. The potentialities of the composed system with differentα45. . . 39

5.14. The ratio of the scattered acoustic power to the incident acoustic power with the default colormap for the complex system. . . 40

5.15. The ratio of the scattered acoustic power to the incident acoustic power with the custom colormap for the complex system. . . 40

List of Tables

4.1. Overview of the developed MATLAB tool. . . 13 4.2. The input transfer or scattering matrix . . . 14 4.3. Optimization options . . . 20 5.1. Computation parameters of a uniform tube with mean flow and thermo-viscous

loss. . . 24 5.2. Computation parameters of a sudden area change with non-viscous mean flow. 29

Notation

Latin characters

A cross-sectional area or surface [m²]

M Mach number

I acoustic intensity [W/m²]

P acoustic power [W]

R specific gas constant [J/(kgK)]

S scattering matrix

T transfer matrix

a amplitude of incident acoustic power [W^1/2] b amplitude of scattered acoustic power [W^1/2]

c speed of sound [m/s]

f frequency of sound wave [Hz]

f, g Riemann invariants

k wave number [m⁻¹]

l length [m]

m⁰ acoustic mass flow density [kg/(m²s)]

p pressure [Pa]

u velocity [m/s]

Greek characters

α cross-sectional area ratio

γ adiabatic index

⁰ specific acoustic exergy [W/kg]

η ratio of acoustic power

λ potentiality of the generated acoustic power ξ potentiality of the dissipated acoustic power

ρ density [kg/m³]

ϕ phase difference between the incident Riemann invariants [rad]

φ phase difference between the Riemann invariants at the inlet [rad]

ω angular frequency [rad/s]

Sub- or superscripts, other symbols

q₁ quantity at inlet of acoustic element q₂ quantity at outlet of acoustic element

q^t transpose of quantity

q^∗ complex conjugate of quantity

q^† complex conjugate and transpose of quantity q⁺ quantity of incident wave

q⁻ quantity of scattered wave

< q > time average of quantity

||q|| norm of vector q q_tot total value of quantity q⁰ fluctuating value of quantity

1. Introduction

This bachelor thesis is based on the work of Y. Aur´egan and R. Starobinski [1]. They presented a method to compute the minimal and maximal ratio of the dissipated or generated acoustic power to the incident acoustic power from the scattering matrix of an acoustic element. The aim of the thesis is to implement and verify this method and additional functionality in a MATLAB tool.

In the next chapter the definitions of used physical quantities are presented to avoid any kind of misunderstanding. The third chapter explains the above-mentioned method. It also shows, how to compute the ratio of the scattered acoustic power to the incident acoustic power from the amplitudes of the incident Riemann invariants and their phase difference. Within a small modification the acoustic power ratio is also expressed by the Riemann invariants at the inlet of an element. The fourth chapter describes the developed MATLAB tool, which contains the derived computation prescriptions. In the fifth chapter the computation tool is applied on the following acoustic elements:

•uniform tube with non-viscous adiabatic flow,

•uniform tube with flow and thermo-viscous loss,

•sudden area change with incompressible non-viscous flow.

The obtained results are compared to the analytical solutions. At the end, a less trivial acoustic system is investigated, which contains several acoustic elements.

The developed tool allows to obtain the potentialities of the generated acoustic power, which represent the limits of the acoustic power generation. Therefore the general behaviour of the system is determined. The manipulation of the entries of the transfer or scattering matrix serves to achieve the required general behaviour. The consideration of the ratio of the scattered acoustic power to the incident acoustic power enables to find the critical ratio the Riemann invariants, for which the system shows undesired behaviour, e.g. becomes unstable.

The obtained knowledge can be used in design of an acoustic system to avoid this critical ratio.

2. Characterisation of Acoustic Elements

2.1. Thermodynamical Variables in Linear Acoustics

In acoustics flow variables consist of a mean value and a superimposed fluctuation. The total value is denoted by an index _tot and the acoustic fluctuation with an apostrophe ’, for example ptot = p+p⁰. For notation convenience from now on this convention is used for pressure, velocity and density.

In linear acoustics the pressure and velocity fluctuations can be expressed in terms of the Riemann invariants¹f and g:

p⁰

ρc =f +g u⁰ =f −g

(2.1)

wheref represents the downstream propagation and g the upstream propagation.Besides the Riemann invariants, there is another common notation with”⁺” and ”⁻”, for example

p⁰ =p⁺+p⁻,

wherep⁺ is pressure of an incident wave andp⁻ the pressure of a scattered wave. Notice that the third notation also exists, where”⁺”denotes the direction of the mean flow and ”⁻” the opposite direction. In this paper the first notation is used:

”⁺” denotes the quantities of the incident (incoming) waves,

”⁻” denotes the quantities of the scattered (outgoing) waves.

The relations between incident and scattered quantities and the Riemann invariants can be derived by considering an one-dimensional acoustic element (see fig. 2.1). At the inlet of the element the relations are

p⁺₁ +p⁻₁

ρ₁c₁ =f₁+g₁ =⇒ p⁺₁

ρ₁c₁ =f₁ and p⁻₁

ρ₁c₁ =g₁ (2.2)

u⁺₁ +u⁻₁ =f₁−g₁ =⇒u⁺₁ =f₁ and u⁻₁ =−g₁. (2.3)

1The derivation of the Riemann invariants can be found e.g. in [3].

Figure 2.1.: The upper illustration shows the Riemann invariants of an acoustic element. Ac- cording to the mean flow direction,f represents the downstream propagation and g the upstream propagation.The lower illustration shows the incident and scat- tered quantities of an acoustic element. The quantities with ”⁺” belong to the incident wave and the quantities with ”⁻” to the scattered wave. The inlet of an element is denoted with index 1 and the outlet with index2.

Here the downstream propagating wave corresponds to the incoming wave and the upstream propagating wave to the scattered wave.

At the outlet the considered relations are different:

p⁺₂ +p⁻₂ ρ2c2

=f₂+g₂ =⇒ p⁺₂ ρ2c2

=g₂ and p⁻₂ ρ2c2

=f₂ (2.4)

u⁺₂ +u⁻₂ =f₂−g₂ =⇒u⁺₂ =−g₂ and u⁻₂ =f₂. (2.5) These relations show that fluctuating quantities of the incident and scattered waves are also related:

p⁺₁ ρ1c1

=u⁺₁ (2.6)

p⁻₁

ρ₁c₁ =−u⁻₁ (2.7)

p⁺₂

ρ₂c₂ =−u⁺₂ (2.8)

p⁻₂

ρ₂c₂ =u⁻₂. (2.9)

There is one more basic equation that should be mentioned here. The pressure fluc- tuation is a function of the density fluctuation. According to [3], its first-order isentropic approximation² is

p⁰ =c²ρ⁰. (2.10)

It also holds for the incident or scattered pressure fluctuation.

In section 3.1 an acoustic power is expressed by two variables: an acoustic mass flow density m⁰ and a specific acoustic exergy ⁰. They are defined by Y. Aur´egan and R. Starobinski [1] as

m⁰ =ρu⁰ +ρ⁰u (2.11)

⁰ = p⁰

ρ +uu⁰. (2.12)

These quantities also can be split up into incident and scattered parts and then expressed by the Riemann invariants at the inlet or at the outlet, for example

⁺₁ = p⁺₁

ρ₁ +u1u⁺₁ =f1(c1+u1) ⁻₁ = p⁻₁

ρ₁ +u₁u⁻₁ =g₁(c₁ −u₁) ⁺₂ = p⁺₂

ρ₂ +u₂u⁺₂ =g₂(c₂ −u₂) ⁻₂ = p⁻₂

ρ2

+u₂u⁻₂ =f₂(c₂+u₂).

(2.13)

2.2. Matrix Representation of Acoustic Systems

The Matrix representation is ideally suited for the analysis of one-dimensional linear acoustic systems. There are different combinations of input and output quantities. In this work three of them are considered:

a) A transfer matrix for the pressure and velocity fluctuations T_pu is defined by





p2

ρc

u





2

=T_pu





p ρc

u





1

. (2.14)

This matrix relates the upstream and downstream pressure and velocity fluctuations (see Figure 2.2).

b) A transfer matrix for the Riemann invariantsTf g is defined by



 f g





2

=Tf g



 f g





1

. (2.15)

2For linear acoustics is the first-order approximation fully justified.

Figure 2.2.: Illustration of the upstream and downstream quantities of an acoustic element.

According to the mean flow direction, the upstream quantities are on the left side of an acoustic element and the downstream quantities on the right side.

It describes the correlation between the upstream and the downstream Riemann invari- ants (see fig. 2.1).

c) A scattering matrix for the Riemann invariantsS_{f g} connects the incident (outgoing) and the scattered (incoming) Riemann invariants (see fig. 2.1)



 g₁ f₂



=S_{f g}



 f₁ g₂



. (2.16)

3. Acoustic Power Dissipation or Generation

3.1. Acoustic Power Dissipation in One-Dimensional Ducts in Presence of Uniform Flow

In this section the dissipated acoustic power is derived for a linear one-dimensional acoustic element¹ in presence of uniform flow. Furthermore, only plane waves are considered. Waves with higher modes are evanescent, if their frequencies are smaller than the cut-off frequency for the first mode².

A generalized definition of an acoustic intensity, which preserves the continuity propagation of sound power in any isentropic and irrotational flow is given by C. L. Morfey [2]

I =< p⁰u⁰ >+u²

ρc² < p⁰² >+u²

c² < p⁰u⁰ >+ρu < u⁰² >, (3.1) where<> indicates the time average. This equation can also be formulated with the specific acoustic exergy η⁰ and the mass flow density m⁰

I =<(ρu⁰+ρ⁰u) p⁰

ρ +uu⁰

>=< ⁰m⁰ >. (3.2) Sound power, crossing any surface, can be found by integrating the acoustic intensity over this surface

Pi = Z

Ai

IidA. (3.3)

According to that, the dissipated power in any one-dimensional acoustic element can be com- puted with following expression:

P_diss=P₁−P₂ = Z

A1

I₁dA− Z

A2

I₂dA (3.4)

whereI₁ and A₁ are the acoustic intensity and the cross-sectional area at the inlet and those values with index ”2” are referenced to the outlet. IfPdiss>0the acoustic power is dissipated,

1More general derivation is given in [1].

2More information about the meaning and the computation of cut-off frequency for annular geometry can be taken from [3] and for two-dimensional problems from [4].

P_diss<0 indicates acoustic power generation. Splitting of the specific acoustic exergy⁰ and the mass flow densitym⁰ into incident and scattered parts leads to

P₁ = Z

A1

<(⁺₁ +⁻₁)(m⁺₁ +m⁻₁)> dA

P₂ = Z

A2

<(⁺₂ +⁻₂)(m⁺₂ +m⁻₂)> dA.

(3.5)

The acoustic mass velocities can be expressed through the specific acoustic exergy by means of the correlation betweenp^±,ρ^± and u^± (eq. 2.6 to 2.10), derived in section 2.1:

m⁺₁ =ρ₁u⁺₁ +ρ⁺₁u₁ =ρ₁u⁺₁ +p⁺₁u₁ c²₁ = ρ₁

c1

p⁺₁ ρ1

+u⁺₁u

= ρ₁ c1

⁺₁ (3.6)

m⁻₁ =ρ₁u⁻₁ +ρ⁻₁u₁ =ρ₁u⁻₁ +p⁻₁u₁ c²₁ = ρ₁

c₁

−p⁻₁

ρ₁ −u⁻₁u

=−ρ₁

c₁⁻₁ (3.7) m⁺₂ =ρ₂u⁺₂ +ρ⁺₂u₂ =ρ₂u⁺₂ +p⁺₂u₂

c²₂ = ρ₁ c₁

−p⁺₂

ρ₂ −u⁺₂u

=−ρ₂

c₂⁺₂ (3.8) m⁻₂ =ρ₂u⁻₂ +ρ⁻₂u₂ =ρ₂u⁻₂ +p⁻₂u2

c²₂ = ρ1

c₁ p⁻₂

ρ₂ +u⁻₂u

= ρ2

c₂⁻₂. (3.9) With this knowledge P1 and P2 can be computed by

P₁ = Z

A1

<(⁺₁ +⁻₁)(⁺₁ −⁻₁)ρ₁

c₁ > dA= Z

A1

< ρ₁ c₁

⁺₁2

− ⁻₁2

> dA=

= Z

A1

ρ₁

c₁ < ⁺₁2

> dA

| {z }

P₁⁺

+ Z

A1

−ρ₁

c₁ < ⁻₁2

> dA

| {z }

P₁⁻

(3.10)

and, respectively,

P₂ = Z

A2

−ρ₂

c₂ < ⁺₂2

> dA

| {z }

P₂⁺

+ Z

A2

ρ₂

c₂ < ⁻₂2

> dA

| {z }

P₂⁻

. (3.11)

These expressions show that the sound power can be split up into the incident and scattered parts.

To get rid of the time average, Y. Aur´egan and R. Starobinski [1] used the following formula to calculate the incident and scattered sound power for plane waves:

P^±= 1 2Re



 Z

Ai

m^±(^±)^∗dA



, (3.12)

where^∗ indicates complex conjugate and Re() is the real part. For one-dimensional systems integration over cross section of any duct can be replaced by multiplication with its cross- sectional area A_i:

P₁⁺= 1 2Re



A₁ρ1

c₁ ⁺₁(⁺₁)^∗

| {z }

real value

dA



= A1ρ1

2c₁ ⁺₁(⁺₁)^∗ (3.13) P₁⁻= 1

2Re

−A₁ρ1

c₁⁻₁(⁻₁)^∗dA

=−A1ρ1

2c₁ ⁻₁(⁻₁)^∗ (3.14) P₂⁺= 1

2Re

−A₂ρ2

c₂⁺₂(⁺₂)^∗dA

=−A2ρ2

2c₂ ⁺₂(⁺₂)^∗ (3.15) P₂⁻= 1

2Re

A₂ρ2

c₂⁻₂(⁻₂)^∗dS

= A2ρ2

2c₂ ⁻₂(⁻₂)^∗. (3.16) Using a new flow parameter³ Y_i = ^ρ_2c^iAi

i and inserting eq. (3.13) to (3.16) into the equation of the dissipated acoustic power leads to

P_diss = (Y₁⁺₁(⁺₁)^∗ +Y₂⁺₂(⁺₂)^∗)−(Y₁⁻₁(⁻₁)^∗+Y₂⁻₂(⁻₂)^∗). (3.17) Introducing the incident and scattered acoustic power amplitudesa and b, where

a=



 a₁ a₂



=





Y₁^1/2+₁ Y₂^1/2+₂



 (3.18)

b=



 b₁ b₂



=



 Y₁^1/2−₁ Y₂^1/2−₂



, (3.19)

leads to the expression for the dissipated energy

P_diss =a^†a−b^†b. (3.20)

The scattered acoustic power amplitude b is computed from the incident acoustic power amplitudea using a scattering matrix S:

b =Sa and b^†=a^†S^† respectively. (3.21) Therefore the dissipated acoustic power is calculated by:

P_diss=a^†a−a^†S^†Sa. (3.22)

The derivation of the scattering matrix for the acoustic power amplitudes a and b from the scattering matrix for the Riemann invariants is given in appendix A.3.

3Notice that Y. Aur´egan and R. Starobinski [1] have different value forYi, because they use instead of the specific acoustic exergy its effective value.

Figure 3.1.: Illustration of the incident and scattered acoustic power amplitudesa andb of an acoustic element.

The incident and scattered amplitudes a and b are illustrated in fig 3.1. It should be noticed that the incident acoustic power amplitude also can be expressed directly by the Riemann invariants using the transformation matrix_aM_{f g} from appendix A.3:

a=_aM_{f g}



 f1

g₂



=





(c1+u1)Y₁^1/2 0 0 (c₂−u₂)Y₂^1/2







 f1

g₂



. (3.23)

3.2. The Scattering Matrix for the Acoustic Power Amplitudes and the Potentiality

Starting from eq. (3.22), the ratio of the dissipated sound power to the incident sound power is considered, wherePinc=a^†a =||a||² and an= _||a||^a :

P_diss

P_inc = 1−a^†S^†Sa

||a||² = 1−a^†_nS^†San. (3.24) Therefore Pdiss = Pinc −Pscat, the ratio of the scattered acoustic power to the incident acoustic power can be computed by

η_si = P_scat Pinc

=a^†_nS^†Sa_n. (3.25)

The matrix S^†S is a hermitian matrix because of the property (S^†S)^† =S^†S. Therefore it can be diagonalized

T^†(S^†S)T =



 λ₁ 0

0 λ2



, (3.26)

where λ_i are the real eigenvalues and the column vectors of the unitary matrix T are the eigenvectors of S^†S. It should be considered that S^†S is a positive semidefinite matrix⁴, therefore allλ_i are non-negative.

4 a^† S^†S

a=b^†b=||b||²>0 for anya.

Introducing the vector d =T^†a_n with d^†d= 1, eq. (3.24) is written as

η_si =a^†_nT



 λ1 0

0 λ₂



T^†a_n =d^†



 λ1 0

0 λ₂



d=

2

X

i=1

λ_i||d_i||², (3.27)

The values λ_max and λ_min represent the maximum and the minimum of the potentially generated acoustic power⁵. The eigenvaluesλi allow to describe the behaviour of an acoustic system. There are three special cases:

a) The acoustic power generation or dissipation is absent if both eigenvalues are equal 1, i.e. S is an unitary matrix (S^†S=I_identity).

b) The system is passive if the acoustic power of the incoming waves is always greater than the acoustic power of the outgoing wavesη_si <1, i.e. both eigenvalues are smaller than 1.

c) The system is active respectively (ηsi >1), if both eigenvalues are greater than 1.

Concluding,λ_maxandλ_minrepresent only limits of the generated or dissipated acoustic power.

The exact value ofη_si can be calculated by eq. (3.27), where eachλ_i is more or less weighted by the vector d. The vector d depends on the ratio between both of the incident acoustic amplitudes a₁ =Y₁^1/2+₁ and a₂ =Y₂^1/2+₂. That means the influencing parameters are the flow and geometry parameters (u,c,ρ,A), the amplitudes of the incident Riemann invariants f₁, g₂ and their phase difference.

3.3. The Ratio of the Scattered Acoustic Power to the Incident Acoustic Power

This section shows, how to compute the ratio of the scattered acoustic power to the incident acoustic power η_si = ^P_P^scat

inc from the scattering matrix for the acoustic power amplitudes and the incident Riemann invariants f₁ and g₂. Starting with eq. (3.24) and replacing S^†S with H leads to

η_si = a^†Ha

||a||². (3.28)

5Y. Aur´egan and R. Starobinski [1] consider the potentialities of the dissipated acoustic power, which are obtained from the potentialities of the generated acoustic power: ξ_max=λ_min,ξ_min=λ_max.

The incident amplitudeacan be expressed by the Riemann invariants according to eq. (3.23).

It leads to

η_si =



 f₁^∗ g₂^∗





t

aM_{f g}H _aM_{f g}



 f₁ g₂





(c1+u1)²Y1||f1||²+ (c2−u2)²Y2||g2||² =

=

c₁Y₁^1/2f₁^∗





(1 +M₁) (1−M₂)^Y

1/2 2 c2

Y₁^1/2c1

g^∗₂ f₁^∗





t

Hc₁Y₁^1/2f₁





(1 +M₁) (1−M₂)^Y

1/2 2 c2

Y₁^1/2c1

g2

f1



 c²₁Y1||f1||²

(1 +M1)² + (1−M2)^2c_c²²2^Y2 1Y1

||g2||²

||f₁||²

=

=





(1 +M₁) (1−M2)^Y

1/2 2 c2

Y₁^1/2c1

g^∗₂ f₁^∗





t

H





(1 +M₁) (1−M2)^Y

1/2 2 c2

Y₁^1/2c1

g2

f1



 (1 +M₁)²+ (1−M₂)²

Y₂^1/2c2

Y₁^1/2c1

2

_||g

2||

||f1||

2 . (3.29) f₁, g₂ and their complex conjugates can be written in polar form:

f₁ =||f₁||e^i(ωt+ϕf1) g₂ =||g₂||e^i(ωt+ϕg2)

)

=⇒ g₂ f1

= ||g₂||

||f1||e^{iϕg2−iϕf1} (3.30) f₁^∗ =||f₁||e^{−i(ωt+ϕf1)}

g₂^∗ =||g₂||e^{−i(ωt+ϕg2)} )

=⇒ g^∗₂

f₁^∗ = ||g₂||

||f₁||e^{iϕf1−iϕg2}. (3.31) Introducing a new flow parameter β = ^Y

1/2 2 c2

Y₁^1/2c1

= q_αρ

2c2

ρ1c1 and using the derived equations for the Riemann invariants leads to η_si, expressed by an amplitude ratio between the incident waves ^||g_||f^2||

1|| and by their phase difference ϕ=ϕ_g2 −ϕ_f1:

η_si =





(1 +M₁) (1−M₂)β^||g_||f^2||

1||e^−iϕ





t

H





(1 +M₁) (1−M₂)β^||g_||f^2||

1||e^iϕ



 (1 +M₁)²+ (1−M₂)²β²_||g

2||

||f1||

2 . (3.32)

Notice that η_si is not influenced by the frequency of propagating waves directly. If the scattering matrix S of an acoustic element is not function of frequency, the generated power in this element would be the same for all frequencies. An example of this kind of element is given in section 5.3.

There are also possibilities to expressη_si by other combinations of the Riemann invariants, e.g. by the Riemann invariants at the inlet f₁ and g₁. At first η_si has to be written as a function of f₁ and g₂ using matrix notation:

η_si =



 f₁^∗ g^∗₂





t

aM_{f g}^† H _aM_{f g}



 f₁ g₂







 f₁^∗ g₂^∗





t

aM_{f g a}^† M_{f g}



 f₁ g2





, where (3.33)

aM_{f g} =





(c₁+u₁)Y₁^1/2 0 0 (c₂−u₂)Y₂^1/2



=Y₁^1/2c₁





1 +M₁ 0 0 (1−M₂)β



. (3.34)

aM_{f g} transforms the incident acoustic power amplitude vector a to the incident Riemann invariants f₁ and g₂. Next, g₂ is expressed by f₁ and g₁ using the scattering matrix for the Riemann invariants S_{f g}:

g1 =Sf g,11f1+Sf g,12g2 =⇒ g2 = 1

S_{f g,12}g1−Sf g,11

S_{f g,12}f1. (3.35)

That leads to the correlation



 f1

g₂



=





1 0

−S_{f g,11} S_{f g,12}

1 S_{f g,12}





| {z }

incMupstr



 f1

g₁



. (3.36)

Finallyηsican be expressed by the amplitudes off1,g1and their phase differenceφ =φg1−φf1:

η_si =



 f₁^∗ g₁^∗





t

incM_{upstr a}^† M_{f g}^† H _aM_{f g inc}M_upstr



 f₁ g₁







 f₁^∗ g₁^∗





t

incM_{upstr a}^† M_{f g a}^† M_{f g inc}M_upstr



 f₁ g₁





=

=



 1

||g₁||

||f₁||e^−iφ





t

incM_{upstr a}^† M_{f g}^† H aMf g incMupstr



 1

||g₁||

||f₁||e^iφ







 1

||g1||

||f1||e^−iφ





t

incM_{upstr a}^† M^t_{f g a}M_{f g inc}M_upstr



 1

||g1||

||f1||e^iφ





. (3.37)

Notice that this function exists only if S_{f g,12}6= 0, i.e. g₁ and g₂ should be related. It always holds for subsonic problems.

4. Description of the MATLAB tool

In this section there is a description of the computation tool, which consists of the following MATLAB functions:

Name Description

potentiality computation.m computes the potentialities λ_min and λ_max and the scattering matrix for the acoustic power amplitudes.

potentiality plot.m plots the potentialitiesλ_minandλ_maxversus the sound wave frequency.

potentiality plot data cursor.m adapts the data cursor of the function po- tentiality plot.m.

find data index.m searches for the index of an structure array element with the specified frequency.

pgen plot f1g2.m computes and plotsη_si as a function of the incoming Riemann invariants f₁ and g₂. pgen plot f1g2 polar data cursor.m

pgen plot f1g2 cartesian data cursor.m

adapt the data cursor of the function pgen plot f1g2.m.

pgen plot f1g1.m computes and plotsη_si as a function of the Riemann invariants at the inletf₁ and g₁. pgen plot f1g1 polar data cursor.m

pgen plot f1g1 cartesian data cursor.m

adapt the data cursor of the function pgen plot f1g1.m.

pdiss fmincon f1g2.m searches for the maximum or the minimum of ηsi in a specified search domain.

Table 4.1.: Overview of the developed MATLAB tool.

An example of using all these functions is shown in appendix A.4. Most of these functions have a variable number of arguments, which are set in italics. The required variables are underlined. If remaining variables are not determined by the user, they are set to their default values, which are chosen adequate for the computation.

4.1. potentiality computation.m

[output struct]=potentiality computation(input struct,u1,c1,u2,c2,alpha,rho1,rho2) This function computes the scattering matrix S for the acoustic power amplitudes, the maximum and the minimum of the potentiality λ_max and λ_min. The computation algorithm is derived in section 3.2. The function requires a structure or a structure array input struct.

Each structure array element should have two fields with fixed names. The first field has the name freq, where the frequency of acoustic wave is stored. The second field keeps the input transfer or scattering matrix and has one of the following names:

Name Description

Transfer pu The transfer matrix for the pressure and velocity fluctuations.

Transfer fg The transfer matrix for the Riemann invariants.

Transfer fg The scattering matrix for the Riemann invariants.

Table 4.2.: The input transfer or scattering matrix

It is important that the field name is given correctly, because for each of these matrices there is a different computation path. The scattering matrix S is computed in steps

T_pu−→T_{f g} −→S_{f g} −→S, (4.1)

which are described in appendixes A.1 to A.3. Each transformation is implemented as an inline function within the main function. If the input matrix isT_pu, the function calculatesT_{f g}, than S_{f g} and finallyS. If, for example,S_{f g} is given, other matrices,T_puandT_{f g}, are not needed, so these computation steps can be skipped. Each element stores the matrixS (under the name

’Scatter pow’) and all of the matrices, computed on the way to it. Besides that, the function adds fields with the minimum and the maximum of the potentially generated acoustic power λmin and λmax to each structure array element. These fields have names lambda min and lambda max. The function returns this new structure array. The illustrations of input struct and output struct can be seen on Figures 4.1 and 4.2.

Other input parameters are velocity u, speed of sound c, cross-sectional area ratio alpha and density rho. Indices ”1” and ”2” denote the inlet and the outlet. If velocities are not specified, they are set to zero and other quantities have the same value at the inlet and at the outlet, because only their ratios are relevant, rather than the values themselves. There is one exception. If u1 is specified, the sound of speed c1 has to be put in for the computation of the Mach number M₁. If c1 is not determined, an error message appears. In other cases all not specified arguments are set automatically to default values.

Figure 4.1.: An example of a structure arrayinput struct.

Figure 4.2.: An example of a structure array output struct.

4.2. potentiality plot.m

[freq,lambda min,lambda max]=potentiality plot(input struct,name)

This function plots the maximum and the minimum of potentially generated power λmax

and λ_min. The function extracts the values of the fields freq, lambda max and lambda min from each element ofinput struct, which is returned bypotentiality computation.m. Than it creates three subplots: λ_max versus the corresponding frequencyf (blue dash-dot line),λ_min versus frequency (green dashed line) and both λ_i together versus f. The first two subplots are made for the case, if λ_max and λ_min have different orders of magnitude, therefore one of them can not be displayed well on the common subplot. Also on each subplot there is an additional red line, which indicates P_scat = P_inc the neutral state of an acoustic element. It serves as a boundary between the dissipated (above the red line) and generated (below the red line) acoustic power. The function use the m-file potentiality plot data cursor.m, which adapts a data cursor to display η_si and the corresponding frequency f. An example of the third subplot is presented in fig. 4.3.

Figure 4.3.: An example of a potentiality plot created by potentiality plot.m. λmax (blue dash-dot line) and λ_min (green dashed line) are plotted versus the corresponding frequency. The red line is a boundary between dissipated (above the line) and generated (below the line) acoustic power over the incident power.

.

The plots are saved automatically in the current folder as MATLAB figure, PDF and EPS.

The name of the file and the title of the plots is specified by the input string name. The default name is ”The minimum and maximum of the potentially generated acoustic power”.

For further use, the function returns three column vectors freq, lambda min and lambda max with the plot data.

4.3. find data index.m

[output]=find data index(input struct,freq min,freq max)

This function is an auxiliary tool, which searches in the structure array input struct for el- ements with frequency in the interval [freq min, freq max]. If only one boundary freq min is given, than the function searches for elements with the frequency within the interval [freq min−0.001Hz, freq min+0.001Hz] frequency with tolerance ±0.001 Hz. output is a matrix with two rows: frequency and index of the corresponding element. This tool allows to find quickly an element, if its frequency is known. It is important that the input structure array has a form, defined in section 4.1, because the tool compares values of the fields with the name freq.

4.4. pgen plot f1g2.m

[Pgen min,Pgen max,PHI deg,R,Pgen]=pgen plot f1g2(Scatter pow,u1,c1,u2,c2,alpha,...

rho1,rho2,cmap switcher,r min,r max,r step,angle number,name)

This function calculates and plots the ratio of the scattered acoustic power to the incident acoustic powerηsias a function of the amplitude ratio between the incident Riemann invariants

||g₂||

||f₁|| and their phase difference ϕ(see eq. (3.32)). The tool creates two figures (see fig. 4.4 and 4.5). The first one shows a polar plot with a radiusr= ^||g_||f^2||

1|| and an angleϕ. The value of η_si on each point is represented by a color. The same variables are plotted on the second created figure, but in Cartesian coordinate system.

The function requires a scattering matrix S for the acoustic power amlitudes as input, which is computed by the functionpotentiality computation.m. The next seven parameters are the same as the input parameters of the function potentiality plot.m. The remaining arguments are plot parameters: r min, r max and r step specify a range and a resolution of the amplitude ratio ^||g_||f^2||

1||. angle number determines an angular resolution. Both of the figures have 64 colors to displayη_si¹. The input argumentcmap switcher chooses the distribution of colors according to the ratio η_si. 0 means the default distribution (see figure 4.5), where the

1If necessary, the color number can be set manually by changing the value of the variablecolor step number in the function file.

Figure 4.4.: An example of a polar plot created by pgen plot f1g2.m. The color represents the ratio of the scattered acoustic power to the incident acoustic power. The amplitude ratio between the incident Riemann invariants ^||g_||f^2||

1|| and their phase differenceϕare equal to the polar coordinates. The colormap is created with the custom color distribution function.

Figure 4.5.: An example of a Cartesian plot created by pgen plot f1g2.m. The color rep- resents the ratio of the scattered acoustic power to the incident acoustic power.

The colormap is created with the default color distribution function.

.

smallest values are blue colored and the highest values are red colored. 1stays for the custom distibution (see figure 4.4), where η_si = 2 is red, η_si = 1 is green and η_si = 0 is blue. All the values bigger than 2 remain red. It means blue colored areas show the domains, where an acoustic system is passive. The behaviour of the system in red colored areas is active.

Nevertheless this distribution is only useful if the acoustic power generation spectrum is large enough. The figures are saved with a default name as MATLAB figure, PDF and EPS, if the stringname is not specified. This name also serves as the title of the figures.

To supply both of the figures with an appropriate data cursor, the function needs the functions pgen plot f1g2 polar data cursor.m and pgen plot f1g2 cartesian data cursor.m.

Besides the figures, this function can return matrices Pgen min and Pgen max, which contain the coordinates of points with minimal and maximal values of computed η_si. Notice that these points don’t represent the maximal or minimal values ofη_si, which can be achieved on a determined area. They are simply the maximal or minimal values achieved on the com- puted grid. Therefore this feature would never return the real maximum or minimum, until any of computed grid points randomly hits an extremum. Other output arguments are matrices PHI deg, R and Pgen with the plot data.

4.5. pgen fmincon f1g2.m

[output]=pgen fmincon f1g2(Scatter pow,phi deg,r,switcher,r min,r max,u1, c1, u2, c2, alpha, rho1, rho2)

This function searches for the minimum or maximum ofη_si in the amplitude ratio interval [r min, r max]. The core of this tool is the inline function, which calculates η_si from the scattering matrix S for the acoustic power amplitudes, ^||g_||f^2||

1|| and ϕ (see eq. (3.32)). For the computation it also requires additional parametersu1, c1, u2, c2, alpha, rho1 and rho2, which are set to their default values or can be determined by user. For the searching for the minimum or the maximum the Optimization Toolbox is needed, because of the using of the function fmincon. Depending of a value ofswitcher (”0” for minimum and ”1” for maximum)fmincon searches for the maximum/minimum of η_si. To find the extremum correctly, the initial point with coordinates r (stands for ^||g_||f^2||

1||) and phi (stands for ϕin degree) is required. The solver is very sensitive to the choice of the initial point. To be sure that the solver found the correct value, the initial point should be put near the expected maximum or minimum location. To do that, the user should investigate a plot of the dissipated or generated power, created by pgen plot f1g2.m. Changing of any optimization options has to be done manually in .m- file. The default values are generated with MATLAB function optimset(@fmincon). The following options are set manually:

Property Value TolX ’1e-006’

LargeScale ’off’

Algorithm ’activate-set’

Table 4.3.: Optimization options of fmincon.

Notice that fmincon is able to find only one value for an extremum. If there are more than one, it usually finds the nearest one to the initial point. The function returns a column vector output withϕ (in degree), ^||g_||f^2||

1|| and ηsi as columns.

It is important to check the reliability of obtained results:

a) The maximum or minimum should be within the range between λ_min and λ_max.

b) The location of the maximum or minimum should be close to the location, expected after investigating an appropriate scattered power plot.

4.6. pgen plot f1g1.m

[Pgen min,Pgen max,PHI deg,R,Pgen]=pgen plot f1g1(Scatter pow,Scatter fg,...

u1,c1,u2,c2,alpha,rho1,rho2,cmap switcher,r min,r max,r step,angle number,name)

This function calculates and plots the scattered acoustic power over the incident acoustic power η_si as a function of the amplitude ratio between the Riemann invariants at the inlet

||g₁||

||f1|| and their phase difference φ (see eq. (3.37)). The function creates polar and Cartesian plots and has the same input arguments as pgen plot f1g1.m and one additional argument Scatter fg, because this tool needs S and S_{f g} as input matrices to computeη_si as a function of f₁ and g₁.

5. Illustrative Examples

5.1. Uniform Tube with Non-Viscous Adiabatic Mean Flow

The transfer matrix of non-viscous adiabatic flow is given by U. Neunert [5]:



 f₂ g₂



=





exp (−ik⁺l) 0 0 exp (ik⁻l)







 f₁ g₁



, (5.1)

where the wave number for non-viscous flow is given by:

k^±= 2πf

c±u. (5.2)

At first, the analytical solution for the value of the dissipated acoustic power is derived. It is expected that such a system doesn’t dissipate or generate acoustic power at all. To compute Pdiss eq. (3.22) is used:

P_diss = ¯a^ta−¯b^tb=||a||²− ||b||²,

where a and b are functions of the Riemann invariants according to eq. A.6 and A.7. That leads to

Pdiss= (c1 +u1)²Y1f¯1f1+ (c2−u2)²Y2g¯2g2−(c1−u1)²Y1g¯1g1−(c2+u2)²Y2f¯2f2, (5.3) where

f¯₁f₁ =||f₁||²

¯

g₁g₁ =||g₁||² f¯₂f₂ =f₁exp (−ik⁺l)f₁exp −ik⁺l

= ¯f₁f₁exp ik⁺l

exp −ik⁺l

=||f₁||²

¯

g₂g₂ =g₁exp (ik⁻l)g₁exp ik⁻l

=||g₁||².

There is no area change, the flow is adiabatic and non-viscous: A₁ =A₂ =A, c₁ =c₂ =c, u₁ =u₂ =u, ρ₁ =ρ₂ =ρ and, consequently, Y₁ =Y₂ =Y.

That leads to

P_diss =Y (c+u)²||f₁||²+ (c−u)²||g₁||²−(c−u)²||g₁||²−(c+u)²||f₁||²

= 0.