Three-Dimensional Analytical Solutions for Acoustic Transverse Modes in a Cylindrical Duct with Axial Temperature Gradient and Non-Zero Mach Number

: Cylindrical ducts with axial mean temperature gradient and mean ﬂows are typical elements in rocket engines, can combustors, and afterburners. Accurate analytical solutions for the acoustic waves of the longitudinal and transverse modes within these ducts can signiﬁcantly improve the performance of low order acoustic network models for analyses of acoustic behaviours and combustion instabilities in these kinds of ducts. Here, we derive an acoustic wave equation as a function of pressure perturbation based on the linearised Euler equations (LEEs), and the modiﬁed WKB approximation method is applied to derive analytical solutions based on very few assumptions. The eigenvalue system is built based on the proposed solutions and applied to predict the resonant frequencies and growth rate for transverse modes. Validations of the proposed solutions are performed by comparing them to the numerical results directly calculated from the LEEs. Good agreements are found between analytical reconstruction and numerical results of three-dimensional transverse modes. The system with both mean temperature proﬁle and mean ﬂow presents a larger absolute value of the growth rate than the condition of either uniform mean temperature or no mean ﬂow.


Introduction
Combustion instabilities are typically present in rocket engines, aero-engines and land-used gas turbine engines, and may cause severe damage due to coupling between unsteady heat release from combustion and the acoustic system within the combustors or even the entire engine [1,2]. It is thus significant to predict the combustion instabilities and optimise the engine geometries and operating conditions to eliminate these instabilities during the design stage of engines [3][4][5].
There are typically two methods to numerically predict and analyse combustion instabilities [5,6]. The first is to directly simulate the coupling mechanisms of the unsteady flow, combustion, and acoustic dynamics based on the complete three-dimensional compressible Computational Fluid Dynamics (CFD) simulations;the compressible large eddy simulation (LES) is typically preferred among these CFD solvers [5,7]. Although great achievements have been made in the development of the LES, this approach remains very expensive and time-consuming, and is difficult to use in the real engine design process [5,8,9]. The second method is to decouple the calculation of the unsteady heat release from the flame and acoustic system [10][11][12]; the first term is characterised by a flame transfer function for linear analysis or a flame describing function for weakly nonlinear analysis, see e.g., the recent reviews [13,14]. The generation, propagation, transmission, and reflection of acoustic waves, or even the entropy and vorticity waves within the combustor with complex geometries, are characterised by low order acoustic network models [15][16][17][18][19][20][21], linearised Euler mass and flux, and conserved mass, momentum, and energy, the Euler equations and equation of state are expressed as follows: where ρ, u = (u x , u θ , u r ), p denote the density, velocity components, and pressure, respectively,q represents heat flux, and the heat capacity ratio γ and universal gas constant R g are assumed to be constant along the cylindrical duct. To linearise the Euler equations, the primitive variables are assumed to be the sum of a time-average steady component that is denoted by() and a small time-varying perturbation in the form of time and angular periodic dependence, expressed as a =ā +â(x, r) exp(iωt + in θ θ), where ω = 2π f + iω i is the complex angular frequency and the circumferential wavenumber n θ is definitely an integer due to the continuity in the θ direction. It should be noted that ω i < 0 denotes that the corresponding thermoacoustic mode is unstable.
The conservation relations for the steady flow variables are deduced, and can be written as follows: 1 With a prescribed mean temperature profileT(x), the consequent inhomogeneous mean flow field can be calculated along the axial direction. For greater convenience of subsequent derivation, the derivatives of the steady flow variables are expressed in the form of the relative density derivative and steady variables, written as follows: wherec = (γR gT ) 1/2 is the local sound speed.
Then, the first derivatives of the steady temperature, speed of sound, mean flow Mach number M x =ū x /c, and local wave number k 0 = ω/c are expressed as follows: It is evident that the axial temperature gradient is associated with the steady heat flux in the presence of non-zero mean flow. Therefore, the mean temperature gradient can be maintained by steady heat addition or extraction when the mean flow propagates through the three-dimensional cylindrical duct.
For the perturbation variables, the linearised Euler equations (LEEs) of the energy and momentum can be obtained: whereŝ denotes the entropy perturbation, which is introduced by the relationŝ/c p =p/γp − ρ/ρ to replace the density perturbation in the axial momentum LEE by the pressure perturbation and entropy wave. Here, c p = γR g /(γ − 1) is the heat capacity at constant pressure. With the unsteady heat fluxq set to zero, only the steady heat communication is taken into account. Considering the coupling term of −ū x (dū x /dx)(ŝ/c p ) in Equation (12), the entropy wave is generated during the acoustic wave travelling through the inhomogeneous mean temperature zone and in return affects the acoustic field in the presence of non-uniform mean velocity. It has been validated that the sound regeneration of the entropy wave convected by the mean flow can be neglected, especially in the high-frequency domain [43,55]. Therefore, the axial momentum LEE (Equation (12)) is further simplified by neglecting the acoustic-entropy coupling term, provided by A pure acoustic problem is derived, facilitating the acoustic wave equation and analytical solutions.
Assuming that the wave number is sufficiently larger than the critical value, |k 0 | |M x α|, and terms with order higher than M 2 x can be ignored, Equation (15) can be divided by (iω + dū x /dx) to eliminate the term (dp/dx)û x in Equation (11), leading to Then, Equation (15) is divided byū x , expressed as Subsequently, Equation (16) and Equation (17) are manipulated according to following formula: Following the above procedure, all terms withû x on the left-hand side (LHS) of the consequent formula are eliminated, and the termsû θ andû r on the right-hand side (RHS) can be replaced with the terms of the pressure perturbation in the linearised radial and azimuthal momentum equations (Equations (13) and (14)), written as follows: where β = d 2ρ /dx 2 /ρ. Dividing both sides by the x-dependent coefficient (1 − M x α/ik 0 ) 2 and assuming that |k 0 | |α| and |k 2 0 | |β|/2, we finally obtain the wave equation in the form of pressure perturbation only, expressed as follows: The general PDE in the form of pressure perturbation is ultimately obtained from the governed equations with the assumption of a sufficiently large wave number and the negligible effect of entropy perturbation in such a high-frequency domain. This allows for the derivation of the analytical expression of the pressure perturbation in the next section and the theoretical reconstruction of three-dimensional acoustic waves in the inhomogeneous background mean field.

Analytical Solutions of Three-Dimensional Acoustic Wave
The partial differential wave equation (Equation (20)) is characterized by partial derivatives and coefficients concerning x and r on the LHS and RHS, respectively. Therefore, separation of variables can be applied by substitutingp(x, r) = X (x)R(r) into the wave equation, resulting in the radial and axial ordinary differential equations (ODEs): and where λ is a constant that is independent of the radial and axial directions. For the radial component, the solution for the Bessel equation can be written as the superposition of two linearly independent Bessel functions: where J n θ and Y n θ are the Bessel functions of the first and second kind. Note that the coefficient c 2 must be 0, as Y n θ diverges as r tends to 0. With the rigid-wall boundary condition at r = r 0 , a series of discrete solutions for the constant λ is calculated and corresponds to different acoustic transverse modes labelled (n θ , n r ), where n r denotes the (n r + 1)th solution for each value of n θ .

Analytical Solutions of the Axial ODE Using a Modified WKB Approximation Method
For the axial component, no known solutions can directly deal with spatially varying coefficients of the axial ODE (Equation (22)), especially in the presence of both the axially arbitrary temperature profile and mean flow. The modified WKB approximation method is then used to resolve the axial ODE, which is essentially based on the large wave number assumption and successfully applied on the one-dimensional acoustic wave equation within varying coefficients [43].
The modified WKB solution uses the form of the separate amplitude and phase factors to represent the properties of the acoustic wave propagating through the axially varying mean field, expressed as follows: where C is a constant coefficient and a and b are x-dependent real variables. The assumption of a larger wave number naturally results in the precondition of |b| |a| within the modified WKB solution.
To obtain real variables a and b, the locally complex wave number k 0 = ω/c = k r + ik i has to be substituted into Equation (22) by assuming |k r | |M x k i |, provided by Then, the axial ODE can be transformed into an algebraic equation using the modified WKB solution, written as The real part is further simplified by combining previous assumptions for the larger wave number, |k 0 | |α|, |b| |a| and |k r | |M x k i |, expressed as: The solutions of b are directly written as Hence, the cut-off frequency f c is provided by It should be noted that this work only considers the condition ω > ω c in order to make sure that the formula under the square root is always positive, as acoustic waves often rapidly dissipate with axial distance when ω < ω c .
Solutions of a are further derived by substitute b ± into the imaginary part, respectively, where k x = ηk r . The analytical solution of X is thus assumed to be the superposition of plane waves propagating downstream and upstream. Then, the solution of pressure perturbationp can be expressed aŝ where X ± = ρk x,1 and the subscript '1' represents variables at the inlet. Coefficients C 1 and C 2 can be calculated by the given axial boundary conditions and the constant λ is determined by the mode numbers n θ and n r .

Solutions of the Three-Dimensional Velocity Perturbations
With the proposed solution ofp in the separated form, the radial and azimuthal momentum LEEs in Equations (13) and (14) can be transformed into non-homogeneous differential equations forû r andû θ . Consequently, corresponding analytical solutions can be derived directly via the variation of constants method, expressed as follows: where Coefficients C 3 and C 4 are determined by the relevant boundary conditions or the initial values of the radial and azimuthal velocity perturbations. For example, C 3 and C 4 equal zero when there are no radial and azimuthal velocity perturbations at the inlet (û r (x 1 ) =û θ (x 1 ) = 0). Furthermore, a zero axial component of the vorticity perturbation (ξ x = 0) in the incoming flow naturally provides the condition of C 3 = C 4 (this is discussed in Section 6).
Compared to sound waves travelling at the speed of sound, radial and azimuthal velocity disturbances involve waves propagating downstream at the axial mean flow velocityū x . Typically, these convection waves are associated with the development of the vorticity wave during the propagation of acoustic waves in the inhomogeneous background field [22,45,46].
Subtracting Equation (16) from Equation (17) yields a formula that does not contain the first derivative ofû x : Then, the axial velocity perturbationû x is easily obtained by substituting the solutions of the pressure and transverse velocity components into Equation (35). It should be noted that the r-derivatives of the Bessel function J n θ (λr) introduced byû r andû θ can be eliminated according to the Bessel equation of the radial component in the Equation (21). Therefore, the explicit expression is as follows: where and the terms with order higher than M 2 x are always neglected in the solutions. It is obvious that the axial velocity perturbation involves both the acoustic and convective waves at the same time. The latter greatly affects the resonant frequencies and growth rate of threedimensional combustors that have a velocity-dependence end, for example, an acoustically closed outlet. There are a total of four constant coefficients from C 1 to C 4 in the solutions of pressure and velocity perturbations that are determined by boundary conditions or initial values.
Finally, Equations (32), (33), and (36) constitute our analytical solutions of the threedimensional acoustic field, which should satisfy the high-frequency assumption |k 0 | |α| (the local wave number is sufficiently larger than the critical value) and that Mach number terms with order higher than M 2 x can be neglected. It should be noted that these proposed analytical solutions are not limited to the distribution forms of continuous axial mean temperature profiles.

Case of Zero Mean Flow
An assumption of zero mean flow is typically applied to the combustion chambers when Mach numbers are sufficiently small (M x ≈ 0). Therefore, the case of no mean flow is treated separately in this subsection and represents a benchmark result to present how the axial mean temperature gradient affects the thermoacoustic properties of the cylindrical combustion chamber containing the mean flow.
Substituting M x = 0 into the Equation (26) leads to the wave equation for the case of no mean flow: This coincides with the axial wave equation directly obtained from the LEEs for no mean flow (Equation (16) in our previous work [54]) if the imaginary part, as stated in the earlier text, is assumed to be sufficiently small, i.e., |k r | |M x k i |. The modified WKB method is used again to deal with spatially varying coefficients of the ODE, and solutions of b and a are derived by assuming |k 0 | |α|, expressed as follows: where the subscript '0' denotes the condition of no mean flow. The solution of the pressure perturbation is subsequently provided by where X ± 0 (x, ω) = ρk x,1 Then, the axial velocity is derived from the momentum LEEs in the condition of zero mean flow, expressed aŝ The radial and circumferential velocity perturbations have degenerated expressions, provided byû andû The analytical expressions for the three-dimensional acoustic field are finally derived in the case of zero Mach number. The two constant coefficients in the proposed solutions, C 1 and C 2 , can be determined by acoustic boundary conditions.

Validation Configuration
The proposed solutions are applied to predict the three-dimensional acoustic field for a straight cylindrical duct with the mean temperature profileT(x) and the inlet flow Mach number M x,1 . A linear temperature distribution is accounted for, with the expression where the axial mean temperature decreases fromT 1 at the inlet, x = x 1 = 0, toT 2 at x = l.
In order to validate the analytical solutions, the following boundary conditions are chosen for the inlet and outlet. An external pressure perturbationp in (x = 0, r, θ) = p 1 J n θ (λr) exp(in θ θ) is prescribed at the duct inlet and a pressure release conditionp out (x = l, r, θ) = 0 is prescribed at the outlet.
A dimensionless frequency is defined by Ω = f l/c 1 . The frequency larger than cut-off frequency, ω > ω c , is re-written as follows: Certain dimensionless transfer functions fromp 1 to acoustic perturbations are defined as functions of geometric positions and the forcing frequency, It is important that the analytical model has the almost same linearised Euler equations with the exception of the negligible entropic effect on the high-frequency acoustic waves. Hence, the proposed analytical solutions can be well verified by comparing them to the numerical results from LEEs. Table 1 presents the parameters used in the analyses in the next subsection.

Reconstructions of Sound Responses
In this subsection, the proposed analytical solution is validated by comparing the predicted three-dimensional acoustic field to the results of the numerical LEEs.
For the first transverse (1T) mode, n θ = 1 and n r = 0, dimensionless transfer functions are calculated analytically on the plane θ = π/2 when Ω = 0.8 and M x,1 = 0.2. As shown in Figure 2, our proposed solutions accurately reconstruct the three-dimensional acoustic field compared to the numerical results from LEEs. Figure 3 presents the real and imaginary parts of F p and F u x along an axial line of (r, θ) = (r 0 /2, π/2) for the 1T mode. The analytical solutions agree well with the numerical results, with the reduced forcing frequency Ω being 0.8, 1.2, and 1.6, respectively.  It is worth noting that the high-frequency condition of the modified WKB method, that is, |k 0 | |α|, can be re-written in the form of a dimensionless frequency: The cut-off frequency Ω c for a transverse mode (λ = 0) often provides a much larger critical value than Ω 0 .
Therefore, the high-frequency assumption has almost no effect on the accuracy of the proposed analytical solutions when applied to the transverse acoustic field. High precision can be achieved when the forced frequency of the perturbation is slightly different from the cut-off frequency.

Boundary Conditions and Eigenvalue Matrix
In this subsection, the proposed three-dimensional pressure and velocity solutions are applied to a cylindrical duct containing an axial mean temperature profile and a mean flow, then used to predict transverse modes. For the incoming flow (x = 0) at the inlet, no vorticity perturbation is assumed upstream of the varying temperature region.
The flow vorticity is derived by taking the curl of the three-dimensional velocity vector u, expressed as follows: The components for the vorticity vector can be calculated by substituting the proposed velocity solutions. Then, the zero vorticity perturbation at the inlet can be expressed in the form of four coefficients C 1 ∼ C 4 , written aŝ where ∆ = ik 0,1 − M x,1 α 1 . These zero vorticity conditions can provide two equations, as the radial and circumferential components have the same expressions when the axial vorticity equals zero. Another two equations can be provided by acoustic boundary conditions at the inlet and outlet, e.g., the open end (p = 0) and acoustically closed end (û x = 0) boundary conditions.
The eigenvalue system is consequently built by combining four boundary conditions of the duct ends, expressed as The complex angular frequency ω = 2π f + iω i can be solved by letting the determinant of the eigenvalue matrix M be zero, where ω i < 0 denotes that the corresponding thermoacoustic mode is unstable.

Effect of Mean Flow and Temperature Gradient on Resonant Frequencies and Instabilities
Predictions for the first transverse mode (1T) are carried out for the eigenvalue system in the presence of different outlet mean temperaturesT 2 . Two kinds of boundary conditions are prescribed at the inlet and outlet, namely, both open end boundary conditions or both acoustically closed ends. The latter is representative of the inlets and outlets of many real combustion chambers.
As shown in Figure 4, the frequencies and growth rates of the thermoacoustic modes vary with the outlet mean temperatureT 2 . Good agreements are found between analytical solutions and numerical LEE results for both kinds of boundary conditions. When the outlet mean temperature increases, the spatially averaged speed of sound increases, which results in a higher resonant frequency. It is clear that when the duct outlet temperatureT 2 differs substantially from the inlet temperatureT 1 , the growth rate does not equal zero. A decrease in axial mean temperature leads to an unstable mode, while temperature increase corresponds to a stable mode. The resonant frequencies for the open-open ends are almost same as those of the closed-closed ends. However, the growth rates for the closed-closed ends have larger absolute values irrespective of whether the thermoacoustic modes are stable or unstable. Figure 5 presents the evolution of the resonant frequency and growth rate as a function of the inlet flow Mach number M x,1 . The analytical results calculated by the proposed solutions present good consistency with the results of the numerical LEEs. The growth rate is zero when the mean flow is stationary. With a baseline of no mean flow, the 1T mode becomes more unstable as the incoming flow Mach number increases, especially for the closed-closed ends.  These results show that the system with both mean temperature profile and mean flow presents different linear stabilities from that with either uniform mean temperature or no mean flow. The presence of mean temperature gradient and the mean flow lead to evident shifts of the resonant frequencies and growth rates. Therefore, both should be taken into account during predictions and simulations for sound propagation through a three-dimensional cylindrical duct.

Transverse Modes with Different Boundary Conditions
In addition to the 1T mode, the first radial mode and the second transverse mode were investigated in order to provide general verifications of the analytical solutions. With both kinds of boundary conditions, as presented in Table 2, analytical predictions of resonant frequencies and growth rates are in good agreement with the numerical LEE results in the case of M x,1 = 0.1 andT 2 = 1200 K. High precision can be achieved by applying the proposed WKB-type solutions in the eigenvalue system. The convection wave generated by the acoustic waves propagating through the region of varying mean temperature and inhomogeneous mean flow results in a larger absolute value of the growth rate when a velocity-dependent boundary is prescribed.

Conclusions
In this work, we obtain an analytical solution for a three-dimensional acoustic field in a cylindrical duct in the presence of an axial temperature gradient and mean flow. This paper extends our previous works involving a one-dimensional straight duct [43] to a three-dimensional cylindrical duct. We first account for the multi-dimensional acoustic perturbations that propagate through the inhomogeneous field; both the axial mean temperature distribution and mean flow are present. A second-order partial differential wave equation is obtained from linearised Euler equations (LEEs) in the cylindrical system. Then, the separation of variables and a modified WKB method are applied to derive the analytical expressions of pressure and velocity perturbations when the high-frequency assumption |k 0 | |α| is satisfied. Based on the proposed analytical solutions, an eigenvalue system is built by providing boundary conditions of both duct ends and assuming zero vorticity perturbation at the inlet.
Validation is conducted by comparing the analytical solutions to the results of numerical LEEs. For a linear axial temperature profile, the proposed analytical solutions are applied to predict the acoustic perturbations in the cylindrical geometry for three transverse modes and two kinds of boundary conditions. The results show that the proposed analytical solutions can provide accurate predictions for the three-dimensional acoustic field and the thermoacoustic modes with both the temperature gradient and the mean flow considered. High precision is achieved when frequencies are larger than the cut-off frequency. Both axially varying mean temperature and mean flow can change the apparent growth rate within a thermoacoustic system. Therefore, the assumption of either no mean flow or zero temperature gradient should be used cautiously in the low order modeling of combustion systems.