Characterizing hydrogen-fuelled pulsating combustion on thermodynamic properties of a combustor

Unlike hydrocarbon fuel, hydrogen is ‘green’ and attracting more and more attentions in energy and propulsion sectors due to the zero emission of CO and CO2. By applying numerical simulations, we explore the physics of how a hydrogen-burnt flame can sustain pulsating combustion and its impact on the thermodynamic properties of a standing-wave combustor. We also explain how implementing a heat exchanger can mitigate such pulsating combustion. The dynamic interactions of the unsteady flow-flame-acoustics-heater are examined by varying the mass flow rate ṁH2 and the heating bands’ surface temperature TH. The frequency and amplitude of the pulsating combustion are shown to depend strongly on ṁH2. In addition, varying TH is shown to lead to not only the molar fraction of the combustion species being changed but also the flame-sustained pulsating oscillations being mitigated somehow. Finally, nonlinearity is observed and identified in the unsteady flow velocity and the two heat sources. Fuel cells based on hydrogen combustion are hoped to provide a more environmentally friendly source of energy but still require a lot of development. The authors numerically investigate the physics of pulsating combustion supported by a hydrogen burnt flame.

C ombustion instability is an undesirable phenomenon that occurs frequently in lean premixed combustion systems for power generation or propulsion [1][2][3] . Typically, it results from dynamic interactions between acoustics and premixed turbulent 4 and laminar flames 5 . Combustion instability could be characterized by a loudly 'singing' and violently 'dancing' flame in longitudinal or azimuthal direction (mode). Longitudinal-type instability has been intensively studied 6 . The dominant thermoacoustic instability is experimentally shown to switch by observing a shift from a low-to a high-frequency 'tone'. When combustion instability occurs 7,8 , it can cause severe vibration, overheating, flame flashback, and decreased performance 9 . In severe cases, the engines, propellant systems and structures could be damaged.
Intensive research has been conducted on hydrogen-involved combustion instability. Matsuyama et al. 10 conducted LES (large eddy simulations) on a single-element atmospheric combustor. It is found that the primary driving mechanism of 1k Hz combustion instability is the acoustically coupled pulsating motion of the inner H 2 /O 2 flame or periodic ignition of H 2 /O 2 mixture. Recently, Hemchandra et al. 11 discovered that there are two different mechanisms driving combustion instability. One is a strong coupling between acoustics and hydrodynamic modes. The other one is a weak coupling resulting in semi-openloop forcing of the flame by a self-sustained hydrodynamic mode. Urbano and Selle 12 analyzed transverse combustion instabilities in a reduced-scale rocket motor. The interaction between acoustics and vorticity is found to be the main damping mechanism for coaxial H 2 /O 2 flame-sustained instability. Injector-driven combustion instabilities are experimentally observed in a hydrogen/oxygen rocket motor 13 . The observed instabilities are a result of the interaction between the injector resonant frequencies and the combustion chamber resonant frequencies.
Combustion instability may be mitigated by applying two classical means. One is passive control. Kim et al. 14 experimentally evaluated the attenuation effect of sponge-like porous medium on combustion instability in a lab-scale swirl combustor. It is found that up to 40% of pressure attenuation is achievable. A perforated plate/liner 15 is analytically studied on damping combustion instability in an annular combustor 16 . Tao et al. 16 show that a cooling/bias flow can enhance the acoustic damping effect of such perforated plates/liners. The damping performance is found to be optimized at Strouhal number near 0.3. Similar stabilization effects are observed on perforated slits on a modelled combustor 17 . Wu et al. 18 actively tuned the acoustic damping of a Helmholtz resonator to mitigate combustion instability in a Y-shaped Rijke tube. Approximately 60 dB sound pressure level (SPL) reduction is achieved by actively tuning the resonator's back-wall. Similar implementation of a Helmholtz resonator 19 is conducted on a Rijke tube 20 .
The other approach is active control using a developed controller. Phase-shift or PI (proportional-integral) controller is the typical one that is applied to suppress combustion instabilities by using loudspeakers 21,22 . An excellent review on active control of combustion instability can be found in refs. 23,24 . In addition to the classical active and passive control approaches, alternative means have been proposed and tested. For example, Jocher et al. 25 experimentally confirm that combustion instability could be controlled in open-loop configuration by applying a magnetic field via reducing the growth rate of flame perturbations. Henderson and Xu 26 experimentally study the effect of applying direct-current electric field on damping thermoacoustic instability in a horizontal Rijke tube. Approximately 20 dB SPL reduction is achieved. Combustion instability is mitigated by actuating a swirler in a premixed methane-fuelled combustor 27 . It is found that creating a higher-intensity turbulence and a swirl flow can alter the flame position and structure, and so the combustor stability behaviours.
Fuel composition is found to play an important role in generating and suppressing combustion instability. Diao et al. 28 experimentally study the flame-acoustics interaction in a rectangular H 2 /O 2 shear-coflow combustor with CH 4 actively blended. It is found that heat release oscillations are significantly dampened, when gaseous methane is blended with gaseous hydrogen. The effect of blending methane with hydrogen on combustion instability is experimentally studied 29 on a swirl combustor. It is found that increasing hydrogen concentration (up to 30% in volume) in the fuel leads to enhanced responses of the flame at non-resonant frequencies. Choi and Lee 30 found that when hydrogen is added into syngas both acoustic pressure, unsteady heat release respond sensitively and nonlinearly. Similar combustion instability studies of H 2 /CO/CH 4 syngases were conducted on a partially premixed gas turbine combustor 31 . H 2 contents are identified to play an important role in producing higher-order (higher frequency) instability mode and modeshifting. Acoustic measurements of a lean premixed combustor with hydrogen and methane or propane fuelled confirm the key role of hydrogen in producing higher sound -pressure-level instability 32 . The effects of fuel composition on a premixed gas turbine combustion are well reviewed in ref. 33 , with emphasis being placed on combustion system stability 34,35 and emissions.
To the best knowledge of the authors, there are few reports in the literature discussing the physics and mechanisms underlying hydrogen-fuelled pulsating combustion, and its impact on the thermodynamic properties of a standing-wave combustor. In addition, there is a lack of alternate but effective means for preventing or mitigating combustion instability. In this work, an acoustically open−open premixed combustor fuelled with hydrogen-air mixture and a constant-temperature heat exchanger is numerically studied to shed light on the dynamic flame−acoustics−heater interaction. The hydrogen mass flow rate _ m H 2 and the surface temperature T H of the heating bands are evaluated. It is found that varying T H can lead to partial or complete change in the molar fraction of the combustion species and the mitigation in the flame-sustained pulsating oscillations. Finally, nonlinearity is identified and observed on these heat sources.

Results
Effect of mass flow rate _ m H 2 . As the mixture of the hydrogen and air is burned, a conical shaped flame is confined in the combustor. Figure 1 shows the turbulent reaction rate as the hydrogen mass flow rate _ m H 2 is set to three different values. It can be seen that the flame length and the surface area are dramatically increased, and so is the heat release. As it is well known that unsteady heat release is an energy-efficient monopole-like sound source, acoustic pressure disturbances are generated with a small amplitude initially. However, whether the acoustic disturbances will grow into limit cycle pulsating combustion oscillations needs to be examined.
The effect of the flame on producing pulsating pressure oscillations is shown in Fig. 2. Figure 2a shows the time evolution of the acoustic pressure fluctuations at the midpoint of the combustor in the axial direction, as the hydrogen mass flow rate _ m H 2 is set to 1.0 × 10 −5 kg s −1 . It can be seen that there are smallamplitude pressure disturbances initially. However, such distances grow quickly into a limit cycle. Figure 2b shows the corresponding phase diagram of the pressure fluctuations. The circle shape characterizes the limit cycle oscillations in phase plot. Figure 2c illustrates the phase diagram of the acoustic pressure, as the mass flow rate _ m H 2 of the hydrogen is set to 5.0 × 10 −6 kg s −1 . It can be seen that a limit cycle is also produced. However, the amplitude is approximately 50% of that at _ m H 2 ¼ 1:0 10 À5 kg s À1 . Further decreasing _ m H 2 leads to a reduction in the amplitude of the pulsating oscillation, as shown in Fig. 2d.
Unlike hydrocarbon fuels, hydrogen is a nonhydrocarbon that exhibits the largest mass diffusivity features. Figure 3 compares the profiles of the reactants and combustion product species in both axial and radial directions, as the hydrogen mass flow rate _ m H 2 is set to three different values. The red colour curves show the major hydrogen distribution. Along the axial direction (see Fig. 3a-c) the mole fraction begins to decay sharply; complete depletion of the hydrogen is achieved after the reaction zone. In the radial direction, this finding is consistent with the variation of _ m H 2 . The mole fraction profile of H 2 O is observed to vary dramatically in the axial direction, as the hydrogen mass flow rate is changed (see the blue colour curves in Fig. 3a-c). As _ m H 2 is decreased from 1.0 × 10 −5 to 1.0 × 10 −6 kg s −1 , the separation between the preheat and reaction zones becomes more sharp. In other words, the wave/peak thickness, i.e. the maximum molar fraction region, is reduced. The profile of O 2 is almost the same in the axial direction, as _ m H 2 is varied. It is injected from the Bunsen burner, and so its mole fraction is increased gradually from the pre-reacting zone to the reacting zone and the after-reaction zone, and is finally 'saturated'. As far as the radial profiles (see Fig. 3d-f) are concerned, they are similar in shape for Fig. 3d, f, and in the growth rate for Fig. 3e. The shape profiles are similar especially near the inner surface of the tube. The maximum mole fraction of H 2 is reduced, and so is the thickness of the wave/peak of the water (H 2 O) profile. The hydrogen mole fraction begins to decrease from the point r/R = 0 and it becomes much sharper, as _ m H 2 is reduced to _ m H 2 ¼ 1:0 10 À6 kg s À1 . In general, H 2 , O 2 and H 2 O are all a function of the spatial distance through the flame. Both water and hydrogen exist in the preheat zone. Note that more complicated 21-step chemical reactions are studied, as discussed in Supplementary Note 1. The new results with the 21 chemical reaction mechanisms/steps are provided in Supplementary Table 1. Comparison is made between the multiple-step and the much-simplified one-step reaction in terms of the temperature and velocity contours, and the combustion-driven pressure and the heat of reaction pulsating oscillations. Interested readers can find more details in Supplementary Figs. 1 and 2, and insightful discussion in Supplementary Note 1.
As the mixture of the hydrogen and air is burned, the instantaneous temperature at the midpoint of the combustor in the axial direction is monitored in real time. This is shown in Fig. 4 with _ m H 2 being set to three different values. It can be seen from Fig. 4a-c that the temperature is decreased dramatically, as _ m H 2 is decreased. The temperature fluctuates dramatically before the limit cycle oscillations are produced by the flame, i.e. t ≤ 0.5 s. However, with limit cycles being completely generated at t ≥ 1.0 s, the temperature seems to be 'saturated'. To gain insight into the fluctuations of the temperature, a phase diagram of the instantaneous temperature T(t) is plotted, as shown in Fig. 4d-f. It can be clearly seen that the gradient of the temperature, i.e. dT/dt, is non-zero. This means that the temperature is periodically oscillating around its mean value. Furthermore, the oscillation amplitude is dramatically reduced with a decrease in _ m H 2 ( Fig. 4d-f). Finally, the fluctuating trend of the temperature reveals that the heat-to-acoustics energy conversion is nonirreversible.
The mode shapes of the acoustic pressure fluctuations are determined, and are discussed in Supplementary Note 2 (shown in Supplementary Fig. 3). Comparison is then made with the case of using CH 4 as the fuel. It is shown that the maximum pressure amplitude is observed near the mid-point of the combustor in the axial direction. A pressure node is expected at the combustor outlet, no matter whether CH 4 or H 2 is fuelled.
Effect of the heating bands. Heating bands with a constant surface temperature are implemented, aiming to dampen and/or prevent the onset of the flame-sustained pulsating oscillations. Time evolution of the pressure fluctuations as the temperature of the heating bands is increased from 300 to 1300 K at t = 2.0 s is shown in Fig. 5a, as _ m H 2 is set to three different values. It can be seen that as the hydrogen mass flow rate is greater than _ m H 2 ¼ 5:0 10 À6 kg s À1 , increasing the temperature of the heating bands leads to a mitigation in the limit cycle oscillations. However, when the hydrogen mass flow rate is reduced to _ m H 2 ¼ 1:0 10 À6 kg s À1 , actuating the heating bands is found to successfully prevent the onset of the limit cycle oscillations, as   denoted by the solid green colour cure in Fig. 5a. To quantify the damping effect of the heating bands, the instantaneous acoustical energy E a (t) per unit volume of the combustor is determined: where p′(t) is the acoustic pressure fluctuation, ω is the oscillation frequency, γ = 1.4 is the specific heat ratio, and p 0 = 1.01325 × 10 5 Pa is the mean pressure. Note that the total acoustical energy E tot a at a given time can be determined by integrating E a (t) over the entire combustor volume. If a 3D cylindrical coordinate is considered, then E tot a ðtÞ ¼ dx. Since the combustor volume is specified with the constant cross-sectional area and the axial length L, E a can be applied to characterize the intensity of combustion oscillations.
Time evolution of the acoustical energy E a is shown in Fig. 5b as _ m H 2 is set to three different values. It can be seen that as _ m H 2 ¼ 1:0 10 À5 kg s À1 , the acoustical energy E a is reduced by 61.7% from t = 2.0 s to t = 4.0 s. Approximately 81% acoustical energy reduction is achieved at _ m H 2 ¼ 5:0 10 À6 kg s À1 . At _ m H 2 ¼ 1:0 10 À6 kg s À1 , the pressure fluctuations are successfully prevented from growing into limit cycles, as observed for _ m H 2 ¼ 5:0 10 À6 kg s À1 and _ m H 2 ¼ 1:0 10 À5 kg s À1 . This confirms that the heating bands are applicable to dampen flamesustained pulsating oscillations.
Time-frequency analysis is conducted to study the heat of reaction of the hydrogen-fuelled flame and the total heat fluxes of the heating bands before and after changing the temperature of the heating bands (Fig. 6). It can be seen that the nonlinearity is associated with the heat release rates of the flame and the heating bands, due to the comparable frequency peaks. Before the heating bands' temperature is changed, i.e. t ≤ 2 s, the heat release of the flame is associated with multiple modes at different frequencies. The dominant mode frequency is approximately 204 Hz (see Fig. 5a). However, the harmonic modes are excited as the heating bands' temperature is increased to 1300 K at 2.0 < t ≤ 4.0 s. Nonlinearity characteristics of the combustion system are further identified by conducting the frequency spectrum analysis of the flow velocity. This is shown in Fig. 7, as the hydrogen mass flow rate _ m H 2 is set to three different values. Figure 7a-c shows the time evolution of the instantaneous velocity. It can be seen that as _ m H 2 is decreased, with no change in the temperature of the heating bands (t ≤ 2.0 s), the mean velocity is reduced from 3.1 to 1.53 m s −1 . The amplitude of the velocity fluctuation is also decreased dramatically. Figure 7d-f illustrates the frequency spectrum of the velocity fluctuation, before and after the temperature of the heating bands is changed. It can be seen that 4% of the dominant frequency shift occurs from t ≤ 2.0 to t > 2.0 s. A similar frequency shift is observed from 203.0 to 177.6 Hz, as the mass flow rate _ m H 2 is decreased from 1.0 × 10 −5 to 1.0 × 10 −6 kg s −1 . However, the dominant frequency shift is approximately 12.7%. This is most likely due to the fact that the H 2 mass flow rate and the total heat release are reduced. The mean temperature and the speed of sound in the combustor are also decreased. To maintain the standing-wave mode shape and to satisfy the inlet and outlet boundary conditions, the dominant frequency of the H 2 -excited oscillations is thus decreased. In general, the flame-excited pulsating oscillation signature is found to depend strongly on the hydrogen mass flow rate.
The effect of the heating bands on combustion species is also examined, as shown in Fig. 3. Comparison is made before and after the surface temperature T H of the heating bands is increased from 300 to 1300 K at t = 2.0 s. It can be seen from Fig. 3a-c that the mole fraction profiles of the reactants of H 2 and O 2 , and the product H 2 O at T H = 1300 K are quite similar to those at T H = 300 K in the axial direction. Furthermore, the thickness of 'the wave/peak' profile is also the same. This is most likely due to the negligible effect of the mean heat release of the heating bands on the hydrogen-fuelled flame, as they are placed far away downstream of the combustor. As shown in Fig. 3d-f, the radial profiles of the mole fraction of O 2 look quite similar, before and after T H is changed. This is most likely due to the fact that the inner burner is associated with a mixture of hydrogen and air. In addition, there is a simultaneous flow of air injected from the bottom of the cylindrical tube (see Fig. 8). The mole fraction profiles of H 2 and H 2 O, especially the maximum molar fraction of the species at r/R = 0, are however different. This may be due to the fact that the periodic combustion oscillations are dampened with the actuation of the heating bands, which affects the hydrogen and H 2 O distribution along the radial direction. The oscillation intensity and frequency may contribute to the radial profiles variation of H 2 and H 2 O.
As revealed by Rayleigh 1,2 , limit cycle pulsating combustion occurs due to the unsteady heat release, and acoustic pressure fluctuations are in phase. In order to shed light on the dynamic flame−flow−heater interaction, Rayleigh index angle θ RI is determined (Supplementary Note 3) to check whether the unsteady heat release from the flame Q′ f and the heating bands Q′ H are in or out of phase with the acoustic pressure fluctuations. θ RI is defined as the phase difference between the heat release and the pressure fluctuation at the frequency ω at any instant. Supplementary Figures 4 and 5 show the time evolution of the phase difference between the acoustic disturbances and the unsteady heat release from the flame and the heater, respectively.

Discussion
We have numerically shown how hydrogen-fuelled flame can sustain self-excited pulsating oscillations, and its impact on the thermodynamic properties of a standing-wave combustor with heating bands confined. The physics behind the flame−acoustics −heater−flow interaction is numerically examined. It is found that the hydrogen mass flow rate _ m H 2 determines not only the pulsating combustion characteristics but also the molar fraction profiles of the combustion species. Approximately 13% of the dominant frequency shift is observed in the dominant mode. Furthermore, the combustor is found to be highly nonlinear in terms of the acoustic velocity and heat of reaction Q f from the flame and the heat flux Q H from the heating bands. Finally, varying the surface temperature T H of the heating bands is shown to be able to successfully prevent the onset of limit cycle pulsating oscillations excited by the hydrogen-fuelled flame. By conducting further analyses on the mode shape and by determining the Rayleigh index angle it has been confirmed that, when limit cycle pulsating combustion occurs, the combustor is a standing-wave one. Rayleigh index angle reveals that the unsteady heat release constructively interacts with the acoustic disturbances in the flame, while destructive interactions are observed on the heating bands. This is consistent with the Rayleigh criterion and explains why the heating bands can be applied to dampen pulsating combustion oscillations. In addition, it is interesting to show that increasing the surface temperature of the heating bands leads to the combustion species profiles being varied. The present work sheds light on the fundamental physics and mechanism underlying hydrogen-fuelled combustion instability. It also opens up an alternative control means to dampen or to prevent the onset of the pulsating oscillations. The overall view of the 2D mesh generated by using Ansys ICEM CFD is shown in Fig. 8b. The mesh consists of 25,868 quadrilateral elements, 372 elements along the axis of the tube and 76 elements along the radius. Near the heating bands and the Bunsen burner, a finer grid is applied, as shown in Fig. 8b. The mesh is chosen on the basis of our mesh-independence study as discussed in Supplementary Note 4. The mesh-independence results are provided in Supplementary Figs. 6 and 7. Note that a 2D combustor is modelled for simplicity and for reducing the computational time and cost. However, a 3D modelled combustor is also investigated and compared with the 2D one for completeness, as discussed in Supplementary Note 5. This is shown in Supplementary Fig. 8. Approximately 16% difference is found between the pressure pulsation amplitudes predicted from the 2D and 3D combustors. However, the computational time is increased by 400%.
With the flame confined, there is an unsteady heat release and a mean flow results from the convection effect. In addition, unsteady heat release is an energyefficient monopole-like acoustics source 36,37 . Acoustic distances are produced and propagated along the combustor. The interaction between the acoustic disturbances and the heating bands leads to further unsteady heat release. When the temperature of the mean flow is higher than the surface temperature of the heating bands, they absorb the heat released from the flame. Otherwise, heat is produced by the heating bands and transferred to the mean flow. The heat transfer process in our modelled combustor is unsteady. Thus, the current modelled combustor involves the dynamic interaction between the premixed flame, the unsteady flow, the acoustics and the heating bands.
Governing equations. The flame−acoustics−flow−heater interaction problem 38,39 is modelled by solving certain equations governing fluid dynamics and the chemical combustion reaction. All these equations need to be solved simultaneously and some of them are coupled, e.g. the fluid velocity will determine how the combustion species 40 will spread through the computational domain and at which particular locations in the numerical domain there might be more or less chemical reaction taking place.
The governing equations of mass, momentum, species fraction, energy conservations in the combustor are given as where χ k is the local mass fraction of the kth species, J k is the diffusion flux of the kth species, R k is the net rate of the kth species by chemical reaction, δ ij denotes the Kronecker delta function and δ ij = 1 when i = j; otherwise δ ij = 0. The turbulent viscosity μ t is introduced to model the effect of turbulence on the flow according to the Boussinesq hypothesis. This unknown needs to be determined via an additional modelling, for which the RNG k−ε model in ANSYS Fluent is employed. The turbulent viscosity is computed from k, the turbulent kinetic energy, and ε, the turbulent dissipation, as given in the following expression: The turbulent kinetic energy k and the dissipation rate ε are determined from the transport equations of similar structures to the governing equations for the fluid flow. They are given as and where G k and G b denote the production of turbulence kinetic energy resulting from the mean velocity gradient and buoyancy effects, respectively. Y M describes part of the overall dissipation rate due to the fluctuation dilatation in compressible turbulence. The eddy-dissipation k-ε model is used to determine the production term R k in Eq. (4). This production term is computed with the following expression: where MW k is the molecular weight of the species k, v′ k denotes the stoichiometric coefficient for the reactant kth species and v′′ j is the stoichiometric coefficient for the product jth species in the chemical reaction. Further details can be found in ref. 41 or the Fluent documentation, which contains more information about the modelling of the governing equations. A more complex turbulence model 42 , e.g. RSM (Reynolds stress model), could have been applied. For comparison, RSM (involving seven equations) is also implemented and evaluated as described in Supplementary Note 6 (Supplementary Figs. 9 and 10). It is shown that the hydrogen-driven pulsating combustion 43 is successfully generated by using either the RSM or the k-ε model. To reduce the computational time and cost, we chose the k-ε turbulence model in the current studies.
Model settings. The essential properties of the computational fluid dynamics (CFD) model are set as follows: (1) a pressure-based coupled solver; (2) unsteady simulations with a time step of Δt = 1.0 × 10 −4 s and a second-order discretization; (3) a spatial discretization: second-order pressure discretization, and third-order MUSCL scheme for all equations, except the turbulence dissipation (second-order upwind); (4) hydrogen-air mixture (air components: nitrogen (76.8%), oxygen (23.2%) in mass); (5) species transport model, volumetric reactions, eddydissipation turbulence-chemistry interaction; (6) a turbulence model: k−ε−RNG; and (7)  temperature (setting T = 300 K for 0.0 ≤ t ≤ 2.0, T = 1300 K for 2.2 s < t, controlled via UDF (user-defined function); (d) the Bunsen burner tip is a mass flow inlet; (e) the bottom end of the open tube is a pressure inlet (i.e. 0 Pa); (f) the top end of the open tube is a pressure outlet (i.e. 0 Pa). The numerical simulation is started with a flow field as the initial solution, which includes a steady flame, i.e. a flame that is developed in a steady-state simulation. The fuel-air ratio 40 is set to 1.0 for all current simulations. Detailed information is summarized in Supplementary  Table 2.
Data post-processing. Short-time Fourier Transform (STFT) is used to analyse the data regarding the acoustic pressure, velocity, and the reaction of heat and heat fluxes. STFT is developed based on the DFT (Discrete Fourier Transform). Let us consider z(n) as an N-point acquired pressure signal in the time domain t(N) = NΔt, and let W N = exp(−j2π/N) and j ¼ ffiffiffiffiffiffi À1 p . The DFT of z(n) is Z(n) = DFT{z (n)} is defined as Frequency-domain figures are obtained with Δt = 1/10,000 s and N = 8192 and a Hamming window.

Data availability
The data that support the findings of this research are available from the corresponding author on request.

Code availability
The code (case files) and data sets generated by using ANSYS FLUENT 15.1 during this study are available from the corresponding author on request.