Ekrem Ekici Department of Engineering, University of Cambridge, Trumpington Street, Cambridge, United Kingdom, CB21PZ email: ee331@cam.ac.uk Matthew P. Juniper1 Department of Engineering, University of Cambridge, Trumpington Street, Cambridge, United Kingdom, CB21PZ email: mpj1001@cam.ac.uk Shape Optimization of an Industrial Aeroengine Combustor to reduce Thermoacoustic Instability We parameterize the geometry of an industrial aeroengine combustor using Free Form Deformation (FFD). We then define the thermoacoustic system parameters and impose the acoustic boundary conditions to calculate the thermoacoustic eigenmode of the model using an open source parallelized Helmholtz solver. We then use adjoint methods to calculate the shape derivatives of the unstable eigenvalue with respect to the shape pa- rameters. First we calculate the sensitivities with respect to surface movements. Second, we calculate the sensitivities with respect to the FFD control points. We modify the FFD control point positions in order to reduce the thermoacoustic growth rate until the mode considered is stable. These findings show how, when combined with other constraints, this method could be used to reduce combustion instability in industrial annular combustors through geometric modifications. Keywords: thermoacoustic instability, design optimization, free form deformation, adjoint method, 1 Introduction Thermoacoustic oscillations endanger the safe operation of gas turbine engines. These oscillations are extremely sensitive to small changes in some system parameters, such as operating point, fuel composition and system geometry [1]. The sensitivities to the boundary conditions, mean temperature gradient, mean Mach number, and interaction index of the thermoacoustic system are examined in [2] through a parametric study of a one-dimensional Galerkin expansion model. Regarding the geometrical sensitivities, a database and design strategies for suppressing acoustic oscillations by changing the geometrical parameters are examined in [3] and [4] by focusing on the burner shape. In addition to the parametric and experimental studies, relatively cheap and accurate tools are desirable to examine the influence of the design parameters. To determine the sensitivity to system parameters, adjoint methods can be used [5]. With adjoint methods, the behaviour of the eigenvalue with respect to the system parameters can be calculated cheaply. In addition to the operating parameters, geometrical parameters of the combustor can be considered. The motivation is that small modifications to the geometry could lead to passive control of the system [6]. A low-order network model is used in [7] to stabilize the unstable modes of a longitudinal combustor by changing its shape, informed by adjoint methods. The same approach is applied to the annular combustor geometry in [8] by making small changes in radii, cross-sectional areas and lengths over the modules within the network model of the thermoacoustic system. However, low-order thermoacoustic network models are limited to relatively simple combustor geometries. For this reason, more complex geometries are better analysed with Helmholtz solvers, at the cost of increased computational effort. In Helmholtz solvers, the mean flow is assumed to be zero, and the fluctuating heat release rate is modelled as a distributed acoustic source. Adjoint methods give the sensitivity of the thermoacoustic eigenvalue to generic geometry changes. It is helpful to parameterize the geometry of a combustor and then to find the sensitivity to changes in those parameters. This creates a smaller and better posed optimization problem, with the loss of some generality. Helmholtz solvers are well-suited to handle complex combustor geometries, such as those in aircraft gas turbines. The parametrization of thermoacoustic system geometries were first examined for a two dimensional Rijke tube with B-Splines in [6]. Next, this method was extended to a relatively complex academic combustor, MICCA, with a Non- Uniform Rational B-Splines (NURBS) parametrization [9]. When dealing with realistic geometries, however, the tool needs to manage both local and global parametrization across several components. Using free form deformation (FFD) allows us to handle any geometry and allows us to control the degree of local or global parametrization complexity. This method is first introduced in [10] to increase the modelling capabilities and the representation accuracy of the surfaces of solid bodies in a free-form manner. With a free-form lattice formed by a few control points, the FFD technique provides the sensitivity to local and global deformations of the embedded geometry through control point displacements. In addition to computer graphics, the effectiveness of FFD in handling complex geometries has made it a strong candidate to improve structural designs in industry. One of the first examples of FFD is found in aerodynamic shape optimization problems in [11] and has recently been applied to airfoil geometries with adaptive parametrization techniques [12]. For turbomachinery applications, FFD has been combined with adjoints in order to optimize a jet-engine fan blade [13]. In gas turbine engines, the applications of FFD can be found for compressor [14] and turbine [15] designs but not for the combustor. In previous work, we applied FFD to relatively simple geometries [16] to reduce combustion instability. In this study, we extend FFD to an industrial aeroengine combustor geometry to reduce the thermoacoustic growth rate of the unstable eigenmode. We use an open-source adjoint Helmholtz solver, helmholtz-x, to calculate the sensitivity of the thermoacoustic eigenmodes 1Corresponding Author. Version 1.23, August 23, 2024 Journal of Engineering for Gas Turbines and Power Copyright © 2024 by ASME PREPRINT / 1 to changes in the FFD control points. We use the finite element method to discretize the system. The nonlinear eigenvalue problem is solved using fixed point iteration. The direct and adjoint eigenfunctions are used to calculate the shape derivatives of the FFD control points in Hadamard-form. Lastly, we change the geometry by moving the FFD control points in the directions informed by shape derivatives and eliminate the thermoacoustic instability for the specific unstable eigenmode. 2 Methodology 2.1 Thermoacoustic Helmholtz Equation. The derivation of the direct and adjoint thermoacoustic Helmholtz equations follows the methodology in [17]. The direct Helmholtz equation and momentum equation are ∇ · (︂ 𝑐2∇𝑝1 )︂ + 𝜔2𝑝1 = 𝑖𝜔(𝛾 − 1)�̂�1 (1a) −𝑖𝜌0𝜔û1 + ∇𝑝1 = 0, (1b) where 𝑐 is the spatially-varying speed of sound, 𝑝1 is the acoustic pressure, û1 is the acoustic velocity, 𝜔 is the complex valued angular frequency, 𝛾 is the heat capacity ratio, �̂�1 is any fluctuating heat release rate and 𝑝0 is the mean pressure. Eq. (1) can be written as L(𝜔)𝑝 = 0, where L is a differential operator that is linear in 𝑝1 but potentially nonlinear in 𝜔. The local 𝑛 − 𝜏 formulation is used to model the heat release rate caused by the unsteady flame response: 𝑞1 (x, 𝑡) �̄�0 = 𝑛ℎ(x) ∫ Ω 𝑤(x)û (x, 𝑡 − 𝜏(x)) · n𝑟 𝑑x �̄�𝑏 (2) where �̄�0 is the mean heat release rate, 𝑛 is the interaction index, �̄�𝑏 is the mean velocity, 𝜏 is the time delay and n𝑟 is the unit normal vector in the reference direction. ℎ(x) and 𝑤(x) are the heat release rate distribution and measurement function fields, which are modified Gaussian functions integrating to 1 over the domain. Introducing ⟨𝑝†1 |L𝑝1⟩ = 0, the adjoint Helmholtz equation and momentum equations are [17] ∇ · (︂ 𝑐2∇𝑝†1 )︂ + 𝜔∗2𝑝†1 = 𝑖𝜔∗ (𝛾 − 1)�̂�1 (𝜔∗) (3a) −𝑖𝜌0𝜔 ∗û1 + ∇𝑝†1 = 0. (3b) where 𝑝1 † is the adjoint eigenfunction and 𝑤∗ is the complex conjugate of the angular eigenfrequency. 2.2 Boundary Conditions. For annular combustors, most of the walls are Neumann (∇𝑝1 · n = 0) and Robin boundaries. Specific impedance, 𝑍 , can be imposed on Robin boundaries: ∇𝑝1 · n − 𝑖𝑤 𝑐𝑍 𝑝1 = 0 (4) for the direct solution and ∇𝑝†1 · n − 𝑖𝑤∗ 𝑐𝑍 𝑝 † 1 = 0 (5) for the adjoint solution. The reflection coefficient of the inlet choked boundary condition is 𝑅𝑖𝑛 = 1 − 𝛾𝑖 �̄�𝑖𝑛/(1 + (𝛾𝑖 − 1)�̄�2 𝑖𝑛 ) 1 + 𝛾𝑖𝑛�̄�𝑖𝑛/(1 + (𝛾𝑖𝑛 − 1)�̄�2 𝑖𝑛 ) , (6) where 𝛾𝑖𝑛 is the heat capacity ratio on the inlet choked boundary and �̄�𝑖𝑛 is the Mach number near the downstream of the inlet choked boundary. Similarly, we write the choked outlet condition as 𝑅𝑜𝑢𝑡 = 1 − (𝛾𝑜𝑢𝑡 − 1)�̄�𝑜𝑢𝑡/2 1 + (𝛾𝑜𝑢𝑡 − 1)�̄�𝑜𝑢𝑡/2 , (7) where 𝛾𝑜𝑢𝑡 is the heat capacity ratio on the outlet choked boundary and �̄�𝑜𝑢𝑡 is the Mach number near the upstream of the outlet choked boundary. These reflection coefficients can be imposed on Robin boundaries (Eq. (4) and Eq. (5)) through the specific impedance, 𝑍 = (1 + 𝑅)/(1 − 𝑅) [18]. 2 / PREPRINT Transactions of the ASME 2.3 Finite Element Formulation. Within the finite element framework, we integrate the terms in (1a) over the domain and multiply by a test function 𝑣 to obtain ∫ Ω ∇ · (︂ 𝑐2∇𝑝1 )︂ 𝑣 dx + ∫ Ω 𝜔2𝑝1𝑣 dx =∫ Ω 𝑖𝜔(𝛾 − 1)�̂�1𝑣 dx (8) We use trial functions 𝜙𝑖 such that 𝑝1 = ∑︁ 𝑖 𝜙𝑖 · 𝑝1,𝑖 and test functions 𝜙𝑗 such that 𝑣 = ∑︁ 𝑗 𝜙𝑗 . Integrating the terms in Eq. (8) by parts gives ∑︂ 𝑖 (︂∑︂ 𝑗 (︁ − ∫ Ω 𝑐2∇𝜙𝑖 · ∇𝜙𝑗 dx + ∫ 𝜕Ω 𝑐2∇𝜙𝑖 · n 𝜙𝑗 dS + ∫ Ω 𝜔2𝜙𝑖 𝜙𝑗 dx )︁ 𝑝1,𝑖 )︂ = ∑︂ 𝑖 (︂∑︂ 𝑗 (︁ ∫ Ω 𝑖𝜔(𝛾 − 1)�̂�1𝜙𝑗 dx )︁ )︂ , (9) where n is the normal vector of the relevant boundary. We can transform the second integral in Eq. (9) into the Robin integral using Eq. (1b). The matrix form of Eq. (9) is [︁ A + 𝜔B + 𝜔2C ]︁ p = D(𝜔)p (10) where A = − ∫ Ω 𝑐2∇𝜙𝑖 · ∇𝜙𝑗 dx, (11a) B = ∫ 𝜕Ω 𝑖𝑐 𝑍 𝜙𝑖 𝜙𝑗 dS, (11b) C = ∫ Ω 𝜙𝑖 𝜙𝑗 dx, (11c) D = FTF (𝛾 − 1) 𝑞0 𝑢𝑏 ∫ Ω 𝜙𝑗 ℎ(x) dx ∫ Ω 𝑤(x) 𝜌0 ∇𝜙𝑖 · n𝑟dx (11d) and p is the direct eigenvector. To derive the adjoint in matrix form, we take the conjugate transpose of Eq. (11) and calculate the right eigenvector to give the adjoint eigenvector. 2.4 Helmholtz solver. We implement an adjoint Helmholtz solver called helmholtz-x2. helmholtz-x uses Gmsh [19] for mesh generation, DOLFINx [20] for building the finite element model and UFL [21] for defining the weak forms. The sparse matrices are assembled with the PETSc package [22], and the PEP solver in the SLEPc package [23] is used to determine the nonlinear quadratic eigenvalue problem (Eq. (10)). The shift-and-invert spectral transformation is exploited to enhance the convergence of the eigenvalue to the targeted guess. Fixed point iteration with relaxation is implemented in order to converge to the eigenvalue [24]. The solver is parallelized with the OpenMPI library [25]. 3 Shape Parametrization The application of FFD to thermoacoustic instability in simple geometries is presented in a previous paper [16]. In this paper we investigate FFD to control thermoacoustic instability in industrial combustor geometries. In this section, we summarize the FFD method [10]. 3.1 Free Form Deformation. Free form deformation links mesh nodes to a handful of control points through trivariate Bernstein basis polynomials. The control points form the arbitrarily shaped control lattice (Fig. 1). In general, cylindrical or cubic control lattices are preferred in order to handle the control points more easily. Any mesh node within the control lattice can be represented by parametric coordinates (𝑠, 𝑡, 𝑢) as in Eq. (12), where X0 denotes the center of the FFD lattice and S, T and U represent the parametric unit vectors in the radial, circumferential and axial directions, respectively. X = X0 + 𝑠S + 𝑡T + 𝑢U (12) Considering the angular lattice in Fig. 1, the mesh nodes are initially transformed into cylindrical coordinates and then their parametric coordinates are calculated with Eq. (12). The range of the parametric coordinates is between 0 and 1 for the radial (𝑟) and axial (𝑧) directions and between 0 and 2𝜋 for the circumferential direction (𝜙). 2https://github.com/ekremekc/helmholtz-x/tree/JEGTP Journal of Engineering for Gas Turbines and Power PREPRINT / 3 https://github.com/ekremekc/helmholtz-x/tree/JEGTP Fig. 1 Free form deformation (FFD) configuration for a mesh using an angular control lattice with control points (red dots). r ,φ and z denotes the radial, circumferential and axial directions. Green lines visualize the connections between the control points. The position of the FFD control points can be arbitrarily defined depending on the application. In this paper, we specify the positions of the control points in a equispaced pattern within an angular lattice using Eq. (13). P𝑖 𝑗𝑘 = X0 + 𝑖 𝑙 S + 𝑗 𝑚 T + 𝑘 𝑛 U (13) where 𝑙, 𝑚, 𝑛 specify the number of control points in the 𝑟, 𝜙, 𝑧 axes. After specifying the positions of the FFD control points in the lattice and calculating Eq. (12), we are then ready to perform mesh deformation. We first displace the positions of the FFD control points and deform the mesh nodes individually, as shown in Eq. (14). X𝑓 𝑓 𝑑 = (︄ 𝑙∑︂ 𝑖=0 (︃ 𝑙 𝑖 )︃ (1 − 𝑠)𝑙−𝑖𝑠𝑖 (︂ 𝑚∑︂ 𝑗=0 (︃ 𝑚 𝑗 )︃ (1 − 𝑡)𝑚− 𝑗 𝑡 𝑗 (14) (︂ 𝑛∑︂ 𝑘=0 (︃ 𝑛 𝑘 )︃ (1 − 𝑢)𝑛−𝑘𝑢𝑘P𝑖 𝑗𝑘 )︂)︂)︄ In shape derivative calculations, the displacement field V𝑖 𝑗𝑘 is required for the control point P𝑖 𝑗𝑘 . We compute this by taking the derivative of the mesh nodes with respect to the control point, as shown in Eq. (15). 𝜕 𝜕P𝑖 𝑗𝑘 (︂ X𝑓 𝑓 𝑑 )︂ = V𝑖 𝑗𝑘 = (︄ 𝑙∑︂ 𝑖=0 (︃ 𝑙 𝑖 )︃ (1 − 𝑠)𝑙−𝑖𝑠𝑖 (︂ 𝑚∑︂ 𝑗=0 (︃ 𝑚 𝑗 )︃ (1 − 𝑡)𝑚− 𝑗 𝑡 𝑗 (15) (︂ 𝑛∑︂ 𝑘=0 (︃ 𝑛 𝑘 )︃ (1 − 𝑢)𝑛−𝑘𝑢𝑘 )︂)︂)︄ Then we input V𝑖 𝑗𝑘 into Eq. (16) to calculate the shape derivative of P𝑖 𝑗𝑘 . 3.2 Shape Derivatives. The shape sensitivities for the thermoacoustic Helmholtz equation are introduced in [17]. With direct and adjoint eigenvectors, we can calculate the shape derivative of the FFD control points in Hadamard-form. The most general expression for the shape derivative is with Robin boundary conditions: 𝜔′ 𝑖 𝑗𝑘 = ∫ Γ V𝑖 𝑗𝑘 · n𝑖 𝑗𝑘 (︄ − 𝑝1 †∗ (︂ 𝜅𝑐2 𝜕𝑐 𝜕𝑛 )︂ 𝜕𝑝1 𝜕𝑛 + ∇ · (︂ 𝑝1 †∗𝑐2∇𝑝1 )︂ 4 / PREPRINT Transactions of the ASME − 2 𝜕𝑝1 †∗ 𝜕𝑛 𝑐2 𝜕𝑝1 𝜕𝑛 )︄ 𝑑𝑆 (16) where 𝜔′ 𝑖 𝑗𝑘 is the complex-numbered shape derivative for the control point P𝑖 𝑗𝑘 , and n𝑖 𝑗𝑘 is its outward normal vector. We consider the FFD control points on Neumann boundaries and impose 𝜕𝑝1/𝜕𝑛 = 0 and 𝜕𝑝 † 1/𝜕𝑛 = 0. We compute gradient information in all three directions for the pertinent control point. The influence of the FFD control point movement on the stability of the system is given by the imaginary component of the shape gradient. In cases where the computed mode is unstable, displacing the control point in the opposite direction of the imaginary component of the shape derivative will make the system more stable. The steps of the shape optimization procedure using free form deformation can be summarized as: • The three dimensional mesh is generated. • The FFD lattice and control points are defined after calculation of the parametric coordinates of the mesh nodes. • Direct and adjoint eigenmodes are calculated with degree 2 finite elements. • The shape derivatives of the FFD control points are calculated and normalized. • The shape is deformed following the direction provided by the shape derivatives, with a certain step size. • The eigenmode for the deformed geometry is computed with degree 1 finite elements to see the outcome of the deformation. 4 FFD Geometry and System Parameters 4.1 Geometry and FFD Setup. We start by simplifying the CFD geometry of the industrial annular combustor by removing cooling holes on the liner walls. In this paper, we only deal with the axial eigenmodes. For that reason, we only consider a single sector of the full annulus. We generate 384, 785 finite elements using Delaunay-triangulation. The global FFD setup for the combustor geometry can be seen in Fig. 2. We place 2, 3 and 4 FFD control points in the radial, circumferential and axial directions, respectively. Fig. 2 FFD configuration of the simplified aeroengine combustor geometry. The red dots show the locations of the FFD control points of the angular lattice. The torus-like geometry in the combustion chamber represents the volumetric heat release rate field, h (x). The choked boundaries are also shown. 4.2 Parameters. The parameters of the thermoacoustic problem are tabulated in Table 1. The imposed volumetric heat release rate field is shown in Fig. 2. This field integrates to 1 over the torus domain. The measurement function distribution (Fig. 3(a)) is implemented by means of a three dimensional Gaussian function in Eq. (17), where 𝛼 controls the width of the function around point 𝑥𝑟 (𝑥0, 𝑦0, 𝑧0). 𝐺 (x) = exp (︂ − (︁ (𝑥 − 𝑥0)2 + (𝑦 − 𝑦0)2 + (𝑧 − 𝑧0)2 )︁ /(2𝛼2) )︂ 𝛼3 (2𝜋)3/2 (17) Journal of Engineering for Gas Turbines and Power PREPRINT / 5 Table 1 Parameters of the thermoacoustic problem for the annular combustor. Parameter value unit 𝑟𝑔𝑎𝑠 287.0 Jkg−1K−1 𝑝𝑔𝑎𝑠 334.809 kPa �̄�𝑡𝑜𝑡 -8.19568E6 W �̄�𝑏 90.98 m s−1 𝑛 -6.912 - 𝜏 0.0039 s 𝛼 0.004 - 𝑀𝑖𝑛 0.03074 - 𝑀𝑜𝑢𝑡 0.0702 - The relationship between the heat release rate and acoustic velocity at the measurement point is defined through an 𝑛−𝜏 model [24]. The temperature data of the combustor is taken from the network model and extrapolated onto the 3D numerical grid. The speed of sound field changes in the axial direction within the combustion chamber. Using the temperature field shown in Fig. 3(b), the distribution of the speed of sound can be determined through the equation 𝑐 = √︁ 𝛾𝑟𝑔𝑎𝑠𝑇 . (a) Measurement function field. (b) Non-dimensional temperature field. Fig. 3 The (a) measurement function field and (b) non-dimensional temperature field of the combustor. The temperature field is obtained from a low order network code. The specific heat capacity is assumed to change linearly with the temperature, 𝑐𝑝 (𝑇) = 973.60 + 0.133𝑇 . We compute the mean density field with the ideal gas model 𝑝𝑔𝑎𝑠 = 𝜌0𝑟𝑔𝑎𝑠𝑇 . We only consider one sector, so the heat release for this sector is �̄�𝑠𝑒𝑐 = �̄�𝑡𝑜𝑡/20. Choked boundary conditions are imposed through the reflection coefficients for the inlet and outlet boundaries using Eqs. (6) and (7). The other boundaries are assumed to be perfectly reflecting (Neumann) surfaces. 5 Results 5.1 Eigenmodes. Using this geometry, we show two unstable and one stable axial eigenmodes in Fig. 4. In this paper, we only calculate shape derivatives of the most unstable eigenmode, shown in Fig. 4(a). We first compute the eigenfrequency in Fig 4(a) without 6 / PREPRINT Transactions of the ASME the flame response and call it 𝑓𝑟𝑒 𝑓 . We then use this to divide the eigenvalues. (a) f = 1.2056 + 0.0321i (b) f = 1.5189 + 0.0006i (c) f = 1.8014 − 0.0014i Fig. 4 Various axial eigenmodes of the combustor with the corresponding eigenfrequencies scaled with fr ef . 5.2 Surface Sensitivities. We start by obtaining the shape sensitivities of the boundaries for the axial thermoacoustic mode in Fig. 4(a). For this, we divide the shape derivative by the surface area of the boundary to obtain local average [17]. As a result, we modify Eq.(16) to obtain 𝜔′ 𝑠 = ∫ Γ 1 𝐴𝑠 (︄ − 𝑝1 †∗ (︂ 𝜅𝑐2 𝜕𝑐 𝜕𝑛 )︂ 𝜕𝑝1 𝜕𝑛 + ∇ · (︂ 𝑝1 †∗𝑐2∇𝑝1 )︂ − 2 𝜕𝑝1 †∗ 𝜕𝑛 𝑐2 𝜕𝑝1 𝜕𝑛 )︄ 𝑑𝑆𝑠 (18) Journal of Engineering for Gas Turbines and Power PREPRINT / 7 where subscript 𝑠 denotes the corresponding surface. We use Eq. (18) to calculate the average shape derivatives of the combustor surfaces. We then normalize them by the surface with the largest imaginary part. The imaginary parts of the computed surface sensitivities are visualized in Fig. 5. We do not show the surface sensitivities for surfaces with areas less than 1.3 ×10−3 m2 because the sensitivity per unit area is too large to be seen on this scale. (a) (b) Fig. 5 Surface sensitivities for the unstable eigenmode from (a) front and (b) rear views. Inspecting Fig. 5, the surfaces near the injector are found to be the most influential surfaces for control of the thermoacoustic growth rate. In addition, the inlet and outlet boundaries should be moved in the outward normal direction to achieve a more stable axial eigenmode. 5.3 Control Point Sensitivities. As we consider Neumann boundary conditions for FFD control points, Eq. (16) simplifies to 𝜔′ 𝑖 𝑗𝑘 = ∫ Γ (︁ V𝑖 𝑗𝑘 · e𝑖 𝑗𝑘 )︁ ∇ · (︂ 𝑝1 †∗𝑐2∇𝑝1 )︂ 𝑑𝑆. (19) where e𝑖 𝑗𝑘 represents the unit vector. We present an example displacement field for the control point P1,1,0 in Fig. 6. We use Eq. (19) to calculate the FFD control points’ shape sensitivities. We then normalize them as in section 5.2. The imaginary components of the normalized shape derivatives are shown in Fig. 7. We see that the FFD shape derivatives show that the thermoacoustic growth rate can be reduced by enlarging the combustor geometry in the axial direction. This is consistent with the area derivatives presented in Fig. 5. 8 / PREPRINT Transactions of the ASME Fig. 6 An example displacement field for the FFD control point P1,1,0. We use Eq. (15) to compute this. 5.4 Mesh Deformation. We now apply free form deformation to the entire geometry. We update the positions of the FFD control points in the directions provided by the imaginary part of the shape derivatives. We deform the mesh by a certain step size equal to 1/40𝑡ℎ of the combustor length. We then update the location of the measurement point 𝑥𝑟 to be positioned at a given distance from the injector end. We also adapt the temperature field such that it fits inside the combustion chamber after deformation. For the deformed geometry, we obtain the eigenfrequency of 𝑓 = 1.1728 − 0.0023𝑖 for the same axial eigenmode. This shows that this free form deformation has reduced the growth rate and thereby stabilized the unstable axial eigenmode. The frequency of the deformed system is reduced because the combustor length is increased in the axial direction. The slice view for the optimized geometry is presented in Fig. 8. 6 Conclusions We have proposed a shape optimization procedure to reduce thermoacoustic instability in an industrial gas turbine combustor geometry. We use the open-source finite element thermoacoustic solver, helmholtz-x. We assume an 𝑛 − 𝜏 formulation for modelling the flame response and we input the mean flow data from an industrial low order network code. We parallelize the calculation of shape derivatives for this industrial geometry. We then parametrize the geometry with free form deformation. We compute the direct and adjoint eigenmodes with degree 2 finite elements. We initially calculate the influence of surface deformations on the unstable eigenmode and calculate the averaged shape derivatives of the boundaries. These derivatives provide some insight as to the deformation directions that will stabilize the thermoacoustic eigenmodes. We then calculate the derivative of the unstable eigenvalue with respect to the positions of FFD control points. Following these gradients, we deform the shape using the FFD control points to reduce the growth rate of the unstable eigenmode. With these deformations, we suppress thermoacoustic instability of this axial mode. The proposed shape changes through the FFD control points agreed with the averaged shape derivatives. The results in the paper show that adjoint-based shape sensitivity can be combined with free-form deformation to stabilize thermoa- coustic modes in industrial combustion chamber geometries. In practice, combustion chambers and plenums have other geometrical constraints that are not considered in this paper. These constraints could, however, be included within the optimization algorithm, while retaining the attractive features of gradient-based optimization outlined in this paper. In conclusion, this paper shows that adjoint Helmholtz solvers can be a useful design tool for optimizing detailed combustor geometries through the FFD technique. Next steps are to improve the accuracy of the components of the Helmholtz solver. Currently, we neglect the acoustic liners and we have removed the swirler geometry within the injector. These influence the thermoacoustic behaviour and have been modelled by another Helmholtz solver [26, 27]. We will implement these within our finite element framework. We will also include experimentally-derived flame-transfer function data. Acknowledgment Ekrem Ekici gratefully acknowledges funding from Türkiye’s Ministry of National Education. The authors would like to thank Simon Stow for supplying the network code and mean flow data. Data availability The open-source helmholtz-x is accessible in the GitHub repository link; https://github.com/ekremekc/helmholtz-x/tree/JEGTP. Nomenclature Journal of Engineering for Gas Turbines and Power PREPRINT / 9 https://github.com/ekremekc/helmholtz-x/tree/JEGTP Roman letters 𝑐 = Speed of sound [m s−1] 𝑢𝑏 = Mean velocity [m s−1] 𝑛 = Non-dimensional interaction index [-] 𝑞𝑡𝑜𝑡 = Total mean heat release rate [W] 𝑝 = Acoustic pressure [N m−2] 𝑝𝑔𝑎𝑠 = Mean pressure [N m−2] û = Acoustic velocity [m s−1] �̄� = Mach number [-] 𝑇 = Temperature [K] 𝑍 = Specific impedance [-] 𝑅 = Reflection coefficient [-] ℎ = Mean heat release rate field [W] 𝑤 = Measurement function field [-] 𝑠 = Parametric coordinate in the radial direction [-] 𝑡 = Parametric coordinate in the circumferential direction [-] 𝑢 = Parametric coordinate in the axial direction [-] 𝑃 = Three dimensional FFD control point vector [-] 𝑋 = Three dimensional mesh node vector [-] 𝑟𝑔𝑎𝑠 = Universal gas constant [J/kg K−1] 𝑥𝑟 = Center of the measurement function field [m] 𝑓 = Eigenfrequency [s−1] Greek letters 𝜔 = Angular eigenfrequency [rad s−1] 𝛾 = Heat capacity ratio [-] 𝜏 = Time delay [s−1] 𝜅 = Curvature [-] 𝛼 = Standart deviation of Gaussian function [-] 𝜙𝑖 = Trial function [-] 𝜙𝑗 = Test function [-] Superscripts and subscripts ()̂1 = perturbed quantities (̂)† = adjoint variables ()∗ = complex conjugate ()0 = mean quantities ()𝑏 = bulk quantities ()𝑢 = upstream ()𝑠 = index of the surface ()𝑑 = downstream ()𝑖𝑛 = inlet ()𝑜𝑢𝑡 = outlet ()𝑡𝑜𝑡 = total ()𝑠𝑒𝑐 = sector ()𝑟𝑒 𝑓 = reference ()𝑖 = the index of the control point in the radial direction ()𝑗 = the index of the control point in the circumferential direction ()𝑘 = the index of the control point in the axial direction References [1] Juniper, M. P., 2018, “Sensitivity analysis of thermoacoustic instability with adjoint Helmholtz solvers,” Physical Review Fluids, 3(11), p. 110509. [2] Zhao, X., Zhao, D., Cheng, L., Shelton, C. M., and Majdalani, J., 2023, “Predicting thermoacoustic stability characteristics of longitudinal combustors using different endpoint conditions with a low Mach number flow,” Physics of Fluids, 35(9), pp. –. [3] Bade, S., Wagner, M., Hirsch, C., Sattelmayer, T., and Schuermans, B., 2013, “Design for thermo-acoustic stability: modeling of burner and flame dynamics,” Journal of engineering for gas turbines and power, 135(11), p. 111502. [4] Bade, S., Wagner, M., Hirsch, C., Sattelmayer, T., and Schuermans, B., 2013, “Design for thermo-acoustic stability: Procedure and database,” Journal of engineering for gas turbines and power, 135(12), p. 121507. [5] Magri, L., 2019, “Adjoint methods as design tools in thermoacoustics,” Applied Mechanics Reviews, 71(2), p. 020801. [6] Falco, S. and Juniper, M. P., 2021, “Shape optimization of thermoacoustic systems using a two-dimensional adjoint helmholtz solver,” Journal of Engineering for Gas Turbines and Power, 143(7), p. 071025. [7] Aguilar, J. G. and Juniper, M. P., 2020, “Thermoacoustic stabilization of a longitudinal combustor using adjoint methods,” Physical Review Fluids, 5(8), p. 083902. [8] Aguilar, J. G. and Juniper, M. P., 2018, “Adjoint methods for elimination of thermoacoustic oscillations in a model annular combustor via small geometry modifications,” Turbo Expo: Power for Land, Sea, and Air, Vol. 51050, American Society of Mechanical Engineers, , , p. V04AT04A054. [9] Ekici, E., Falco, S., and Juniper, M. P., 2024, “Shape sensitivity of thermoacoustic oscillations in an annular combustor with a 3D adjoint Helmholtz solver,” Computer Methods in Applied Mechanics and Engineering, 418, p. 116572. [10] Sederberg, T. W. and Parry, S. R., 1986, “Free-form deformation of solid geometric models,” Proceedings of the 13th annual conference on Computer graphics and interactive techniques, , , pp. 151–160. [11] Samareh, J., 2004, “Aerodynamic shape optimization based on free-form deformation,” 10th AIAA/ISSMO multidisciplinary analysis and optimization conference, , , p. 4630. 10 / PREPRINT Transactions of the ASME [12] He, X., Li, J., Mader, C. A., Yildirim, A., and Martins, J. R., 2019, “Robust aerodynamic shape optimization—from a circle to an airfoil,” Aerospace Science and Technology, 87, pp. 48–61. [13] Giugno, A., Shahpar, S., and Traverso, A., 2020, “Adjoint-based optimization of a modern jet-engine fan blade,” Turbo Expo: Power for Land, Sea, and Air, Vol. 84096, American Society of Mechanical Engineers, , , p. V02DT38A026. [14] John, A., Shahpar, S., and Qin, N., 2017, “Novel compressor blade shaping through a free-form method,” Journal of Turbomachinery, 139(8), p. 081002. [15] Li, L., Yuan, T., Li, Y., Yang, W., and Kang, J., 2019, “Multidisciplinary design optimization based on parameterized free-form deformation for single turbine,” AIAA Journal, 57(5), pp. 2075–2087. [16] Ekici, E. and Juniper, M., 2024, “Adjoint based shape optimization for thermoacoustic stability of combustors using free form deformation,” Turbo Expo, Apollo - University of Cambridge Repository, , , p. , doi: 10.17863/CAM.107364. [17] Falco, S., 2022, “Shape optimization for thermoacoustic instability with an adjoint Helmholtz solver,” Ph.D. thesis, University of Cambridge. [18] Rienstra, S. W. and Hirschberg, A., 2004, “An introduction to acoustics,” Eindhoven University of Technology, 18, p. 19. [19] Geuzaine, C. and Remacle, J.-F., 2009, “Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities,” International journal for numerical methods in engineering, 79(11), pp. 1309–1331. [20] Barrata, I. A., Dean, J. P., Dokken, J. S., HABERA, M., HALE, J., Richardson, C., Rognes, M. E., Scroggs, M. W., Sime, N., and Wells, G. N., 2023, “DOLFINx: The next generation FEniCS problem solving environment,” , p. . [21] Alnæs, M. S., Logg, A., Ølgaard, K. B., Rognes, M. E., and Wells, G. N., 2014, “Unified form language: A domain-specific language for weak formulations of partial differential equations,” ACM Transactions on Mathematical Software (TOMS), 40(2), pp. 1–37. [22] Balay, S., Abhyankar, S., Adams, M., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Dener, A., Eijkhout, V., Gropp, W., et al., 2019, “PETSc users manual,” , p. . [23] Hernandez, V., Roman, J. E., and Vidal, V., 2005, “SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems,” ACM Transactions on Mathematical Software (TOMS), 31(3), pp. 351–362. [24] Nicoud, F., Benoit, L., Sensiau, C., and Poinsot, T., 2007, “Acoustic modes in combustors with complex impedances and multidimensional active flames,” AIAA journal, 45(2), pp. 426–441. [25] Gropp, W., Lusk, E., Doss, N., and Skjellum, A., 1996, “A high-performance, portable implementation of the MPI message passing interface standard,” Parallel computing, 22(6), pp. 789–828. [26] Ni, F., Miguel-Brebion, M., Nicoud, F., and Poinsot, T., 2017, “Accounting for acoustic damping in a Helmholtz solver,” AIAA Journal, 55(4), pp. 1205–1220. [27] Ni, F., Nicoud, F., Méry, Y., and Staffelbach, G., 2018, “Including flow–acoustic interactions in the Helmholtz computations of industrial combustors,” AIAA Journal, 56(12), pp. 4815–4829. [28] Morgans, A. S. and Stow, S. R., 2007, “Model-based control of combustion instabilities in annular combustors,” Combustion and flame, 150(4), pp. 380–399. Journal of Engineering for Gas Turbines and Power PREPRINT / 11 https://doi.org/10.17863/CAM.107364 (a) (b) (c) Fig. 7 Imaginary parts of the FFD shape derivatives for the eigenfrequency 1.2056 + 0.0321i . The direction of the derivatives show the direction of changes for stabilizing the system. 12 / PREPRINT Transactions of the ASME Fig. 8 The initial (gray) and deformed (green) geometries with initial (red dots) and final (green dots) positions of the FFD control points. The eigenfrequency of the system changed from f = 1.2056 + 0.0321i (gray) to f = 1.1728 − 0.0023 (green). Journal of Engineering for Gas Turbines and Power PREPRINT / 13 List of Figures 1 Free form deformation (FFD) configuration for a mesh using an angular control lattice with control points (red dots). 𝑟, 𝜙 and 𝑧 denotes the radial, circumferential and axial directions. Green lines visualize the connections between the control points. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 FFD configuration of the simplified aeroengine combustor geometry. The red dots show the locations of the FFD control points of the angular lattice. The torus-like geometry in the combustion chamber represents the volumetric heat release rate field, ℎ(x). The choked boundaries are also shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3 The (a) measurement function field and (b) non-dimensional temperature field of the combustor. The temperature field is obtained from a low order network code. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 (a) Measurement function field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 (b) Non-dimensional temperature field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 4 Various axial eigenmodes of the combustor with the corresponding eigenfrequencies scaled with 𝑓𝑟𝑒 𝑓 . . . . . . . . . . . 7 (a) 𝑓 = 1.2056 + 0.0321𝑖 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 (b) 𝑓 = 1.5189 + 0.0006𝑖 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 (c) 𝑓 = 1.8014 − 0.0014𝑖 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 5 Surface sensitivities for the unstable eigenmode from (a) front and (b) rear views. . . . . . . . . . . . . . . . . . . . . . 8 (a) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 6 An example displacement field for the FFD control point P1,1,0. We use Eq. (15) to compute this. . . . . . . . . . . . . 9 7 Imaginary parts of the FFD shape derivatives for the eigenfrequency 1.2056 + 0.0321𝑖. The direction of the derivatives show the direction of changes for stabilizing the system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 (a) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 (c) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 8 The initial (gray) and deformed (green) geometries with initial (red dots) and final (green dots) positions of the FFD control points. The eigenfrequency of the system changed from 𝑓 = 1.2056 + 0.0321𝑖 (gray) to 𝑓 = 1.1728 − 0.0023 (green). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 List of Tables 1 Parameters of the thermoacoustic problem for the annular combustor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 14 / PREPRINT Transactions of the ASME 1 Introduction 2 Methodology 2.1 Thermoacoustic Helmholtz Equation 2.2 Boundary Conditions 2.3 Finite Element Formulation 2.4 Helmholtz solver 3 Shape Parametrization 3.1 Free Form Deformation 3.2 Shape Derivatives 4 FFD Geometry and System Parameters 4.1 Geometry and FFD Setup 4.2 Parameters 5 Results 5.1 Eigenmodes 5.2 Surface Sensitivities 5.3 Control Point Sensitivities 5.4 Mesh Deformation 6 Conclusions Appendices