A semi-analytical solution for acoustic wave propagation in varying area ductswith mean flow
S R YEDDULA(1) ; A S MORGANS(2)(1) Imperial College London, United Kingdom, s.yeddula18@imperial.ac.uk(2) Imperial College London, United Kingdom, a.morgans@imperial.ac.uk

Abstract
A semi-analytical solution for the propagation of plane acoustic waves in a varying area duct, sustaining a meanisentropic flow with axial gas property gradients, is developed. A second order differential equation (ODE)
for acoustic pressure is derived from the linearised Euler equations in the frequency domain, by neglecting
the communication between acoustic and entropy disturbances. This ODE has axially varying coefficients andis solved using an adaptive WKB approximation method, so that the obtained wave-like solution consists ofseparate amplitude and phase factors, expressed as the superposition of downstream and upstream propagating
wave components. This solution is valid at sufficiently large frequencies, meaning the frequencies consideredmust be low enough to ensure a predominantly one-dimensional acoustic field with negligible diffusive effects,
but large enough for WKB method validity. A test case representing a typical industrial gas turbine combustor
outlet converging section with a high mean-flow acceleration is considered for validation. Analytical results
are compared to numerical solutions of linearised Euler equations neglecting acoustic-entropy coupling (2LEE).
When the limiting frequency criteria are satisfied, the semi-analytical solution agrees well with the 2LEE results
at even flow Mach-numbers around 0.5.
Keywords: Semi-analytical, 1-D acoustic field, varying-area duct, WKB method

1

INTRODUCTION

Thermoacoustic instabilities, also known as combustion instabilities, occur due to a two-way coupling between
unsteady heat release and acoustics. If acoustic pressure fluctuations are in-phase with unsteady heat releaseat the flame, and if the resulting energy gain of perturbations exceeds the acoustic energy efflux[1], the amplitude of acoustic oscillations grows. This causes high levels of noise and system vibration until the effectsof acoustic non-linearities leads to limit-cycling or other non-linear oscillations. Thermoacoustic instabilities
are highly undesirable. They are also difficult to simulate, the disparity between the acoustic length scales,
which are O(metre), and the flame wrinkling, O(millimetre), and flame reactions O(micrometre), makes it very
computationally costly.
Computational methods which couple separate treatments for the acoustic waves and the flame have been shownto have much promise [2]. The flame is normally treated using large eddy simulations, while the acoustic waves
can be treated using low order network modelling[2, 3] or Helmholtz solver[4]. This work focuses on low order
network modelling, which employs acoustic waves analytical solutions that have the ability to account for mean
flow effects. However, acoustic wave models for ducts with mean flow, temperature and area gradients, are currently limited to specific simple cases. It is desirable for network models to have an analytical solution, which
can be represented as a superposition of upstream and downstream propagating acoustic wave components, thisoffering more physical insight and being computationally effective. Analytical solutions offer the benefits ofguaranteed convergence, irrespective of grid size, make parametric studies much easier and are therefore helpful
in identifying efficient ways of altering the system to attenuate thermoacoustic oscillations during the design
stage.
The WKB method, named after Wentzel, Kramers, and Brillouin[5], is a technique traditionally used in wavemechanics to obtain approximate wave-like solutions for second order differential equations with variable coeffi-

5633

Actuator
Q̇

T (x)

Anechoic
boundary

M(x)

50mm

250mm

A(x)
x0

x1

x2

x

500mm

Figure 1. (Left) Sketch of varying area duct (x1 to x2 ). A(x), M (x), and T (x) represent the duct cross-sectionalarea, flow Mach number and mean-temperature for given x. For the validation test cases, actuation upstream of
the area variation and an anechoic boundary upstream of actuation are considered. (Right) Sketch of conicalduct test-case with inlet & outlet diameters of 250mm, 50mm, and length of 500mm.

cients. Cummings[6] used a slightly modified WKB method with separate amplitude and phase factors, insteadof asymptotic expansions, to determine the analytical solution for an acoustic field in a uniform area duct withtemperature gradient, but without mean-flow, for sufficiently large frequencies. This approach was later extended
to non-isentropic flows in uniform area duct with temperature gradient and practical mean flow velocities by Li
and Morgans[7]. Although Dokumaci[8] presented a semi-analytical solution for a varying area duct with nonzero mean flow, using a numerical matrizant method, the limiting frequency criteria were developed for a zeromean flow case only, leading to lack of clarity at higher Mach numbers. Furthermore, only the lower limiting
frequency criterion was considered; the upper limit beyond which acoustic field becomes two dimensional was
not looked into.
The present study extends the one-dimensional semi-analytical acoustic solution presented by Li and Morgans[7]for a uniform area duct using the WKB method, to isentropic flows in a duct of varying area as shown in
Fig.1(left). A two step procedure is adopted. Firstly a new differential equation only in terms of acoustic pressure is derived for such a duct for the first time. This is then solved using an adaptive WKB method. Terms
with flow Mach-number of order O(M 3 ) and higher are neglected, a reasonable assumption because combustors
and after-burners sustain low Mach number flows, typically of the order O(0.01) to O(0.2), to maintain flame
attachment.
Further sections are organized as follows, the derivation of the acoustic wave equation and its solution usinga modified WKB approximation technique for isentropic flows is presented in Section 2. Section 3 discusses
the frequency limits which apply to the approximate solution. Section 4 defines the terminology adopted fortransfer functions, and explains the test domain considered in the validation study. Results for the proposed
method against numerically solved linearised Euler equations are presented in Section 5. Conclusions are drawn
in the final Section.

2

ISENTROPIC VARYING AREA FLOWS

2.1 Acoustic wave equation
A flow that is not affected by fluid viscosity, thermal conductivity or heat transfer can be described as an
isentropic flow. Based on principle of linearisation, all flow variables can be separated into a mean time¯ and a perturbation component, denoted by ()0 . The conservation equations
averaged component, denoted by (),for mass, momentum, and energy in an isentropic flow, linearised in the time domain, were given by Marbel andCandel[9]. These are converted to the frequency domain by assuming a time-dependance of the form y0 = ŷeiωt

5634

to the fluctuating quantities, where ω represents the complex angular frequency, and i =√−1, as follows,



d p̂
d û
1 dp
2 du
iω + γMp̂ + u+û + γ p=0dx
dx M 2 dx
dx(1)

negligiblez }| {


d û
du p̂
du
du
1 d p̂iω +û + u + u+= u σ̂ ≈ 0dx
dx
dx γ p ρ dx
dxiω σ̂ + u

d σ̂
=0
dx

(2)

p = ρRg T ;and(3, 4)

where ρ, p, T , u, σ̂ , Rg , denote mean-density, mean-pressure, mean-temperature, mean-axial velocity, entropy
disturbance normalized with specific heat at constant pressure, and universal gas-constant respectively. Rg andγ (ratio of the specific heats), are assumed to be constants throughout the domain, with values of 287J/kgK and
1.4 respectively.
Normalisation, shown in Eqs.(5, 6), is applied to the area and mean-temperature gradients. Notably, these gradients cannot be arbitrary as the net heat-efflux from the system should be zero for an isentropic flow, as given
by Eq.(7). This leads to a relation between the gradients given by, Eq.(8).
α=

1 dA
;
A dx



du
dT
Q̇ = ρu c p
+u
=0
dx
dx

β=

1 dT
,
T dx

⇒

β=

(5, 6)
αM 2 (γ − 1)
1 − M2(7, 8)

On time-averaging conservation equations, and using Eq.(8), other mean flow gradients can be expressed as,
1 du
−α
=
;
u dx (1 − M 2 )

1 dp
αγM 2
=
;
p dx
(1 − M 2 )

1 dρ
αM 2
=
;
ρ dx
(1 − M 2 )1 dc βαM 2 (γ − 1)
= =
;
c dx
2
2 (1 − M 2 )(9, 10, 11, 12)

As the energy equation, Eq.(3), doesn’t have a source-term, the entropy disturbance at the inlet is simply
convected through the nozzle, its magnitude unchanged. The communication between the acoustic field andentropy, as implied by term on right hand side of momentum equation Eq.(2), is of order O(M 2 )[10]. Acousticentropy coupling is then neglected in this analysis, since the flow Mach number is low downstream of the flame.
By dividing Eq.(1) by γ p, and Eq.(2) by u, coefficients of d û/dx are normalized in both equations, and theequations become synonymous with Eq.(13 & 14). Now these equations are subtracted so that û is explicitly
written in terms of p̂ and d p̂/dx, as shown in Eq.(15).
A p̂ + Bd p̂
d û
+Cû += 0;dx
dx

D p̂ + Ed p̂
d û
+ F û +=0dx
dx(13, 14)

where the coefficients A, B,C, D, E, F are functions of space coordinate, x, and wavenumber, ko . Also,û = G p̂ + H

d p̂
,
dx


where G =

D−A
C−F




and H =

E −B
C−F


(15)

Eq.(15) and it’s differential with x can replace û and d û/dx of Eq.(13), in terms of p̂. On neglecting terms oforder higher than M 2 and using relations Eq.(9) to Eq.(12), an O.D.E only in terms of p̂ is obtained,






2αM d 2 p̂
2M dα d p̂
2
2
2 2
2 dα
1 − M2 +
−
2iMk
−
α
+
2αM
−
+
k
+
3α
M
−
3M
p̂ = 0
(16)
o
o
iko
dx2
iko dx dx
dx5635

To summarize, Eq.(16) is the differential equation for the acoustic pressure disturbance, applying to isentropic
flow in a varying-area duct. The assumptions upon which it is based are: (a) |ko |2 >> |αM|2 , (b) neglecting
terms of O(M 3 ) and higher, (c) constant gas properties γ, Rg throughout the flow domain, (d) weak acousticentropy coupling. The above equation can readily be verified for simpler cases of no mean flow, or uniform
area duct. When both of these apply, the classic one-dimensional Helmholtz equation is recovered.
2.2 Semi-analytical solution using adaptive WKB method
The adaptive WKB technique used by Li and Morgans[7] is extended to a varying-area duct with isentropic
mean flow by expressing pressure perturbation as,
Z x
p̂ = Iexp(a + ib) dx(17)x1

where, I is an arbitrary constant, a and b are two real x-dependant variables that are determined using O.D.E,
Eq.(16). The obtained acoustic field solution can be represented in terms of upstream and downstream propagating waves at high frequencies, i.e. where the length scale over which the area change occurs is larger thanthe length scale of acoustics. As the imaginary part of wave-number, ko , is typically much smaller than the realpart in thermo-acoustic systems[11], ko is considered to be real. Upon substituting Eq.(17) into Eq. (16) andequating real and imaginary parts to zero, a solution for a and b is obtained. Firstly upon equating the real partto zero, we arrive at Eq.(18),
 




1 dα
da
2aα
2
2
2
2
+ 2
+ a +
+ aα = O M 2(18)
ko − b 1 − M +2Mko b + 2Mko b
2
dx
ko
ko dx
When the frequency is sufficiently large, a crucial requirement for validity of the WKB method, b, which is oforder ko , is therefore much larger than a, terms with ko in the denominator, and terms of order M 2 . Hence terms
1, 2, and those of O(M 2 ) in Eq.(19) are neglected. These assumptions are further dealt with while discussing
low frequency limits in Section 3.
Term 1

Term 2

z
}|
{ z
{

 }|

da
2aα
1 dα
2
2
2
+ aα + 2Mko b+ O M2∴ b ≈ ko >> a +
+ 2
2
dx
ko
ko dx

(19)

For acoustic systems that satisfy Eq.(19), Eq.(18) would take the form of a simple quadratic equation in b.
ko 2 − b2 (1 − M 2 ) + 2Mko b ≈ 0(20)

Now with possible solutions as, b+ = −ko /(1 + M), and b− = ko /(1 − M), respective a± values can be determined by equating the imaginary part to zero,
Negligible

2ab 1 − M

2



z


}|
{
 db 2αM 2
2aM dα 2αM da
2− 2aMko + 1 − M+
b + αb 1 − 2M −−= O M2
dx
ko
ko dx
ko dx
2

(21)

Terms with ko in the denominator of Eq.(21) and those of the order M 2 and above are small in magnitude ascompared to other terms, hence neglected. Thus we obtain,ααM 2 (γ − 3)
± αM +(22)
2
4
Results for a and b are substituted into Eq.(17) and using superposition, the overall analytical solution can be
expressed as sum of the two simplified solutions,
a± = −p̂ (x, ω) = C1 + P1 + (x, ω) +C1 − P1 − (x, ω)

5636

(23)

with
P1 ± =



A1
A



 

1/2 exp M 2 γ−3 ± M


Z x
1
1
8
dx


 

 exp ∓iω
x1 c ± u
exp M 2 γ−3
±M
8

(24)

where, C1 + and C1 − are constants that can be determined using boundary conditions, M1 represents the Machnumber of the flow and A1 is area of the nozzle, both at inlet (x = x1 ).
The velocity perturbation, û, is determined from Eq.(15), and can be written as,û (x, ω) =


1
+ +
− − −
B+
1 C1 P1 − B1 C1 P1
ρc

(25)

where,
B±
1 =

3



 α
1
αM 2 (γ − 3)
2
±
ik
1
−
M
−
αM
∓
o
iko (1 − M 2 ) − 2αM
2
4

(26)

VALIDATION TEST-CASE

The accuracy of the semi-analytical solution is evaluated by applying it to the test-case shown in Fig.1(Right),
which consists of a conical converging-area duct. For the WKB method, the frequency must be “high”, yet the
acoustic field must remain 1-D, meaning that the frequency must not reach such high values that higher orderacoustic modes, for example circumferential or radial modes, become cut-on. Thus the frequency must lie in anintermediate range, high enough for WKB-method validity yet low enough for the acoustic field to remain 1-D.
3.1 Lower frequency limit
For a uniform mean-flow, in the absence of temperature and area gradients, the strengths of the upstream
and downstream propagating waves are uncoupled. Coupling of these wave strengths mainly arises due to thetemperature and area gradients along the flow[8]. The flow Mach number distribution that results from these
gradients, as a function of x, has only a secondary contribution to the coupling of wave components.
The lower frequency limit, or decoupling criterion, for acoustics in a varying-area duct with temperature gradient, is framed from the terms that are neglected in Eq.(19), while approximating the wave strengths as independent. Using a± given by Eq.(22) and neglecting terms of order O(M 2 ) and higher we get,s⇒ b ≈ ko >>dα−α 2 1 dα−± α 2M ∓ M+ O (M 2 )
4
2 dx
dx

(27)

For a conical duct geometry, the radius changes linearly with x, thus dα/dx = −α 2 /2. Hence Eq.(27) becomes,r3α 2 Mko >>+ O (M 2 )(28)
2
Numericallyp
this criterion is not as severe as for other geometries, for example in an exponential area changecase ko >> |α + O(M)| would be the limit. For the current test-case as shown in Fig.1(Right), the frequency
limit determined from Eq.(28), at the outlet for a flow Mach-number of 0.5 is approximately fl∗ >> 465Hz .
3.2 Higher frequency limit
The frequency at which the first non-axial modes occur in a conical section sets the high frequency limit, andis given by,
0.293cx
fh∗ =
(29)
rx

5637

where, cx and rx are local speed of sound and radius of the duct, respectively. Considering the test case inlet
diameter, shown in Fig.1(Right), the first higher order mode occurs at fh∗ ≈ 1880Hz and so the acoustic field
frequency should be less than fh∗ .
3.3 Setting initial conditions and defining transfer functions
At the duct inlet, x1 , of Fig.1(Right), we enforce an acoustic impedance condition by specifying the relation
between the pressure perturbation, p̂1 , and velocity perturbation, û1 , as given by, p̂1 = ρ 1 c1 Z1 û1 , where Z1 represents the non-dimensional acoustic impedance at the inlet. Dowling’s method[12] is adopted in which acoustic
waves are considered to have not reached the boundary at the outlet, x2 , (semi-infinite acoustic boundary). Thus
for a prescribed/enforced frequency, the eigenvalue problem is converted into an initial value problem. Thisfacilitates numerical solution, allowing higher-order schemes, like 4th order Runge-Kutta, to be readily implemented, and also facilitates comparison of the full frequency response with models. Transfer functions which
capture the axial evolution of the field are used, the same as those considered in [7], given by,
Fp (x, f ) =4p̂ (x, f );ρ 1 c1 û1 ( f )Fu (x, f ) =û (x, f )û1 ( f )(30, 31)

RESULTS

The proposed analytical solution is used to predict the acoustic field for the test geometry of Fig.1(Right), across
two different flow Mach numbers and two different frequencies. The magnitude of the normalised area gradient,α, given by Eq.(5), changes from 2.5 at the inlet to 6 at outlet for this geometry. Frequencies of 1600Hz and
1040Hz are chosen as they satisfy both Eq.(28) and Eq.(29) for the range of flow Mach numbers considered.
In Fig.2, results are presented for an inlet flow Mach number M1 = 0.005, corresponding to an outlet Mach
number of 0.03. As discussed in Section 3, when the flow Mach number is very low, the decoupling criterion
limit given by Eq.(28) approaches zero, making the final result exact.
3

3

2

2

1
1

1
1

0

0

-1

-1
0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

Figure 2. Evolution of F p and Fu with x. T 1 = 1600 K, M1 = 0.005, M2 = 0.03, f = 1600Hz, and Z1 = -1.
(Top left): transfer function gain |F p |, (Bottom left): transfer function phase ∠F p ,(Top right): transfer function gain |Fu |, (Bottom right): transfer function phase ∠Fu .
When the inlet flow Mach number is increased to 0.07, the corresponding outlet Mach number is 0.51. The

5638

results for this case, shown in Fig.3, confirm that that the analytical solution accurately matches the 2LEE
solution, with very small deviations only towards the exit, due to higher flow Mach numbers and α values inthis regime. Predictions for the same Mach number profile, but at a reduced frequency of 1040Hz is presented
in Fig.4, with the solution remaining accurate.
4

4

3

3

2

2

1
1

1
1

0

0

-1

-1
0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

Figure 3. Evolution of F p and Fu with x. T 1 = 1600 K, M1 = 0.07, M2 = 0.51, f = 1600Hz, and Z1 = -1.
(Top left): transfer function gain |F p |, (Bottom left): transfer function phase ∠F p ,(Top right): transfer function gain |Fu |, (Bottom right): transfer function phase ∠Fu .
5

5.5

4
4
3
2.5

2
1
1

1
1

0

0

-1

-1
0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

Figure 4. Evolution of F p and Fu with x. T 1 = 1600 K, M1 = 0.07, M2 = 0.51, f = 1040Hz, and Z1 = -1.
(Top left): transfer function gain |F p |, (Bottom left): transfer function phase ∠F p ,(Top right): transfer function gain |Fu |, (Bottom right): transfer function phase ∠Fu .

5639

5

CONCLUSIONSThis paper presents an analytical solution for one-dimensional acoustic field in a varying area duct sustaining an
isentropic mean flow along with temperature gradients. An acoustic wave equation is derived from the linearised
conservation equations by neglecting acoustic-entropy coupling, and higher order flow Mach number terms. Thissecond order ODE in terms of acoustic pressure is solved using an adaptive WKB method, which assumes awave like solution with separate amplitude and phase factors. The final solution thus obtained is a superpositionof upstream and downstream propagating waves and is computationally easy to evaluate. The frequency regime
within which the proposed solution works accurately is also derived, needing to be high enough for WKBmethod validity but low enough for a 1-D acoustic field. By applying the proposed analytical solution to aconical converging duct, it was shown that even for flow Mach numbers of around 0.5, the acoustic field is
captured accurately and is comparable to the 2LEE numerical solution when the frequency criteria are satisfied.

ACKNOWLEDGEMENTS
The authors would like to acknowledge the European Research Council (ERC) Consolidator Grant AFIRMATIVE (2018-2023) for supporting the current research.

REFERENCES[1] Chu, B.T., On the energy transfer to small disturbances in fluid flow (Part I). Acta Mechanica, 1(3), 1965,
pp.215-234.
[2] Han, X.; Li, J.; Morgans, A. S. Prediction of combustion instability limit cycle oscillations by combining flame describing function simulations with a thermoacoustic network model, Combustion and Flame,162(10), 2015, 3632-3647.
[3] Li, J.; Yang D.; Luzzato, C.; Morgans, A. S. Open Source Combustion Instability Low Order Simulator(OSCILOS–Long) Technical report, 2015.
[4] Laera, D.; Schuller, T.; Prieur, K.; Durox, D.; Camporeale, S. M.; Candel, S., Flame Describing Function analysis of spinning and standing modes in an annular combustor and comparison with experiments,
Combustion and Flame, 184, 2017, 136-152.
[5] A. Schlissel, The initial development of the WKB solutions of linear second order ordinary differential
equations and their use in the connection problem, Hist. Math. 4, 1977, 183–204.
[6] Cummings, A. Ducts with axial temperature gradients: an approximate solution for sound transmission andgeneration. Journal of Sound and Vibration, 51(1), 1977, 55-67.
[7] Li, J.; Morgans, A. S. The one-dimensional acoustic field in a duct with arbitrary mean axial temperature
gradient and mean flow, Journal of Sound and Vibration, 400, 2017, 248-269.
[8] Dokumaci, E. An approximate analytical solution for plane sound wave transmission in inhomogeneousducts, Journal of sound and vibration, 217(5), 1998, 853-867.
[9] Marble, F. E.; Candel, S. M. Acoustic disturbance from gas non-uniformities convected through a nozzle,
Journal of Sound and Vibration, 55(2), 1977, 225-243.
[10] Li, J.; Yang, D.; Morgans, A. S. The effect of an axial mean temperature gradient on communication
between one-dimensional acoustic and entropy waves, International Journal of Spray and Combustion Dynamics, 10(2), 2018, 131-153.
[11] T. Lieuwen, Unsteady Combustor Physics, Cambridge University Press, 2012.
[12] Dowling, A. P.; Mahmoudi, Y. Combustion noise. Proceedings of the Combustion Institute, 35(1), 2015,65-100.

5640