Adjoint-based computation of shape sensitivity in a Rijke-tube
Georg A. MENSAH(1) , Alessandro ORCHINI(2) , Jonas P. MOECK(3)(1) CAPS Lab Department for Mechanical and Process Engineering, ETH Zürich, Switzerland, gmensah@ethz.ch(2) Institut für Strömungsmechanik und Technische Akustik, TU Berlin, Germany(3) Department of Energy and Process Engineering, NTNU Trondheim, Norway,

Abstract
Geometrical shape is the most important parameter defining acoustic spectra. While adjoint-based sensitivity
analysis of shape alterations is a well established technique in other engineering disciplines such as airfoil
design, it is less common in acoustic and especially thermoacoustic contexts. Only recently the thermoacoustic
community has begun to deploy such methods for stability assessment. The present paper aims at continuing
this effort by discussing the challenges of domain shape deformations to adjoint-based gradient calculations
and a parameter-free implementation into a state-of-the-art thermoacoustic Helmholtz solver. Deformations of asimple thermoacoustic system – a Rijke tube – will be studied as an example case to validate the method. Despite
this practical focus on thermoacoustics, the results of the paper are fruitful for other acoustic topics involving
stability analyses, as the underlying concepts readily apply to any eigenvalue problem linear or nonlinear in its
eigenvalue. Indeed, the gradient information calculated from first-order theory is the first building block for
efficient algorithms to optimize the shape of acoustic devices.
Keywords: adjoints, thermoacoustics, shape derivatives, eigenvalue problems

1

INTRODUCTION

The most influential parameter to acoustic resonance is shape and this paper discusses how to efficiently compute the influence of shape deformations to thermoacoustic stability based on an adjoint method.
In gas turbine engineering a critical aspect is thermoacoustic stability of the combustor. This is the robustness ofthe combustor against self-excited pressure fluctuations that arise from a constructive feedback of the fluctuationsof the rate of heat release in the flame and the acoustics of the combustion chamber. A linear modeling strategy
for this phenomenon is the use of a Helmholtz equation that is equipped with a term modeling the acoustic
feedback due to heat release [1]. It reads:∇ · c2 ∇ pb + ω 2 pb + F(ω)∇ pb|~xref ·~nref = 0
in Ω.
(1)
Here, pb is the Fourier transform of the pressure fluctuation, c the local speed of sound, ω the complex eigenfrequency, F(ω) is a function modeling the flame response linking the velocity fluctuation in the direction ~nrefat the point ~xref to the fluctuation in heat release rate in the flame and accounting for the thermodynamic properties of the flame. Ω denotes the considered domain. It is confined by its boundary ∂ Ω at which impedance
conditions are set. Mathematically these are mixed-type boundary conditions
D(ω) pb + ∇ pb ·~n = 0

on ∂ Ω.

(2)n denotes the outward-pointing
where D is typically defined from a normalized impedance Z as D = iωc
Z and ~
unit normal. The equations (1) and (2) present an eigenvalue problem. In contrast to an acoustic Helmholtz
equation, the thermoacoustic Helmholtz equation is highly non-linear in its eigenvalue, and non-self-adjoint due
to the flame response term and possibly due to the boundary conditions if they are not energy conserving). The
review [2] contains a recent overview on nonlinear eigenvalue problems and their mathematics.
How the eigenvalue ω depends on the domain Ω is the key question discussed in this paper. In particular, how
the eigenvalue changes with infinitesimal changes in the domain shape at first order is discussed. Such a change

5489

Ω(0)

distortΩ(ε)

Figure 1. Illustration of a distortion of a domain. The domain for intermediate values of ε is not uniquely
defined. However, linear interpolation between 0 and ε is in many cases the most convenient way of defining a
distortion mapping.
is commonly referred to as shape derivative:ω(Ω(ε)) − ω(Ω(0))
∂ω
= lim
.
ε→0
∂ε
ε

(3)

Here, Ω(ε) describes a family of domains smoothly governed by some parameter ε. An example of such a
distortion is depicted in Fig. 1.
Shape derivatives of eigenvalue problems have a long tradition. A classic example is the buckling of a column
problem (e.g. [3]). A general mathematical introduction to shape gradients and shape optimization can befound in [4]. Sensitivity of eigenvalue problems is discussed in [5]. In the recent years, the thermoacoustic
community has applied adjoint methods to compute the sensitivity of the eigenvalue with respect to many design
parameters [6, 7]. Only little attention is paid to the domain shape. For instance, Aguilar et al. used low-order
network models for modeling longitudinal modes, i.e., configurations that are essentially one-dimensional. In
these models the length of a tube appears as an explicit modeling parameter, such that shape may be treatedlike any other modeling parameter [8]. This approach can be extended to three dimensions by introducing
a parameterized discretization mesh. The parameters of which will then explicitly appear in the consideredequations. This is an idea that is used for instance in airfoil optimization using parameterized airfoil shapes[9]. However, the current study follows an approach which does not require to incorporate the considered
deformations on the level of a computational grid; in that sense the approach is independent of a parameterized
mesh. The adjoint formalism presented can be, hence, easily implemented in an existing code base. The
mathematics follows the derivation in [10], but is simplified to be independent of the mode normalizations.

2

THEORYIn general, a linear operator L consists of differential operators defined in a domain Ω, with corresponding
boundary conditions valid on ∂ Ω. In this report the notation(LΩ pb in Ω
L pb :=(4)L∂ Ω pb on ∂ Ωwill be used to distinguish between these two components of L.
Studying the eigenvalue problem

L(ω, ε) pb = 0(5)where ε is a parameter determining the domain shape, i.e., Ω = Ω(ε), is straightforward by means of perturbation theory. Introducing the ansätzeω = ω0 + ω1 ε + h. o. t. and pb = pb0 + pb1 ε + h. o. t.

(6)and the Taylor series expansion of L

L(ω, ε) = L(ω0 , 0) + ∂ω L(ω, 0) ∆ω + ∂ε L(ω0 , 0) ε + h. o. t.
| {z }| {z } | {z }:=L0,05490:=L1,0:=L0,1(7)

into the eigenvalue problem (5), yields(L0,0 + L1,0 ω1 ε + L0,1 ε) ( pb0 + pb1 ε) + h. o. t. = 0.

(8)

Sorting for like powers of ε gives the first-order equation

L0,0 pb1 + L1,0 pb0 ω1 + L0,1 pb0 = 0.

(9)This is the standard first-order equation for the eigenvalue problem. The unknowns ω1 and pb1 are found by
applying the inner product with the adjoint baseline solution pb†0 .
The adjoint operator of L is commonly defined to satisfy
ED
ha|LΩ bi = L†Ω a b + const. ∀a, b,(10)
where const. originates from boundary integrals that result from an integration by parts. This constant is calledboundary conjunct in [10] and denoted in the following as
const. := [a|b]LΩ .

(11)

The boundary conjunct is typically used to specify the boundary conditions L†∂ Ω for the adjoint equation by
demanding that the integration constant is vanishing. Such adjoint boundary conditions allow for a convenient
usage of the adjoint operator as subsequent calculations do not need to consider the constant. However, as thediscussion below shows, the boundary conjunct is required when the boundary itself is subject to a perturbation,i.e., displacements of the boundary are to be studied using perturbation methods.
2.1 Expansion of the boundary conditions at first order
To see the effect of boundary displacement, we assume that for all points on ∂ Ω(0) there is a mapping to∂ Ω(ε). This mapping is assumed to be one-to-one. Hence, locally we have for each point on ∂ Ω(0) a trajectory~x0 =~x0 (ε), such that ~x0 (ε) will be on ∂ Ω(ε). These trajectories are not unique. Yet, with this notation1 we can
expand a function value for a point ~x0 :pb(~x0 (ε)) = pb(~x0 (0)) + pb0 (~x0 (0))

d~x0
ε + h. o. t.
dε ε=0(12)

Using the expansion of a boundary value (12) at first order for a mixed boundary condition (2), we get thefollowing first-order expansion in terms of ε and ω:


 


d~x0
d~x
d~x0
0
0
0
0 =D(ω0 ) pb(~x0 (0)) + pb (~x0 (0))ε + ∇ pb(~x0 (0)) + pb (~x0 (0))ε · ~n(~x0 (0)) +~n (~x0 (0)) εdεdεdε+ D0 (ω0 ) pb(~x0 (0))∆ω.
(13)
At first order, this amounts to:
0 = D(ω0 ) pb(~x0 (0)) + ∇ pb(~x0 (0)) ·~n(~x0 (0))|{z}=L∂ Ω 0,0 pb


 



d~x0
d~x0
d~x0
0
0
0
+ D(ω0 ) pb (~x0 (0))+ ∇ pb (~x0 (0))·~n(~x0 (0)) + ∇ pb(~x0 (0)) · ~n (~x0 (0))
ε
dε
dε
dε
{z
}
|

(14)

=L∂ Ω 0,1 pb

0+ D (ω0 ) pb(~x0 (0)) ∆ω.
|{z}=L∂ Ω 1,0 pb

1 Please note the usage of ()0 for a derivative, which is considered to be a row vector in cases of spatial arguments.

5491

From (9) we find the boundary condition to the first-order correction pb1 (~x):
0 =D(ω0 ) pb1 (~x0 (0))∇ pb1 (~x0 (0)) ·~n(~x0 (0)) + D0 (ω0 ) pb0 (~x0 (0))ω1




 
d~x0
d~x
d~x0
·~n(~x0 (0)) + ∇ pb0 (~x0 (0)) · ~n0 (~x(0))+ ∇ pb00 (~x0 (0)).
+ D(ω0 ) pb00 (~x0 (0))
dε
dε
dε
|
{z
}

(15)

=r1 (~x0 (0))

All x0 are evaluated at ε = 0 which means this is a condition holding on the unperturbed surface:D(ω0 ) pb1 + ∇ pb1 ·~n = −r1 − ω1 D0 (ω0 ) pb0

on ∂ Ω(0)(16)

Because the perturbed boundary condition is formulated on the unperturbed surface, the approach presented inthis paper operates on the unperturbed domain only. This is the very reason why no parameterized mesh thatevolves with the perturbation parameter ε needs to be created. The deformation enters the problem statement interms of d~x/dε included in r1 .
2.2 The boundary conjunct of wave equations
For the thermoacoustic Helmholtz equation, the boundary conjunct arises from an integration by parts of the
Laplacian. It reads
I


(17)
[a|b]∇(·c2 ∇) =c2 a ∇b ·~n − ∇a ·~n b dS.
∂ Ω(0)

If a = pb† and b = pb, we findI
c2 pb† ∇ pb ·~n − ∇ pb† ·~n pb dS = −I † 
pb pb ∇(·c2 ∇) =∂ Ω(0)∂ Ω(0)


!
c2 pb† D(ω) pb − D† pb† pb dS = 0(18)

Here D† (ω) denotes the function determining the adjoint mixed-type boundary condition. In order to make theboundary conjunct vanish, the adjoint boundary condition needs to be specified by D† (ω) = D(ω).
With this specification, we find a particular formula for the boundary conjunct of pb†0 with pb1 :h

pb†0 pb1i∇(·c2 ∇)I=
∂ Ω(0)I=∂ Ω(0)I
c2 pb†0 ∇ pb1 ·~n − ∇ pb†0 ·~n pb1 dS
 I
c2 pb†0 ∇ pb1 ·~n − (−D(ω0 ) pb†0 ) pb1 =∂ Ω(0)∂ Ω(0)

c2 pb†0DE== − pb†0 r1

∂

0

I



−r1 − ω1 D (ω0 ) pb0 dS = −∂ Ω(0)
DE− ω1 pb†0 D0 (ω0 ) pb0 .

c2 pb†0 (∇ pb1 ·~n + D(ω0 ) pb1 )c2 pb†0 r1 dS − ω1I∂ Ω(0)(19)

c2 pb†0 D0 (ω0 ) pb0 dS∂

Note, that the boundary conjunct consists of surface integrals of quantities that are derived from the known
baseline solution pb0 and pb†0 . Hence, the computation of the boundary conjunct is numerically cheap.

5492

2.3 A first-order formula for shape sensitivity
Taking the inner product of (9) with the adjoint solution, we get
E
E
D
E D
D
pb†0 LΩ 0,0 pb1 + pb†0 LΩ 1,0 pb0 ω1 + pb†0 LΩ 0,1 pb0 = 0
*
⇒

+

(20a)

E
E
D
D
+ [ pb†0 | pb1 ]LΩ + pb†0 LΩ 1,0 pb0 ω1 + pb†0 LΩ 0,1 pb0 = 0

(20b)

E
E
D
E
D
D
− ω1 pb†0 D0 (ω0 ) pb0 + pb†0 LΩ 1,0 pb0 ω1 + pb†0 LΩ 0,1 pb0 = 0(20c)

L†Ω 0,0 pb†0 pb1| {z }=0−D

pb†0 r1

E∂∂

Unlike for perturbations of equation parameters, pb1 is not completely vanishing from the equation as it still
appears in the boundary conjunct [ pb†0 | pb1 ]LΩ . However, on the boundary ∂ Ω(0), the first-order drift pb1 can beinferred from the boundary displacement – see eq. (19) – and is thus known. Hence, the above equation allows
determining the eigenvalue drift at first order:E
E D
Dpb†0 LΩ 0,1 pb0 − pb†0 r1∂ E
E D
ω1 = − D
.
(21)
† 0
†
pb0 LΩ 1,0 pb0 − pb0 D (ω0 ) pb0∂

2.4 Implementational questions
The term r1 contains first and second order derivatives, which are to be integrated over the unperturbed domain
surface. In this study, we use second-order tetrahedral finite elements, which allow computing these quantities
on the triangles forming the discretized domain surface. For the Bubnov–Galerkin finite elements utilized,
the discrete and the continuous adjoint formulation, i.e. first deriving the discretization and then adjoining theoperators or vice versa, amount to the same discretization matrices for the adjoint equations. However, theutilized computational scheme applied to obtain the examples below is a continuous adjoint approach, because
the sensitivity is computed from the continuous formulation (21). Indeed, any approach utilizing the boundary
conjunct must be classified as continuous adjoint, as after the discretization the differential operators become
matrices and there would not be any concept of a boundary conjunct. This is opposed to the discrete adjoint
approaches using parameterized meshes to enable shape sensitivity analysis. An obvious major benefit fromthe continuous strategy is that there is no need to incorporate the deformations on the mesh-level, which would
require specialized software, and restrict the variety of applicable mesh types. We emphasize that the Rijke tube
example below uses an unstructured tetrahedral mesh, which has not been processed any further to conduct theanalysis, see Fig. 2.

3

ExampleAs an example, we consider a tubular duct that features a zone of heat release in its interior. Such Rijke-tube
configurations are commonly used as minimum working examples in the thermoacoustics community, e.g. in[1]. The Rijke tube studied here features a length of 0.5 m and a radius of R = 0.025 m. All boundaries except
the outlet are considered closed (i.e. soundhard)2 the outlet is open (i.e. sound soft). The speed of sound in thefirst half of the tube is set to 347 m/s; 695 m/s are specified in the second half of the tube. An n–τ-model
describes the flame dynamics
F(ω) = n exp(−iωτ)(22)where n is chosen to model the temperature jump for a medium with the same properties of air and τ is 1 ms.
The zone of heat release is set in the middle of the tube with the reference point being at the center axis
0.005 m upstream of the middle of the tube.
2 A sound-hard boundary condition is a simple model for sonic flow at the inlet.

5493

Figure 2. Computational grid used for the example. The lower row shows the inlet (left) and outlet (right)
surface.
Two eigenmodes, #M1 and #M2, are considered in this article; they are depicted in Fig. 3. As examples,
distortions of the inlet and the outlet surface are studied. Two different distortion mappings are considered. The
first one simply stretches the tube in the axial direction while the second one transforms the flat surface into aparaboloid. The two mappings are mathematically expressed as~x 7→ ~x + εR~ez(23)~x 7→ ~x + ε(R − (x2 + y2 ))~ez(24)

for the shift andfor the paraboloid. The two mappings are of comparable order, as the paraboloid is designed such that a pointon the center axis of the tube is displaced exactly in the same manner as when applying the planar shift.
Table 1 lists the shape sensitivities obtained from the adjoint method for various combinations of mode, distortion and considered surface. For comparison, shape sensitivities computed using a first-order forward finite
difference scheme are also listed in the table. To obtain these results, new meshes were generated from the oldones by slightly displacing the grid points on the considered surfaces according to the distortion mappings (23)and (24). The finite difference step size was ∆ε = 1E − 9. Because the displacements are defined to shift theconsidered surfaces in the direction of the z-axis, application to the inlet and the outlet will shorten and lengthen
the tube, respectively. This is in accordance with the computed sensitivities in the frequencies (real(ω)) and explains the opposite trends for the sensitivities obtained at the inlet and the outlet. Because the two displacement
types are similar at first order, their sensitivities are also very similar.
However, for all reported cases the finite difference and the adjoint method show excellent agreement. The
minor difference in the two methods can be attributed to a discretization error as the adjoint method discussed
here follows a continuous adjoint approach. Yet, the results validate the approach.

CONCLUSIONS
Adjoint-based computation of the shape gradient in a Rijke tube has been presented. Comparison of the resultswith finite-difference-approximations show good agreement. Compared to such traditional methods the adjoint-

5494

Figure 3. Unperturbed modes considered in the example.
Table 1. Comparison of adjoint-based results and finite difference estimates of the shape sensitivity in 1/s
case
#T1
#T2
#T3
#T4
#T5
#T6
#T7
#T8

domainoutletoutletinlet
inlet
outletoutletinlet
inlet

distortparaboloid
paraboloid
paraboloid
paraboloid
shift
shift
shift
shift

mode
#M1
#M2
#M1
#M2
#M1
#M2
#M1
#M2

adjoint−6.1556 − 4.054i−7.162 + 19.980i
6.554 + 7.103i
61.136 + 18.361i−6.236 − 4.107i−7.255 + 20.240i
6.638 + 7.195i
61.923 + 18.597i

finite difference−6.095 − 4.060i−7.138 + 19.926i
6.611 + 7.199i
61.904 + 18.598i−6.195 − 4.086i−7.218 + 20.187i
6.71 + 7.286i
62.753 + 18.888i

based formulation requires only one solution of the direct and one solution of the adjoint problem, thereby
considerably saving computational effort for the solution of nonlinear eigenvalue problems and costly remeshing.

Acknowledgments
A. Orchini thanks the Alexander von Humboldt Foundation for financial support via a Humboldt Research
Fellowship for Postdoctoral Researchers.

REFERENCES[1] Nicoud, F.; Benoit, L.; Sensiau, C.; Poinsot, T., Acoustic Modes in Combustors with Complex Impedancesand Multidimensional Active Flames, AIAA J. (45), 2007, 426–441.
[2] Güttel; S., Tisseur, F., The Nonlinear Eigenvalue Problem, Acta. Numer. (26), 2017, 1–94.
[3] Cox, S.J.; Overton, M.L., On the optimal design of columns against buckling, SIAM J. Math. Anal. (23),1992, 287–325.

5495

[4] Delfour, M.C.; Zolesio, J.-P., Shapes and geometries metrics, analysis, differential calculus, and optimization, Smith, R.C., Philiadelphia (USA), 2nd Ed., 2011.
[5] Kato, T., Perturbation Theory for Linear Operators, Springer, Heidelberg (Germany), 2nd Ed., 1980.
[6] Juniper, M.P., Sensitivity analysis of thermoacoustic instability with adjoint Helmholtz solvers, Phys. Rev.
Fluids (3), 2019, 110509.
[7] Magri, L, Adjoint methods as design tools in thermoacoustics, Appl. Mech. Rev. 71(2), 2019, 020801.
[8] Aguilar, J.; Magri, L.; Juniper, M.P., Adjoint-based sensitivity analysis of low-order thermoacoustic networks using a wave-based approach, J. Comput. Phys., 2017, 163–181.
[9] Giles, M. B.; Pierce, N. A., An Introduction to the Adjoint Apporach to Design, Flow Turbul. Combust(65), 200, 393-415[10] Parker, R.G.; Mote, C.D. Jr., Exact Boundary Condition Perturbation Solutions in Eigenvalue Problems, J.
Appl. Mech. (63), 1996, 128–135.

5496