An Eﬃcient Method to Reproduce the Eﬀects of Acoustic Forcing on Gas Turbine Fuel Injectors in Incompressible Simulations

Previous studies have highlighted the importance of both air mass flow rate and swirl fluctuations on the unsteady heat release of a swirl stabilised gas turbine combustor. The ability of a simulation to correctly resolve the heat release fluctuations or the flame transfer function (FTF), important for thermoacoustic analysis, is therefore dependent on the ability of the method to correctly include both the swirl number and mass flow rate fluctuations which emerge from the multiple air passages of a typical lean-burn fuel injector. The fuel injector used in this study is industry representative and has a much more complicated geometry than typical premixed, lab-scale burners and the interaction between each flow passage must be captured correctly. This paper compares compressible, acoustically forced, CFD (computational fluid dynamics) simulations with incompressible, mass flow rate forced simulations. Incompressible mass flow rate forcing of the injector, which is an attractive method due to larger timesteps, reduced computational cost and flexibility of choice of combustion model, is shown to be incapable of reproducing the swirl and mass flow fluctuations of the air passages given by the compressible simulation as well as the downstream flow development. This would have significant consequences for any FTF calculated by this method. However, accurate incompressible simulations are shown to be possible through use of a truncated domain with appropriate boundary conditions using data extracted from a donor compressible simulation. A new model is introduced based on the Proper Orthogonal Decomposition and Fourier Series (PODFS) that alleviates several weaknesses of the strong recycling method. The simulation using this method is seen to be significantly computationally cheaper than the compressible simulations. industrial context, this improved computational efficiency allows for greater exploration of the design space and improved combustor design.


Introduction
Combustion instabilities continue to be a major hurdle in the development of efficient, low emission, gas turbines [1]. These instabilities involve interaction between the acoustic field, flow field and combustion. The acoustic field is sensitive to the entire combustor and surrounding component geometry and can be analysed using a number of different methods such as a low order acoustic network model [2], or Helmholtz solver [3]. However, a crucial part of the process is the response of the flame to the unsteady air supply. This is often characterised by a linear flame transfer function (FTF) or non-linear flame describing function (FDF) (see for instance [1,4]). Studies on premixed turbulent swirl flames by Palies et al. [5] have showed that acoustic fluctuations have two main effects on heat release. When approaching the fuel injector, the acoustic wave induces an acceleration of the air within the injector body. This acceleration induces vorticity at the exit of the injector leading to a vortex ring that is convected towards the flame. The approaching vortex ring causes the flame to roll up, modifying the flame area and therefore the heat release rate. The same axial velocity fluctuations that cause the production of the vortex ring also interact with the swirl vanes within the injector creating an azimuthal velocity fluctuation. The azimuthal fluctuation can impact upon the swirl number within the combustor. The swirl number [6] is formally defined as the the ratio between azimuthal momentum flux and axial thrust, however, this requires both the pressure and axial velocity to be known and as such the pressure term is usually excluded leading to: The swirl number also contributes to heat release through its contribution to turbulent flame speed and the cone angle of the flame [7].
It was therefore shown by Palies et al. [5] that fluctuations in turbulent burning velocity could be linked to both the axial and azimuthal velocity fluctuations. The importance of the azimuthal and axial velocity fluctuations on the FTF has also been investigated by Acharya et al. [8].
It should be noted that the phase as well as magnitude of the response to velocity and swirl are important and they could combine either constructively or destructively. In the case of a lean-burn aero engine fuel injector, the geometry is considerably more complicated than in [5] due to the presence of multiple airflow passages. The relative response of these passages will be important in determining the overall heat release.
Given this, it becomes crucial when studying the thermoacoustic response of a fuel injector that the axial velocity fluctuations and the swirl number at the exit of the injector are captured accurately. As an important part of understanding the theromoacoustic cycle, the FTF may be found using experiments [5,9,10] but is often studied with numerical simulations [11][12][13] due to the greater detail that can be extracted. As the FTF is independent of fluctuation amplitude, there is a temptation to use an incompressible numerical method to do this due to its greater computational speed and flexibility of combustion modelling.
Four such examples are the works of Febrer et al. [14], Chatelier et al. [15], Biagioli et al. [16] and Han et al. [17] where the inlet velocity is assumed to fluctuate sinusoidally around a mean value. Unfortunately, this approach may not include the true response of the injector to acoustic forcing. This is due to the inability of incompressible CFD to predict correctly the acoustic impedance of the injector as a whole or its individual passages. This is defined as This complex valued quantity includes the resistance (R) that defines the in phase response of the fuel injector and the reactance (X) which is the out-of-phase response. The acoustic impedance can be found using experiment [18][19][20] and this experimental data has been used to validate the ability of compressible CFD to predict the correct response of fuel injectors to acoustic forcing for a range of designs across the frequency range of interest [18,[21][22][23]. These references show that the magnitudes of reactance and resistance are similar at the frequencies of interest and therefore a compressible methodology is required. Each flow passage in this injector geometry has a different acoustic impedance and each swirler will respond differently to the incoming acoustic wave. Capturing these individual passage responses accurately and how all the passage jets interact further downstream is crucial for predicting how the flame will respond to acoustic perturbations. This paper presents a direct comparison between the predicted injector response to correct acoustic forcing in compressible CFD and mass flow forcing in incompressible CFD as might typically be used in FTF predictions. This will be used to show that simply forcing the mass flow does not give the correct mass flow and swirl variation for the different passages, which would in turn affect the accuracy of FTF calculation. As mentioned previously, there are some clear advantages to using incompressible CFD methods for FTF calculations and to this end a method is presented to use the response predicted by a compressible simulation to provide accurate boundary conditions for an incompressible simulation. By truncating the computational domain at the exit of the injector flow passages, 2D flow field data extracted from the acoustically forced simulations may be applied to the incompressible domain at each timestep as in the strong recycling method [24]. This procedure has distinct disadvantages however as the data for each inlet at each timestep must be loaded during the simulation runtime. If the timestep is altered than the data must be interpolated and if the simulation is extended then a blending region must be used to prevent discontinuity when the time series of data is reused. The POD and Fourier Series (PODFS) method is introduced to do this in a far more efficient fashion by first compressing the inlet data using the proper orthogonal decomposition (POD) and converting the temporal components into a Fourier series that guarantees continuity in time and flexibility of the timestep. Section 2 covers the geometry and test conditions, Section 3 describes the numerical methodology for the incompressible and compressible solvers, and the different forcing strategies, including the new PODFS method. Section 4 gives the results of different forcing methods and clearly shows that the PODFS is the best compromise between accuracy and speed by comparing the fluctuations of mass flow rate and swirl number at the exit of the injector passages and at the injector exit plane.

Injector Geometry and Test Conditions
The injector chosen in the study is representative of a modern radially staged lean burn fuel injector for large civil engines (see Fig. 1). It consists of three flow passages (A, B and C) with passage A feeding air to a pilot flame that is maintained at a relatively rich condition to ensure flame stability. The main staged flame is fed by an air-blast atomiser located between the B and C flow passages which provide air to the mains system. Swirl is induced in passages A,B and C to induce recirculation zones for the flame however the swirl direction is reversed for passage A to increase mixing.
The injector flow was forced at St = 0.3 at an equivalent magnitude to a p = 0.003p 0,in acoustic fluctuation using three different methods: i) correct acoustic forcing using a compressible solver and appropriate characteristic based boundary conditions, ii) bulk mass flow forcing in an incompressible simulation and iii) an incompressible simulation using a new, PODFS method, to introduce the correct relative mass flow and swirl fluctuation at each passage.
Lengths, volumes and areas have been normalised by the injector diameter D while velocities have been normalised by U ref defined as: where p 0,in is the upstream total pressure, p out is the downstream static pressure, and ρ is the upstream air density. Time scales have been non-dimensionalised using: In this study only air is considered as the working fluid. Under reacting conditions fuel is injected from both the pilot atomiser and mains airblast atomiser as a liquid spray. Due to the time it takes for the fuel to evaporate, the flow through the injector has a mixture fraction very close to zero and so no reaction takes place within the injector passages.

Acoustic forcing
The full computational domain is shown in Fig. 2. The acoustic forcing is achieved by applying characteristic boundary conditions to the inlet and outlet using the methods further described in [22] using solvers provided by OpenFOAM. The mean pressure drop across the injector was set to 3 % and first iterated to a pseudo-steady state condition with the PIM-PLE algorithm using a timestep of 1 × 10 −1 t * and a first order backwards Euler method, before being switched to a compressible PISO solver [25] with a timestep of 1 × 10 −3 t * using second order backwards time stepping and characteristic boundary conditions. The Mach number based on U ref is Ma = 0.2 and the Reynolds number based on the same velocity and the injector diameter is Re = 3.8 × 10 5 . The incoming characteristic wave magnitude was set to zero at the hemispherical upstream boundary resulting in a nonreflective condition while the downstream exit incoming wave amplitude was set such that it imposed a St = 0.3 acoustic wave into the domain. This frequency along with the operating conditions were chosen as they are representative of real combustion instabilities in industrial geometries.The forcing resulted in an area averaged fluctuation of axial velocity of |u x | = 0.05U ref . When the forcing is applied using an acoustic signal, the direction of forcing is important [26] because acoustic waves can propagate either upstream or downstream through the fuel injector leading to different fluctuations of axial velocity and swirl number. In this study and in [22], the injector is forced from downstream. Spatial discretisation was provided by a second order TVD scheme for velocity with a first order upwind scheme being used for the turbulent quantities. This method has been validated against experimental measurements of the acoustic impedance, as described in [22], where it was seen to give excellent agreement for a range of frequencies from St = 0.2 to St = 0.5. As such this can be used as the reference case. All simulations use a URANS methodology employing the k − ω SST turbulence closure [27].
The mesh is unstructured and generated using ICEMCFD, it contains 9.2 million cells. The mesh density is set to 0.027D in the first 2.5D of the duct and 0.0064D in the main swirl passages, the cells in passage A were further limited to 0.0032D. The mesh was first generated of the hemispherical plenum, injector and the first 2.5D of the duct using the Delaunay method. After smoothing, the mesh was extruded from the outlet to form the downstream duct. The extrusion was performed using a geometric expansion ratio of 1.1, however in order to ensure that the acoustic waves are resolved accurately in the duct, the duct mesh size was limited to 0.046D. This numerical method has been tested over a range of geometries and test conditions as detailed in [18,21,22] where it was shown to accurately predict the acoustic impedance for each geometry tested.

Mass flow forcing
Mass flow forcing is achieved simply by scaling the inlet velocity in an incompressible simulation by an appropriate factor such that the mass flow through the domain increases or decreases appropriately. This method is often used to approximate the effects of acoustic forcing on a flame to derive the FTF or FDF [14,16,17]. The computational domain is altered such that the hemispherical upstream plenum used in the compressible acoustically forced simulation (see [22]), is replaced with a cylinder of diameter 3.9D, which allows a planar mass flow boundary condition to be used, and the outlet duct shortened to 6.5D. Figure 3 shows the updated domain. At each timestep the inlet mass flow is set by: The values ofṁ in andṁ in can be extracted from an acoustically forced compressible simulation or from an experiment. In this case they are taken from the total mass flow The mesh is designed to be as close as possible to the original compressible mesh in the vicinity of the injector with the same control parameters and procedure. As this mesh does not require the resolution of acoustic waves in the duct, and only serves to represent the mean pressure field, the mesh spacing can be relaxed significantly away from the injector. This was achieved by allowing the geometric expansion of the cells to continue in the xdirection along the duct. This results in a mesh size of 0.175D at the exit.
The velocity equations were solved using the incompressible PRECISE-UNS [28] solver with a second order accurate linear upwind flux calculation with a van Albada limiter. The pressure equation was solved using second order central flux calculation and the turbulent quantities solved using a linear upwind method. The solution was stepped forward in time using a second order backward scheme with a time step of 2 × 10 −2 t * . The k − ω-SST model was used for turbulence closure. These values are considered current best practice for this solver and were also used for the recycled and PODFS simulations described below. The interest in using this solver is that it has access to specialised combustion models that have been used to model academic and industrially relevant combustors [28].

PODFS model
The PODFS (Proper Orthogonal Decomposition Fourier Series) model has been developed here to overcome several weaknesses of the strong recycling method. The strong recycling method extracts the velocity fields of a donor simulation and applies it directly to the inlet of a new simulation. If the two simulations do not have the same time step then the data must be interpolated in time and if the second simulation is longer than the first then the inlet data must be repeated usually with a blending region to prevent discontinuities and to reduce non-physical periodicity [24]. In the PODFS method, the snapshot POD method [29] is applied to the data collected on the same planes within a donor simulation and used as a filter in order to compress the raw data. This filtered data then consists of a set of spatial modes that are independent of time (φ j (x)) and a set of temporal modes that are independent of space (a j (t)).
The periodic nature of the fluctuating boundary conditions may make a model based on the dynamic mode decomposition (DMD) [30] more suitable than POD. POD has been chosen here because it is energetically optimal in the sense that the largest amount of energy is captured with the smallest number of modes, however as DMD modes fluctuate at a single defined frequency, they may be more suitable in cases where the inlet fluctuates at one frequency only. Models based on DMD modes do contain challenges: It may be the case that a single DMD mode is capable of representing a very large proportion of the fluctuating flow. Such a mode will have an oscillation frequency equal to the forcing frequency and a growth rate of zero. Whilst it is obvious that the growth rate should be zero, computing a mode with a growth rate of zero (or sufficiently close to zero) will take a very large number of snapshots because the convergence of modes in DMD is typically slower than for POD. The mode growth rate can be artificially set to zero but then the model accuracy could suffer as multiple computed modes would be contributing to the inlet fluctuation (some with positive and some with negative growth rates). POD also has the advantage that even if the modes are badly converged, using a Fourier series to represent their temporal evolution guarantees that the model solution will remain bounded whereas only DMD modes with zero growth rate can be used for long integration times.
The POD spatial modes in the PODFS model are constructed from the velocity, k and ω fields such that The raw extracted data is normalised leading to U * in defined as: where denotes a spatial average, U in = [u in , k in , ω in ] and the over line represents the temporal average. These values are then used to calculate the covariance matrix according to the snapshot POD method. The eigenvalues of the covariance matrix (λ) describe the relative magnitude of each POD mode. If the POD modes are ordered according to the size of their eigenvalue the eigenvalue energy (E λ (N P )) is defined as: with N P ≤ N S and is chosen such that E λ is sufficiently close to 1. The POD approximation of the normalised inlet variables is then given by: Where k is the kth timestep. The spatial modes are invariant in time and ready for use in the model, but, the temporal modes are of a finite length of time and at a set time interval. In order to overcome these weaknesses of the strong recycling method, the temporal modes can be interpolated using a periodic analytical function. In this case a Fourier series is chosen to reproduce each temporal POD mode. An ordinary Fourier representation would result in the following approximation of each temporal mode: with the coefficients defined by: N F is the number of Fourier terms to be used in the approximation and as N F approaches N S the error between the Fourier representation and the original temporal POD modes will decrease, however not all Fourier terms contribute equally to the error. The standard Fourier reconstruction is therefore modified to enable an arbitrary selection of whatever Fourier terms reduce the error the fastest: The function f (j) allows for an arbitrary ranking of each term in the Fourier series for each temporal mode. The method of ranking each term in the Fourier series is borrowed from the snapshot POD method. If the Fourier coefficient is a representation of how much each Fourier series term contributed to the error between the original signal and the approximation, then the Fourier coefficients are sorted by their magnitude and a Fourier series energy (E b ) can be defined as: In order to reduce the chance of spurious oscillations, the Fourier series energy was chosen to be lower than the eigenvalue energy. The Fourier series representation imposes periodicity and therefore continuity of the POD temporal mode. Spurious oscillations are driven by high frequency components of the Fourier series that appear due to the discontinuity of the POD temporal mode between the last and the first snapshots. The approximation for the inlet variables (Ũ in (x in , t, E λ , E b )) is therefore given by: The reconstructed inlet variable vector is written here as a function of the eigenvalue and Fourier coefficient energies as a reminder that both of these quantities will affect the accuracy of this method. Two parameters that also affect accuracy not considered in this equation are the resolution of the 2D inlet planar mesh and the order of interpolation between the original 3D mesh and the 2D inlet planar mesh that the POD is computed over. As such, α is a scaling parameter used to compensate for inaccuracies in the interpolation process. The inlet planar mesh used here was 200 x 200 points with a resolution of 0.008D. Because the PODFS simulation is run with a timestep twenty times larger than the acoustic simulation, the PODFS was computed from the acoustically forced results using every twentieth time step. The choice of sampling rate is up to the user however the same considerations as are used with constructing Fourier series should be respected. The highest frequency reproduced by the method is determined by the Nyquist criterion.The parameters used in the model reconstruction are available in Table 1. Once the PODFS was computed over the three planes of interest (see Fig. 8) using the data from the acoustically forced simulation an incompressible simulation was run on a domain which only includes the flow downstream of the injector passages (see Fig. 3). The exits of the injector flow passages now become the inlets of this truncated domain and the velocity, k and ω are generated at each timestep using the PODFS models and applied to each of these three inlets.
In the case of this study N p was set such that more than 95% of the fluctuating energy was captured by the POD approximation. As more modes are included, the accuracy of the model improves however the computational cost of generating the model from the donor simulation and the cost of applying the model to the new domain increases. The number of modes chosen here seem sufficient but it is easy to choose fewer or more modes dependent on the problem of interest. The number of POD modes generated is equal to the number of snapshots used to generate the model. As the number of modes used approaches the number of snapshots (N P → N S ) the accuracy of the model approaches a maximum (E λ → 1.00) however any efficiency gains with respect to the strong recycling method disappear. It is therefore up to the user to determine the level of accuracy required. Alpha (α) controls the mean mass flow and was obtained by checking the mass flow rates through each passage in the case of the donor simulation and the PODFS simulation and scaling the mass flow rates accordingly. Figure 4 shows the fluctuation of mass flow rate through the injector passages as a function of time for the acoustically forced flow, mass flow forced case and the PODFS model, the PODFS model clearly outperforms the mass forced case. All mass flows are normalised by the total mean mass flow rate from all passages in the compressible simulation. The time varying total mass flow through the incompressible mass forced simulation has been set equal to that in the compressible simulation. The amplitude of mass flow fluctuation is approximately 5% of the mean mass flow rate. It can immediately be seen that the mass flow amplitude for each passage is different for the acoustic and incompressible mass forced cases. The fluctuations in the incompressible simulation are approximately in line with the bulk mass flow, showing that the mass flow split is constant throughout the forcing. By contrast the amplitude changes between passages in the compressible simulation due to the differing impedances of the passages. For example mass flow fluctuation is significantly reduced in the outermost passage 'C.' This passage to passage variation would have significant effects on the response of a forced flame that simple mass flow forcing is not able to capture. It can also be seen that mean mass flow split is not the same for the two methods. To test whether this is due to the forcing method or simply due to different solvers the incompressible and compressible solvers were run without forcing using the mean mass flow and mean pressure drop respectively. The change in discharge coefficients for the three  Table 2. Very little difference is seen for the incompressible case; further evidence that mass forcing does not affect the mass split. The acoustic forcing on the other hand can be seen to change the mass flow and mass splits. Previous work [22] has shown that the application of acoustic forcing at 300 Hz augments the strength of large scale helical structures downstream of the injector. This effect can be seen in Fig. 5, which shows a phase averaged azimuthal velocity fluctuation through the swirl vanes of passage C for two points in the phase for the acoustically forced simulation and the mass flow forced simulation. Large helical flow structures pass from a region between passages A and B radially outwards, rotating with the swirling flow. These regions of high velocity have an associated lower pressure region which helps accelerate the flow locally leading to variations of axial and azimuthal velocity through the passages in the azimuthal direction. Because the strength of these structures are modified by the forcing the net effect is to modify the mean pressure field and change the mass flow splits between passages. Figure 5 also shows that the velocity fluctuations through the injector are higher in magnitude and have a longer axial wavelength in the case of acoustic forcing. The point of maximum velocity in the azimuthal direction can also be seen to rotate in the direction of the mean flow as the helical structures are convected downstream. Figure 6 shows the fluctuating swirl number for the three passages normalised by the total integrated mean swirl found downstream of the whole injector. It can clearly be seen that the swirl predicted by the acoustic and incompressible mass forced cases are very different at the injector exit planes with higher swirl fluctuations found in the incompressible case. This again would have significant consequences for any FTF calculated using such a forcing  method. With the incompressible solver the axial velocity responds throughout the passage simultaneously (as the speed of sound is infinite) whereas it passes through the compressible simulation as a pulse at the speed of sound. Meanwhile the circumferential fluctuation moves from the vane at the local flow velocity. This leads to a phase shift between axial and circumferential velocity changes, which manifests as a swirl fluctuation, which is different for the two cases. Figures 4 and 6 also show the mass flow and swirl for each passage for the incompressible PODFS based method. It can be seen that this method performs excellently at reproducing the correct mass and swirl inputs into an incompressible simulation. In comparison with the strong recycled method [24] the PODFS method does not require the loading of data at each time step. Instead, all of the data is loaded at the beginning of the simulation and at each time step each processor that has some of the inlet domain must compute the Fourier series that in this case had hundreds of coefficients per mode. In order to test the efficiency of the PODFS method, the same grid was used to run an identical simulation except that the strong recycling method was employed instead of the PODFS method. Despite the need to compute the Fourier series for each mode and at each timestep, the significantly reduced data input of the PODFS method resulted in computations that were only a fraction of the cost of the other methods tested (see Table 3) due to the reduction of the required mesh size which no longer includes swirl vanes. As has been shown above, compressible effects are important in predicting the response of swirl passages, but if the exit velocity from the passages is correct can an incompressible simulation predict the correct downstream flow development? Fig. 7 shows the time history of the total swirl in the duct at the injector exit plane as shown in Fig. 8. The results from the acoustically forced simulation show a complex pattern which is the result of the mixing of the co-and counter-rotating swirl passages and their interaction with the hydrodynamic instability of the flow [22]. This is used as a convenient single measure of the development of the flow. As expected the simple mass flow forced case does not capture the development well. However, the PODFS method is seen to give a much closer representation, particularly of the high frequency fluctuations. The method only reproduces an approximation of the strong recycling method and therefore it is expected that there is a small penalty in accuracy. Figures 4 and 6 show that the method captured the majority of the fluctuations in both mass flow rate and swirl number with a small reduction in the smoothness of the mass flow rate. The PODFS reconstructed axial velocity field for passages A and B at an instant in time is compared to the original data in Fig. 9. It shows that the method can reproduce the interesting flow features such as the vane wakes.

Swirl number and downstream flow development
To further compare the development of the downstream flow the mean axial velocity fields are shown in Fig. 10 for the three simulations. The mass forced and acoustically forced simulations show a similar velocity magnitude at the exit of the injector but as the  Fig. 8 The location of the injector and passage exit planes superimposed over the mean axial velocity from the donor simulation flow develops downstream the mains jet flares outwards towards the wall and the pilot jet opens slightly. One might be tempted to attribute this change to either the numerical method or chosen time step however it can be seen that using exactly the same numerical method, the PODFS simulation reproduces the flow features more accurately suggesting differences in magnitude and phase of the individual passage responses can even change the mean field. This can be explained by the changed mass flow split seen in Fig. 4 and in Table 2. Figures 11 and 12 show the RMS of the resolved velocity fluctuations in the axial and azimuthal directions in the x-y plane for the donor simulation, mass flow forced simulation and the PODFS simulation. The donor simulation includes the full injector geometry, contains 7.5 times the number of cells and is run with a timestep 20 times smaller than the PODFS simulation however the fluctuation velocity fields are almost identical. The fluctuating velocity fields in the case of the mass forced simulation show that this method causes the arms of high axial fluctuation to thicken and flare radially outwards unlike in the other two simulations. As seen in Fig. 5, the azimuthal fluctuations in this case show a short wavelength coherence which leads to a second set of high magnitude fluctuations to form downstream which is not seen in the acoustically forced or PODFS forced simulations. Figures 13 and 14 show the mean and RMS fluctuations of axial velocity from a y-z plane extracted 0.2D downstream of the injector exit plane. It shows that the acoustically and PODFS forced simulations result in a high intensity squaring of the mains jet due to the The reason that an incompressible method with appropriate upstream boundary conditions can accurately capture the effects of acoustic forcing is because the flame in a swirl  stabilised burner is affected primarily by convective and not acoustic effects. The Mach number is also sufficiently low that the mean flow can be considered incompressible (see for example [15]). The convective effects are, in the case of premixed swirl burners, primarily   [1,5]. This vorticity can only be generated in the vicinity of solid boundaries. These effects are much stronger than any localised acoustic effects because the velocities are often an order of magnitude higher than the acoustic particle velocity. The magnitude of the acoustic particle velocity is given as: Whilst the acoustic impedance (Z) of the injector gives the relationship between the fluctuations of the spatially averaged mean axial velocity across the injector exit plane: Where A is the injector exit area. Combining these relations gives the relative strength of these velocity fluctuations: In this study ρc/A|Z| ≈ 4 showing that the average magnitude of axial velocity fluctuations at the injector exit is four times larger than the acoustic particle velocity in the region of the flame. This ratio will also be augmented in the case of a real engine as both the value of ρ and c will be increased at higher pressures and temperatures. The local peak axial velocity fluctuation magnitude may be several times higher, especially in the case were the acoustic forcing increases the strength of large scale helical modes downstream of the injector [22,23]. This analysis shows that the most significant action of the acoustic waves is through their interaction with the fuel injector and once this is taken into account, only convective phenomena are important to the downstream flow.
The accuracy of the PODFS method can be improved by considering a higher proportion of POD modes or a larger number of Fourier coefficients, however increasing E b will often lead to undesirable high frequency oscillations. The POD was also calculated over a set of data with a lower resolution in time than was available from the acoustically forced simulation, increasing this resolution would likely increase the accuracy of the method. The interpolation strategy and order are also significant factors in the accuracy and more sophisticated strategies could improve results.
In this study URANS simulations have been run because of previous work that has shown that URANS is sufficient for capturing the acoustic impedance of this fuel injector [22]. In fully reacting simulations aimed at resolving the flame transfer function, large eddy simulation (LES) is the preferred tool. LES would require a considerably finer mesh in the fuel injector, with a further computational penalty due to the reduced time step, although wall modelled LES would remove the prohibitive requirement to resolve turbulence down to the wall. The PODFS method is also expected to work equally well for LES simulations provided that the mesh resolution is increased to capture the smaller scale structures and the number of modes increased to capture the turbulence that forms inside the injector body. The future aim of this work is to combine PODFS models developed from forced compressible URANS simulations with additional models capable of reproducing the finer scale turbulence such that incompressible and reacting LES calculations can be used to resolve the flame transfer function with a reduced computational cost.

Conclusion
This study highlights the importance of the forcing method used in simulating the FTF of a lean-burn fuel injector where the mass flow and swirl fluctuations of the multiple passages are crucial to correct prediction of flame response. This is especially important in the case of industrially representative fuel injector geometries that are significantly more complicated than typical lab-scale burners. In the case of these industrial designs, the interaction of the flow between each passage needs to be captured accurately. The incompressible mass flow forcing method has been shown to be insufficient to capture these effects compared to a reference acoustically forced compressible simulation validated in previous work leading to significant differences in both mean and RMS of the downstream flow. However, accurate and efficient incompressible simulations are possible by truncating the domain at the exit of each injector flow passage and applying known velocity fields to each inlet extracted from acoustically forced simulations using a novel PODFS method. This model has been introduced and used to overcome several limitations of the strong recycling method with a small reduction in accuracy when compared to the donor simulation. In the future this method will be used to provide inlet boundary conditions to reacting two phase computations such that the FTF of a lean-burn injector can be found in a computationally efficient manner. The PODFS method not only captures the swirl number and mass flow rates better than the mass forced method, it also captures the downstream flow development much more accurately. This shows that provided the flow inside the passages is captured correctly, the downstream flow field will also be correctly predicted. This implies that compressibility effects are only important in the immediate vicinity of the fuel injector. In the future the PODFS method could be used to apply both the fluctuations due to acoustic forcing and large scale turbulent structures to the inlet of two-phase combusting LES simulations. In this case such a simulation over a truncated domain will allow the FTF to be resolved at different engine power settings and different operating conditions without having to resort to expensive compressible LES which would need to include the swirl vane passages and the upstream geometry.