NOMENCLATURE
-
$c$
, -
sound speed,
${\rm{m}}/{\rm{s}}$
- E,
-
total energy, J
-
${E_a}$
,
-
acoustic energy, J
- L,
-
length, m
- P,
-
pressure, Pa
-
$\bar{p}$
,
-
average pressure, Pa
-
$p^{\prime}$
,
-
pressure perturbation, Pa
-
$Q$
,
-
power,
${\rm{W}}/{{\rm{m}}^3}$
-
$\bar{Q}$
,
-
average power,
${\rm{W}}/{{\rm{m}}^3}$
-
$Q^{\prime}$
,
-
power perturbation,
${\rm{W}}/{{\rm{m}}^3}$
- R,
-
gas constant,
- t,
-
time, s
- T,
-
temperature, K
-
$\vec u$
,
-
velocity,
${\rm{m}}/{\rm{s}}$
-
${\vec u^{\prime}}$
,
-
velocity perturbation,
${\rm{m}}/{\rm{s}}$
- V,
-
volume,
${{\rm{m}}^3}$
- x,
-
coordinate, m
Greek Letters
- α,
-
dissipation coefficient
- β,
-
source parameter
-
${{\gamma }}$
,
-
heat capacity ratio
-
${\rm{k}}$
,
-
thermal conductivity, W/m-K
-
${{\mu }}$
,
-
dynamic viscosity,
-
${{\rho }}$
,
-
density,
${\rm{kg}}/{{\rm{m}}^3}$
-
$\bar{\tau}$
,
-
stress tensor, kg/m-s2
- Θ,
-
azimuthal coordinate
1.0 INTRODUCTION
Development of future hypersonic flight vehicles will require advancements in ground test facility technology, including the capability to conduct aeropropulsion testing with air that properly simulates the stagnation conditions associated with high-speed flight through the Earth’s atmosphere(Reference Fetterhoff, Kraft and Laster1–Reference Saari, Jauch and Chu4). In general, liquid rocket engine technology is applied to meet these requirements by burning little fuel and heating up a large amount of air. High-frequency combustion instability is considered as a major challenge. When it occurs, the pressure fluctuates, which can introduce large thermal and mechanical stresses in the combustor, resulting in decreased performance, unacceptable vibrations or even engine failure(Reference Harrje and Readon5–Reference Culick7).
The fundamental mechanisms leading to high-frequency combustion instability in air heaters are the coupling between the sub-processes and the acoustic oscillation in combustor, such as injection, atomisation, vaporisation, mixing and chemical reactions. Pressure oscillations amplify and combustion instabilities self-sustain when heat is released near a region of maximum pressure. This is known as Rayleigh criterion. Only high-fidelity models can capture the above complex coupling process. Reliable prediction of high-frequency combustion instability is still a challenging problem. Modelling approaches used for combustion instability analysis have ranged from acoustic wave equations to detailed computational fluid dynamics. Harry et al.(Reference Harrje and Readon5) and Yang et al.(Reference Yang and Anderson6), for example, summarised the state of art in the 1990s. A more recent review on modelling is given by Pieringer et al.(Reference Pieringer, Sattelmayer and Fassl8) and Culick(Reference Culick7). Such models have contributed immensely to the understanding of high-frequency combustion instability.
Wave-equation approaches can provide a great deal of insight into the coupling of combustion heat release with the acoustic modes. With these methods, characteristics of oscillations can be investigated, such as resonant frequencies, mode shapes and pressure growth rates. On the other hand, these methods have the advantage of low numerical cost, and the investigations can be carried out on a personal computer with a single processor. Crocco(Reference Harrje and Readon5), Zinn(Reference Zinn9), Culick(Reference Culick10) and Yang et al.(Reference Yang and Anderson6) modelled some pioneering studies with these approaches. Recently, researchers from Purdue University conducted computational investigation of acoustics and instabilities in a longitudinal-mode rocket combustor with wave-equation approaches(Reference Portillo, Sisco, Yu and Anderson11–Reference Sisco, Yu, Sankaran and Anderson14). Their results show that the resonant frequencies and the unsteady modes of their predictions agree well with the experimental data. However, in the derivation of the wave equation, effects of body force and viscosity have been neglected, and it has also been assumed that the pressure disturbance is infinitesimally small and the process is isentropic. In general, pressure oscillations do not exceed 30% of the mean pressure(Reference Quinlan, Kirkpatrick, Milano and Mitchell15). In addition, the combustion response function is required to represent the flame response to incident fluctuation. Because of the inherent simplifications and assumptions mentioned above, the wave equation cannot simulate large-amplitude pressure oscillation, entropy fluctuations, mean flow and viscosity or dissipative flow near the wall surfaces. Moreover, wave-equation models use classic acoustic boundary conditions such as open- or closed-end tubes or admittance functions that are different to the actual physical boundary condition. To match experimental observation, these methods need several parameters, for instance, the position of mean heat addition, the position of unsteady heat addition and the time lag. Depending on which mode attention is focused, the above input parameters may be changed(Reference Frezzotti, Terracciano, Nasuti, Hester and Anderson16). Air heaters are characterised by several aspects, such as entropy waves, mean flow, viscosity and non-linear pressure oscillations. In addition, a converging-diverging nozzle is used at the end of the chamber, so the idealised boundary cannot be used. Therefore, it is unreasonable to investigate acoustics and instabilities in air heaters with wave-equation methods.
The highest-fidelity modelling is the Large Eddy Simulations (LES) of the multi-phase turbulent reacting flow field in the combustor(Reference Harvazinski, Talley and Sankaran17–Reference Yuan and Shen22). These methods include effects of entropy waves, non-linear oscillations, viscosity and physically realistic boundary conditions. Moreover, these approaches have the advantage that they do not necessarily require response functions to represent the dynamic response. However, given the limitations of present-day computers or even the most sophisticated supercomputers, they are practical only for the simple combustor with few injectors. Air heaters are characterised by an injector–combustor–nozzle complex configured with several hundreds of shear-coaxial injector elements. Thus, it is expensive and time consuming to predict combustion instability using LES. Also, it is very complicated, using an intermediate level of modelling, which solves the unsteady Reynolds-averaged Navier–Stokes equations, coupled with appropriate response functions.
In general, air heaters are equipped with so many injectors that fuel is uniformly injected into the main flow. On the other hand, the main purpose of air heaters is to burn little fuel (no more than 10% of the total mass flow rate) and heat up a large amount of air. Based on the characteristics of the air heater and firing tests, a new computational model was developed in this paper, which involved solving the unsteady Reynolds-averaged Navier–Stokes equations, coupled with special source terms. With these new methods, acoustics and instabilities in air heaters were investigated numerically in this paper. Compared with wave-equation methods, these new approaches considered the effects of viscosity, entropy waves and non-linearity. Also, the actual physical boundary condition was included automatically(Reference Yu, Sisco and Anderson13,Reference Sattelmayer, Kathan and KÖglmeier23–Reference Qin, Zhang and Wang26) . When compared with LES, it was practical for air heaters with hundreds of injectors, complex configurations and geometries. The numerical investigations could be carried out on a single computer with eight cores. To the authors’ knowledge, it was the first time to try to predict and analyse combustion instability in air heaters using these novel ideas and methods. For this reason, the objective of these studies was to show the general feasibility of the method and to explain the general behaviour of the model. In a conventional analysis, combustion is assumed to mainly take place in the so-called rapid combustion zone, and this zone is usually modelled as a disc in the combustor near the injection head. However, in practice, the flame has a spatial distribution. A model that treats the flame as a disc is not close to physical reality. For instance, Kato et al.(Reference Kato, Fujimori, Dowling and Kobayashi27) developed a one-dimensional linear model and demonstrated that stability boundaries can be predicted with higher accuracy when spatial distribution of heat release accounted for. Thus, an attempt was made in this article to describe the effect of flame spatial distribution on predictions of oscillation frequency and mode.
The article was organised as follows. First, a brief description of the experimental configuration and sample results were presented. This was followed by a description of the process of modelling and verification studies using analytical solutions of acoustics in a closed rectangular cavity. We then presented the computational acoustics and instabilities for the experimental air heater. In the final section, we conclude with summaries and discussions of future research work.
2.0 AIR HEATER CONFIGURATION AND FIRING TEST RESULTS
The detailed geometry of the present air heater is shown in Fig. 1. It consists of a cylindrical chamber with a chocked nozzle and injectors installed on the faceplate. The chamber consists of two cylindrical sections with different diameters and lengths. The injector faceplate consists of hundreds of shear-coaxial injector elements. All injector elements have no recess and same dimensions. The operation conditions are summarized in Table 1. To measure the high-frequency pressure oscillations, four pressure transducers were used and the signals were sampled at 5kHz. As shown in Fig. 1, one high-frequency measurement was located at the middle of the combustor. The other three pressure transducers were installed in the air, liquid oxygen and ethanol manifolds respectively.
Operation conditions
Schematic of the air heater.
Pressure time history for a typical test.
Nineteen firing tests were conducted in all. Sixteen of the firing tests showed unstable combustion characteristics, that is, the peak-to-peak fluctuation amplitudes of chamber pressure at a large oscillation level were more than 10% of mean value. Two of the firing tests were forced to be shut down because of the large oscillation level. One of the firing tests showed an intermediate pressure oscillation with peak-to-peak fluctuation amplitudes of about 10% of mean value. A typical pressure time history is given in Fig. 2. As seen from this figure, the chamber pressure showed large peak-to-peak fluctuation amplitudes of about 34.7% of mean value. For the present hot firing test, air, liquid oxygen and ethanol manifold pressure showed stable characteristics, which indicated intrinsic high-frequency combustion instability. Intrinsic instability is dependent on the combustor acoustics and independent of the supply system. Thus, it is expected to analyse the combustor acoustics and instabilities in the present air heater. Figure 3 depicts a waterfall plot of chamber pressure for the analysis window from 23.5s to 24.7s. As shown in this figure, there existed four dominant frequencies after ignition. Because of the experimental facility limitation, pressure transducer sampled at 5kHz, the oscillation below 2.5kHz could be distinguished. Thus Fig. 3 could only give a Fast Fourier Transform (FFT) spectrum below 2.5kHz, over the time interval from 23.5s to 24.7s. As seen from the FFT result, multiple unstable peaks were observed in the test. It showed strong peaks around 427.2Hz, 1,112Hz, 1,322Hz and 1,570Hz. Acoustics and instabilities in the air heater will be investigated with a new computational model in the following sections.
Waterfall plot of chamber pressure over the time interval from 23.5s to 24.7s.
3.0 MODELLING STRATEGY
3.1 Simplifications of physical processes
As previously stated, the ratio is just 6% of fuel mass flow rate to total mass flow rate. Fuel is injected uniformly into the main flow through shear-coaxial injector elements installed on the faceplate. According to previous results of spray combustion, the combustion heat release mainly takes place within a short distance of the front of the chamber. The main physical processes inside the air heater are to heat up a large amount of air by consuming little fuel. Thus, the main physical processes can be assumed as follows (Fig. 4). Firstly, air is injected uniformly from the faceplate with room temperature (300K) and constant mass flow rate. Then, it is heated up to the expected temperature when it flows through a heat release zone, i.e. combustion is represented as a distributed heat release source in space. Finally, the heated air is choked through a convergent-divergent nozzle. These assumptions are acceptable simplifications of the physical realities due to the aspects mentioned above.
Simplified model of air heater.
3.2 Model of air heater
Based on the above assumptions and significations, the control equations of physical processes in the present air heater can be written as follows:
conservation of mass,
\begin{equation}\frac{{\partial \rho }}{{\partial t}} + \nabla \cdot (\rho \vec u)=0\end{equation}
conservation of momentum,
\begin{equation}\frac{{\partial (\rho \vec u)}}{{\partial t}} + \vec u \cdot \nabla (\rho \vec u) = - \nabla p + \mu {\nabla ^2}\vec u + \frac{\mu }{3}\nabla \left( {\nabla \cdot \vec u} \right)\end{equation}
conservation of energy,
\begin{equation}\frac{{\partial (\rho \vec u)}}{{\partial t}} + \vec u \cdot \nabla (\rho E) = - \nabla \cdot (p\vec u) + \nabla \cdot \left( {k\nabla T} \right) + \nabla \cdot (\bar \bar \tau \cdot \vec u) + Q\end{equation}
state equation,
\begin{equation}p = \rho RT\end{equation}
3.3 Modelling heat release
The heat release rate consists of a mean part and a fluctuating one, as given in Equation (5). According to chemical reaction (Equation (6)), the total power can be taken as
$2.39 \times {10^8}\;W$
, if the ethanol is completely consumed. We assumed that the heat was released uniformly within a straight section. As shown in Fig. 5, we investigated instabilities under different heat release zones. In this context, all dimensional quantities of air heater were non-dimensionalised using the radius of nozzle throat. Volume-weighted average power for different heat release zones is summarized in Table 2.
\begin{equation}Q = \bar Q + Q^{\prime}\end{equation}
Chemical reaction:
\begin{equation}{C_2}{H_5}OH+3{O_2} \to 2C{O_2}+3{H_2}O\end{equation}
Volume-weighted average power for different heat release zones
Schematic of heat release zone for different cases.
In present calculations, the effects of turbulence were neglected and the spray combustion processes were simplified, so they did not include the mechanism for coupling the heat release to the acoustic modes; the instability could not be incited and sustained without this coupling mechanism. In general, such a coupling may be introduced in the form of response functions. However, we cannot derive a general theory to design such response functions; the models would need to be redesigned and calibrated to match the particular test conditions. In this article, a simple and preliminary model for combustion response was introduced to model the coupling of heat release to the acoustic modes(Reference Smith, Xia, Anderson and Merkle28,Reference Smith, Ellis, Xia, Sankaran, Anderson and Merkle29) (Equation (7)), which was expressed as a function of the instantaneous pressure with a time tag. In the present model,
${\rm{\xi }}$
is a scaling function that is dependent on space, which is used to introduce the unsteady heat release continuously in space. In our calculations,
${\rm{\xi }}$
is set to be non-zero within the heat release zone and zero everywhere else. The following parameters were used: the normalised interaction index between heat release and acoustic modes, n = 0.2, and the phase lag,
$\tau =0$
. As seen from Equation (7), the relationship between fluctuating heat release and pressure oscillations shows positive feedback; thus, both oscillations grow unbounded until balanced by the damping of the system. Hence, we artificially controlled the oscillations by stipulating that the maximum
$Q^{\prime}$
at any location did not exceed
$b\bar{Q}$
(Equation (8)). As stated by Hantschk(Reference Hantschk and Vortmeyer30), the limit cycle amplitude is mainly determined by non-linearities in the heat flux from the heat source to the fluid. In our calculations, variation of the maximum
$Q^{\prime}$
was investigated to match the pressure limit cycle amplitude of firing tests. We could get a good agreement for b = 0.2; thus, in our model, b was set to be 0.2.
\begin{equation}Q^{\prime} = \bar Q \cdot n \cdot \xi (x) \cdot \frac{{p^{\prime}(x,t) - \bar p(x,t - \tau )}}{{\bar p(x,t - \tau )}}\end{equation}
\begin{equation}\max \left( {Q^{\prime}} \right) = b\bar Q\end{equation}
3.4 Analysis of heat-to-sound energy conversion in the model
By linearising the control equations (Equations (1), (2), (3)) and neglecting the viscous losses and the conduction effect we got in Equation (9), the derivation was provided in Refs.(Reference Culick7,Reference Poinsot and Veynante31,Reference Zhao, Li, Yang and Zhang32) . Acoustic energy
${E_a}$
is defined as Equation (10). Equation (9) illustrates how sound is produced from unsteady heat release or the process of heat-to-sound conversation. When the heat release fluctuations are in phase with the pressure oscillations, acoustical energy is increased. In our calculations, the time lag has been set to zero. This causes the unsteady heat release to be exactly in phase with the pressure oscillations, which corresponds to enforcing Rayleigh’s criterion.
\begin{equation}\frac{{\partial {E_a}}}{{\partial t}} = \frac{{\gamma - 1}}{{\gamma \bar p}}\int_V {p^{\prime}Q^{\prime}dv}\end{equation}
\begin{equation}{E_a} = \int_V {\left( {\frac{1}{2}\rho {{\vec u{\prime}}^2} + \frac{1}{2}\frac{{{{p{\prime}}^2}}}{{\gamma \bar p}}} \right)} dv\end{equation}
3.5 Mesh and boundary conditions
Figure 6 depicts the computational mesh, and Fig. 4 shows boundary conditions employed in the present simulation. The domain was meshed with 324,576 structured hexahedral elements, 982,116 faces and 333,185 nodes. Characteristic inflow with constant mass flow rate and temperature was employed at inlet, and it was found that the acoustically closed-end condition was nearly realised(Reference Pieringer, Sattelmayer and Fassl8,Reference Shimizu, Morii and Hori24) . Supersonic outflow boundary condition was employed for the choked nozzle, i.e. an extrapolated boundary condition was used at the supersonic exit, and the parameters were not imposed but followed from the solution. No-slip and adiabatic conditions were applied at all wall surfaces.
Mesh for the current calculation.
3.6 Numerical verification
A second-order implicit formulation in time along with second-order accurate spatial discretisation was used in the solution process. The governing equations of continuity, momentum and energy were solved simultaneously. To verify the numerical methods, here, we presented the verification of the numerical predictions for acoustics in a closed rectangular cavity by comparing them with the analytical solutions. The analytical solutions have been used previously for acoustic and instability validation(Reference Smith33–Reference Richecoeur35). As shown in Fig. 7, the closed cavity is 1m long and 1cm wide, and the origin of coordinates is located in the middle of the cavity. The pressure probe was located in the entrance of the cavity. Waves propagate along the x direction of the cavity, and they reflect back and forth between the left and right side-walls to form a standing longitudinal mode.
Schematic of the cavity.
The cavity is initially filled with quiescent air. The operating conditions for these cases are presented in Table 3. The origin of coordinates is located in the middle of the cavity. At
$t=0$
, the first longitudinal mode in the cavity was obtained by applying an initial axial velocity component,
$U=0.5{\rm{sin}}\left[ {\left( {x+0.5} \right)\cdot\pi } \right]$
. The pressure evolutions at the entrance of the cavity obtained by analysis are
\begin{equation}p=101325 - \rho Uc \cdot \cos \left( {\pi \cdot \left( {x + L/2} \right)} \right) \cdot \sin \left( {c \cdot \pi \cdot t} \right)\end{equation}
The period of the first longitudinal mode is
$T=2L/c=5.76{\rm{ms}}$
, and the corresponding frequency is
$f=173.5{\rm{Hz}}$
.
The simulation was conducted for
$\sim$
250ms, after which the pressure and velocity evolutions were obtained. Figure 8 depicts a comparison of the computed and analysed pressure time histories at a point located in the cavity entrance. As seen from the figure of pressure time history, the simulated result shows good agreement with analytical data. Spectra of the pressure signals presented in Fig. 9 indicate that 170.9Hz is the dominant frequency, which is approximately equal to the analytical value (173.5Hz). Moreover, the simulation generated the harmonics of the fundamental frequency, such as 341.8Hz, 512.7Hz and 680.6 Hz. Thus, the present calculations also captured the non-linear effects.
Operating condition of the closed rectangular cavity
Time history of pressure oscillation at the entrance of the cavity.
Power spectra of the simulated pressure signals.
Figure 10 shows the simulated pressure evolutions during a period corresponding to the first longitudinal oscillation. The dots on the solid line in Fig. 8 correspond to the times of each snapshot. Figure 11 gives the distribution of non-dimensional wall pressure plotted over an entire cycle. Both Figs 10 and 11 show that the dominant mode of pressure fluctuation is the first longitudinal mode.
Simulated result of the first longitudinal oscillation. The dots on the solid line in Fig. 8 correspond to the times of each snapshot.
Distribution of non-dimensional wall pressure plotted over an entire cycle.
4.0 RESULTS AND DISCUSSION
Combustion instability is caused by the interaction between the acoustics in the combustor with the multi-phase turbulent reacting flow dynamics. Significant computational resources are required to model these effects, and it is beyond the scope of the present paper. Consequently, the present simulations were accomplished by two procedures. Firstly, the computational model was operated in steady-state mode to obtain the mean-flow solution that was relevant for a given test condition. Secondly, based on the solutions from steady-state mode, the code was run in unsteady mode to obtain the acoustic response. The second procedure was accomplished by introducing a combustion response function in the energy equation. The convergence of grid and time step was confirmed because the same acoustic wave form was obtained using simulations with different grids and different time steps (see Appendix for details).
4.1 Mean flow
To determine the flow field inside the air heater obtained by the present methods, in this section the flow field without the combustion response function is shown (
$Q^{\prime}=0$
). Figure 12 gives the temperature contours in the X–Y plane for different cases. The cold air is heated gradually when it flows through the heat release zone, and then the temperature decreases gradually when the flow is accelerated by the converging-diverging nozzle.
Temperature contours in X–Y plane for different cases.
Temperature profiles obtained for the steady state solution are plotted in Fig. 13. Together with Fig. 12, it is clear that the shorter the heat release zone, the steeper the temperature profile. However, the temperature exhibits a same value at the end of the heat release zone for all cases, which matches the design conditions.
Table 4 summarises the comparison of designed values and numerical results. Except for temperature, computational results show good agreement with the designed values. In our calculations, the effects of multi-species and thermal pyrolysis were not considered, so the computational chamber temperature was over-predicted by about 207K compared with the value determined by the thermodynamic calculation program (1558K). Because of the over-predicted chamber temperature, the chamber pressure was about 4% higher than the designed value (5MPa).
Comparison of designed values and numerical results
Geometry profile and steady-state solution obtained with current model for different heat release zones.
4.2 Unsteady results
The computational results presented so far have involved mean flow in the chamber. In this section, we present results using the combustion response function described in Equation (7). In all the calculations shown here, the unsteady calculations were first initiated by using broadband excitation, and this was done by patching a pressure peak in the form of a spatially spherically decaying pressure pulse with a Gaussian profile onto the initial steady solution (shown in section 4.1). This spherical Gaussian pressure is defined by the position of its maximum
${\vec x_c}$
, the pressure amplitude at the maximum
$p_{max}^{\prime}$
, and its half-value radius
$b$
, and it is given by Equation (12). A frequency spectrum with a Gaussian frequency distribution is introduced by this Gaussian pressure distribution(Reference Pieringer, Sattelmayer and Fassl8). If the larger half-value radius of the pressure distribution is chosen, the fewer frequencies are excited. Specifically, a Dirac peak would contain an infinite number of frequencies. In the present study, the extension of the spherical Gaussian pressure was chosen so that the first five eigenfrequencies could be excited. In the present calculations, the following parameters were used:
$p_{max}^{\prime}=20\% {\rm{p}}$
;
${\vec x_c} = \left( { - 0.82,\;\;0.08,\;\;0} \right)$
;
$b=0.09$
. Figure 14 shows the initial pressure field for the unsteady calculations and also demonstrates the location of the initial spherical Gaussian pressure pulse. Once the above patching process was completed, the combustion response function was turned on. Then, the limit cycle was formed gradually for all cases, as the simulation time passed.
\begin{equation}p^{\prime}(\vec x,0) = {p^{\prime}_{\max }} \cdot \exp \left[ { - \ln 2\frac{{{{\left( {\vec x - {{\vec x}_c}} \right)}^2}}}{{{b^2}}}} \right]\end{equation}
Location of the Gaussian pulse in the computational domain.
4.2.1 Mode
To analyse the simulation results conveniently, we set three planes in the domain perpendicular to the chamber axis shown in Fig. 15, with each plane having four probes. Probes of the 1st to 4th were located counter-clockwise on plane 1, probes of the 5th to 8th were located counter-clockwise on plane 2 and probes of the 9th to 12th were located counter-clockwise on plane 3.
Schematic of the analysis planes and the probes.
Section 4.1 has demonstrated the mean flow field. Here, mode shapes are investigated. After the limit cycle has formed, the waveforms at probe 1, probe 2, probe 3 and probe 4 were similar, and they exhibited a phase difference for the adjacent probes. To show details, plots of unfiltered pressure time histories taken at probes 1, 2, 3 and 4 are shown in Fig. 16, which gives the waveforms of limit cycle. The travelling compression wave can be seen by examining the peaks in the pressure data. Combining the positions of the probes, it is expected that the spinning first tangential (1T) mode applies to cases (a), (b) and (c). Note that the spinning direction may be arbitrary. For case (d), as shown in Fig. 16(d), the signals taken at probe 2 and probe 4 are out of phase with signals taken at probe 1, the reference signals. Likewise, the signals taken at probe 3 are in phase with signals of probe 1. This indicates that the dominant instability has the features of a standing second tangential (2T) mode for case (d). These features will be highlighted later.
Pressure time histories through probe 1 to probe 4.
Time snapshots of the entire combustion flow field at specific instants of time are useful for deducing the global character of the entire flow field. When these time snapshots are spaced sufficiently closely together in time, they can be converted to flow movies, which are helpful for understanding dynamics. Figure 17 shows the temporal behaviour of the pressure field over an entire cycle for different cases. As discussed earlier, 1T oscillations are observed for cases (a), (b) and (c), and their modes are found to be spinning, while a standing 2T wave is clearly observed for case (d). It is worth noting that azimuthal modes can appear as standing or spinning modes, and both are observed in gas turbines. Bifurcations between standing and turning modes were explained previously. Schuermans et al.(Reference Schuermans, Paschereit and Monkiewitz36–Reference Wolf, Staffelbach, Gicquel, MÜller and Poinsot38) attribute these bifurcations to non-linear effects. They propose a non-linear theoretical approach showing that standing wave modes can be found at low oscillation amplitudes, but only one spinning mode is found for large-amplitude limit cycles. Similarly, we simulated many different cases having unique heat release zones but different maximum
$Q^{\prime}$
(Equation (8)). Unfortunately, these bifurcations cannot be observed in our calculations. Some researchers explained these bifurcations using linear approaches: spinning modes would appear only in perfectly axisymmetric configurations, while any symmetry modification would lead to standing modes(Reference Wolf, Staffelbach, Gicquel, MÜller and Poinsot38). Results from Ref. (Reference Evesque, Polifke and Pankiewitz39) suggested that the initial flow conditions might trigger either standing or spinning modes. Some researchers observed that turbulence can cause random mode switching between standing and spinning structures(Reference Schuermans, Paschereit and Monkiewitz36,Reference Noiray, Bothien and Schuermans40) . According to present calculations, it was found that modulation of heat release length might cause mode switching between standing and spinning structures. Increment of heat release length tended to excite higher-order standing modes.
Temporal behaviour of the pressure field over an entire cycle on the slice (
$x = - 9.375$
). Arrowhead denotes the spinning direction of the pressure wave. Conditions correspond to Table 2.
Transverse oscillation modes have been discussed. In what follows, attention is focused on the longitudinal modes. A clear picture of the underlying longitudinal modes emerges in the present simulation, and the resulting wall pressure is plotted at different instants in time over an entire cycle. Such longitudinal mode shapes are shown in Fig. 18. Superimposed on Fig. 18 are the labels of the heat release zone. Examining this figure, it is apparent that the mode shape in the heat release zone appears to be approaching the half-wave mode for cases (a), (b) and (c). As the heat release zone is extended to
$x = - 0.625$
from the faceplate (case (d)), the first longitudinal mode becomes the dominant mode. For the rest of the chamber region, the case (a) mode is the third longitudinal mode and the case (b) mode is the second longitudinal mode, whereas the half-wave mode is clearly observed for case (c) and case (d). All cases considered in this context showed a complex wave shape at any given instant in time, and the pressure reached a peak at different nodal locations at different instances in time. As Smith(Reference Smith, Ellis, Xia, Sankaran, Anderson and Merkle29) explained, it is the combined effects of the forward-running
$u + c$
and the backward-running
$u - c$
waves that cause this inherent phase difference in the wave shape. In all, when the heat release zone locates near the faceplate where it is the pressure antinode of the longitudinal mode, increment of heat release length tends to excite lower-order longitudinal mode.
Computed mode shapes for the longitudinal mode. Conditions correspond to Table 2.
Overall, combined longitudinal-transverse acoustic modes occurred in the present air heater, illustrated in Fig. 19.
Pressure mode shapes for different heat release zones.
When the system reached the limit cycle, the combustion response function was turned off. Figure 20 shows the time histories of pressure oscillations for case (a) after its combustion response function was turned off. A full view of the decay process for pressure fluctuation is shown in the inset of Fig. 20. The fluctuations are observed to be damped, and the flow field gradually becomes stable. In the absence of combustion response function, it was found that the numerical and real viscosity was sufficient to damp out any oscillations within a short time. Accordingly, the unsteady results presented in this paper were predominately acoustic in nature.
Time histories of pressure oscillations for case (a). The combustion response function is turned off at 0.577s.
4.2.2 Frequency
We next turn our attention to the frequency response predictions of the present air heater. Figure 21 depicts a waterfall plot for the pressure signals at probe 1 after the limit cycle has developed, which shows the transition of the frequencies and amplitudes as time proceeds. The sampling duration is 80ms, and the frequency resolution is 12.5Hz. As anticipated, in Fig. 21 the first three harmonics can be clearly observed, and these frequencies do not change considerably during the period of the limit cycle. The higher peaks represent harmonics of the primary mode occurring at integer multiples of its frequency, but with a steadily decreasing amplitude.
Frequency comparisons using the classical acoustic model and the present computational model
Waterfall plot for the pressure signals at probe 1 for
$\sim$
80ms.
The resonant mode frequencies of the air heater may be estimated using Equation (13) based on the following assumptions: the chamber of the air heater was viewed as a closed cylindrical duct with no nozzles and injectors; inviscid uniformly distributed ideal gas was considered; the condition in the combustor was assumed to satisfy linear acoustics, namely small-amplitude oscillations; and there was low-speed mean flow. Table 5 summarises the frequency using the classical acoustic model and the present computational model. Note that the chamber consists of two cylindrical sections with different diameters and lengths. In addition, there exists a notable temperature gradient along the flow direction in the first domain of the combustor. Thus, only the conditions in the last domain of the combustor satisfy the above assumptions. The results obtained from the present model were extracted from Fig. 21. It is evident that the computed frequency agrees fairly well with the theoretical value.
\begin{equation}{f_{mnq}} = \frac{c}{2}\sqrt {{{\left( {\frac{{{\alpha _{mn}}}}{{{R_c}}}} \right)}^2} + {{\left( {\frac{q}{{{L_c}}}} \right)}^2}}\end{equation}
where
${\alpha _{mn}}$
is
${n_{th}}$
root of the
${J_n}$
Bessel function of the first kind, c is the speed of sound in the medium,
${L_c}$
is the length of cylinder,
${R_c}$
is the radius of the cylinder and m, n and q are the wavenumbers for the tangential, radial and longitudinal modes, respectively.
It is interesting to compare the oscillation characters found in experiment with those obtained from computation and theoretical analysis. As already noted, strong peaks around 427.2Hz, 1,112Hz, 1,322Hz and 1,570Hz are observed on hot firing test (Fig. 3). The theoretical frequency of the first longitudinal mode is about 448Hz, estimated by Equation (13). Interestingly, the first experimental dominant frequency (427.2Hz) matches this value. The second dominant frequency (1,112Hz) shows good match with the frequency found in case (c). The third dominant frequency (1,322Hz) matches well with the frequency found in case (b). Also, the match between the last dominant frequency (1,570Hz) and that found in case (a) is good. Note that, according to our calculations, every case has its own oscillation modes and the corresponding frequencies, and these modes and frequencies, were insensitive to the changes of interaction index between heat release and acoustic modes, and to the maximum
$Q^{\prime}$
(Equations (7) and (8)), but sensitive to the heat release length. It was logical to assume that the flame might exhibit a significant unsteady motion; thus, the heat release length varied during the experiment, which led to different oscillation modes and frequencies. A waterfall plot of chamber pressure may support this assumption (Fig. 3), the peak frequency exhibiting maximal amplitude shifts among the dominant frequencies as the experiment time progressed, indicating that the heat release length is changing. However, this assumption needs more tests. Overall, the present model successfully predicted the dominant frequencies observed on hot firing test, except the first pure longitudinal mode.
4.3 Theoretical analysis of the modal behaviour in air heater combustor, i.e. standing or spinning waves
Considering the geometry of the present air heater combustor, the wave equations can be written in polar coordinates:
\begin{equation}\frac{{{\partial ^2}p^{\prime}}}{{\partial {\theta ^2}}} - \alpha \frac{{\partial p^{\prime}}}{{\partial t}} - \frac{{{\partial ^2}p^{\prime}}}{{\partial {t^2}}} = - \frac{{\partial Q^{\prime}}}{{\partial t}}\end{equation}
The Galerkin method is used to solve Equation (14), assuming that the pressure field can be expressed as the sum of a series of orthogonal basis functions
$p^{\prime}\left( {{\rm{\theta }},{\rm{t}}} \right) = \mathop \sum \nolimits_1^\infty {\eta _n}\left( t \right){\psi _n}\left( \theta \right)$
. When considering the situation where the thermoacoustic coupling in the air heater results from the combination of the eigenmodes sharing the nth eigenvalue, the pressure field can be approximated by Equation (15).
\begin{equation}p^{\prime}(\theta ,t) = {\eta _1}(t) \cdot \cos (n \cdot \theta ) + {\eta _2}(t) \cdot \sin (n \cdot \theta )\end{equation}
It is assumed that the modal amplitudes can be expressed as Equation (16), in which A and B vary slowly in time.
\begin{equation}\left\{ {\begin{array}{*{20}{c}}{{\eta _1} = A(t)\sin (nt)}\\\\[-8pt] {{\eta _2} = B(t)\cos (nt)}\end{array}} \right.\end{equation}
Thus, the expression for pressure field
$p^{\prime}$
can be decomposed as Equation (17).
\begin{equation}\begin{array}{l}p(\theta ,t) = {\eta _1}(t) \cdot \cos (n \cdot \theta ) + {\eta _2}(t) \cdot \sin (n \cdot \theta )\\{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = A\left( t \right)\sin \left( {nt} \right) \cdot \cos (n \cdot \theta ) + B(t)\cos \left( {nt} \right) \cdot \sin (n \cdot \theta )\\{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \underbrace {B\left( t \right)\sin \left( {nt + n\theta } \right)}_{Spinning{\kern 1pt} {\kern 1pt} {\kern 1pt} mode} + \underbrace {\left[ {A\left( t \right) - B(t)} \right]\sin \left( {nt} \right)\cos \left( {n\theta } \right)}_{Standing{\kern 1pt} {\kern 1pt} mode}\end{array}\end{equation}
From Equation (17), it can be concluded that, when the modal amplitudes A and B are different, the result is a combination of standing and spinning waves. When the modal amplitudes A and B are identical, the situation corresponds to a spinning wave.
The heat release response function is written as Equation (18).
\begin{equation}Q^{\prime} = F\left[ {p^{\prime}\left( {t - \tau } \right)} \right] = \beta \cdot p^{\prime}\left( {t - \tau } \right) - k \cdot {p{\prime}^3}\left( {t - \tau } \right),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \beta = {\left. {\frac{{dF}}{{dp^{\prime}}}} \right|_{p^{\prime}=0}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k{\kern 1pt} {\kern 1pt} = - \frac{1}{6}{\left. {\frac{{{d^3}F}}{{d{{p{\prime}}^3}}}} \right|_{p^{\prime}=0}}\end{equation}
Considering a general case in which the heat release response is azimuthal non-uniform, the heat release distribution parameter β can be expressed as Equation (19).
\begin{equation}\beta = {\beta _0} \cdot \left\{ {1 + \sum\limits_{m=1}^M {\left[ {{C_m}\cos \left( {m\theta } \right) + {D_m}\sin \left( {m\theta } \right)} \right]} } \right\}\end{equation}
Differential equation for
${\eta _1}$
can be obtained by substituting Equations (15) and (18) into Equation (14), multiplying both sides with
${\rm{cos}}\left( {n\theta } \right)$
and integrating over
$\theta$
ranging from 0 to 2Π. A similar expression for
${\eta _2}$
is obtained by multiplying both sides with sin(n
$\theta$
) before integration. These two differential equations are in the following form. For details of the mathematical development, the reader should consult Ref. (Reference Noiray, Bothien and Schuermans40).
\begin{equation}\begin{array}{l} {\ddot{\eta }_{1} \left(t\right)+\alpha \dot{\eta }_{1} \left(t\right)+n^{2} \eta _{1} \left(t\right)=\beta _{0} \left(1+\frac{C_{2n} }{2} \right)\cdot \dot{\eta }_{1} \left(t-\tau \right)} \\ \\[-7pt]{-3k{\left\{\dot{\eta }_{1} \left(t\right)\left[3\eta _{1}^{2} \left(t\right)+\eta _{2}^{2} \left(t\right)\right]+2\dot{\eta }_{2} \left(t\right)\eta _{1} \left(t\right)\eta _{2} \left(t\right)\right\}\mathord{\left/ {\vphantom {\left\{\dot{\eta }_{1} \left(t\right)\left[3\eta _{1}^{2} \left(t\right)+\eta _{2}^{2} \left(t\right)\right]+2\dot{\eta }_{2} \left(t\right)\eta _{1} \left(t\right)\eta _{2} \left(t\right)\right\} 4}} \right. \kern-\nulldelimiterspace} 4} } \end{array}\end{equation}
\begin{equation}\begin{array}{l}{\ddot{\eta }_{2} \left(t\right)+\alpha \dot{\eta }_{2} \left(t\right)+n^{2} \eta _{2} \left(t\right)=\beta _{0} \cdot \left(1-\frac{C_{2n} }{2} \right)\cdot \dot{\eta }_{2} \left(t-\tau \right)} \\ \\[-7pt] {-3k{\left\{\dot{\eta }_{2} \left(t\right)\left[3\eta _{2}^{2} \left(t\right)+\eta _{1}^{2} \left(t\right)\right]+2\dot{\eta }_{1} \left(t\right)\eta _{1} \left(t\right)\eta _{2} \left(t\right)\right\}\mathord{\left/ {\vphantom {\left\{\dot{\eta }_{2} \left(t\right)\left[3\eta _{2}^{2} \left(t\right)+\eta _{1}^{2} \left(t\right)\right]+2\dot{\eta }_{1} \left(t\right)\eta _{1} \left(t\right)\eta _{2} \left(t\right)\right\} 4}} \right. \kern-\nulldelimiterspace} 4} } \end{array}\end{equation}
The MATLAB solver was used to solve the delay differential equations (Equations (20, 21)). Figure 22 shows the calculation results of the limit cycle amplitudes
${\eta _1}$
and
${\eta _2}$
for different
${C_{2n}}$
and τ. It can be seen that when
${C_{2n}}=0$
, i,e. the heat release distribution is azimuthally uniform, there will be a spinning mode, since the modal amplitude A is equal to B. As the index
${C_{2n}}$
of the heat release azimuthal distribution increases, the modal amplitudes A and B become different; the situation is a combination of standing and spinning waves. When
${C_{2n}}$
becomes larger than its critical value, the standing wave is established. According to this analysis, when the length of the heat release zone is increased, the azimuthally inhomogeneous degree of heat release is increased. Thus, the spinning mode is established for the shorter heat release zone case, while the standing mode is established for the longer heat release case.
The limit cycle amplitudes
${\eta _1}$
and
${\eta _2}$
for different
${C_{2n}}$
and τ.
5.0 CONCLUSIONS
On the basis of air heater characteristics, a new computational model was developed herein, which investigated the acoustics and instabilities of air heaters. Calculation results showed that present simplification ideas and computational models were capable of simulating the critical parameters, acoustics and instabilities inside the air heater. Moreover, these methods had the advantage of numerical cost lower than that of the traditional spray combustion simulation. With appropriate choice of the constants in the present model, expected instability characteristics could be generated, such as limit cycle and grow rate. The present three-dimensional simulations have yielded insight into the consequences of thermo-acoustic coupling. We conclude that the present methods provide a sound basis for the development and calibration of more sophisticated combustion response models.
In a conventional analysis, the combustion is assumed to mainly take place in the so-called rapid combustion zone, and this zone is usually modelled as a disc in the combustor near the injection head. However, in practice, the flame has a spatial distribution. A model that treats the flame as a disc is not close to physical reality. This paper described the effect of spatial heat release distribution on predictions of oscillation frequency and mode. It was found that the frequency and mode shape of oscillations closely depended on the length of the heat release zone. Heat release length had an appreciable effect on longitudinal oscillations. When the heat release zone was located near the faceplate where it was the pressure antinode of the longitudinal mode, the increment of heat release length exhibited an increased tendency toward lower-order modes. Also, as expected, present coupling models could lead to instability for both standing and spinning waves. Heat release length had an appreciable effect on transverse oscillations. Modulation of heat release length may cause bifurcations between standing and turning modes. Lower-order tangential mode occurred for the case having shorter heat release length, and this tangential wave was unsteady, whereas the direction of rotation of these travelling waves was indeterminate, and waves could exist that rotate in one or both directions. It showed a noticeable tendency toward higher-order standing tangential mode with increasing heat release length. Finally, the theoretical analysis of the modal behaviour, i.e. standing or spinning waves, was performed.
Although the present pressure-coupled response function model is preliminary, to match the experimental observations, the calculation results were almost in the linear range. Future work includes investigation of non-linear oscillation using these new methods. We also intend to investigate effects of response model constants on the instabilities.
Acknowledgements
Great thanks to our colleagues from Science and Technology on Scramjet Laboratory for their valuable discussions and comments. This work is supported by the National Natural Science Foundation of China (Nos. 12002356, 12072367, 11802323). The authors acknowledge the financial support of National Natural Science Foundation of Hunan Province, China (No.2020JJ4666).
SUPPLEMENTARY MATERIAL
To view supplementary material for this article, please visit https://doi.org/10.1017/aer.2021.17.
