Reduced Order Modeling Framework for Combustor Instabilities Using Truncated Domain Training

Reduced Order Modeling Framework for Combustor Instabilities Using Truncated Domain Training

Jiayang Xu111Graduate Student Research Assistant, Aerospace Engineering Department, AIAA Student Member, Cheng Huang222Assistant Research Scientist, Aerospace Engineering Department, AIAA Fellow, Karthik Duraisamy333Associate Professor, Aerospace Engineering Department, AIAA Fellow,
University of Michigan, Ann Arbor, MI, 48109

A multi-fidelity framework is established and demonstrated for prediction of combustion instabilities in rocket engines. The major idea is to adapt appropriate fidelity modeling approaches for different components in a rocket engine to ensure accurate and efficient predictions. Specifically, the proposed framework integrates projection-based Reduced Order Models (ROMs) that are developed using bases generated on truncated domain simulations. The ROM training is performed on truncated domains, and thus does not require full order model solutions on the full rocket geometry, thus demonstrating the potential to greatly reduce training cost. Geometry-specific training is replaced by the response generated by perturbing the characteristics at the boundary of the truncated domain. This training method is shown to enhance predictive capabilities and robustness of the resulting ROMs, including at conditions outside the training range. Numerical tests are conducted on a quasi-1D model of a single-element rocket combustor and the present framework is compared to traditional ROM development approaches.

Reduced Order Modeling Framework for Combustor Instabilities Using Truncated Domain Training


 = amplification factor
 = cross sectional area
 = stoichiometric fuel-to-oxidizer ratio
 = cumulative energy in singular values
 = pressure oscillation growth rate
 = characteristic invariant for incoming acoustic wave
 = number of modes in proper orthogonal decomposition (POD) basis
 = chamber length
 = starting location of fuel injection
 = ending location of fuel injection
 = fuel mass flow rate
 = fuel consumption rate
 = oxidizer consumption rate
 = unsteady heat release rate
 = singular value
 = time lag
 = POD basis
 = complementary basis of
 = variable for fuel injection shape description
 = oxidizer mass fraction

I Introduction

Combustion instabilities have long been a hurdle in the development of modern rocket engines. These instabilities are characterized by the coupling between acoustics, hydrodynamics and heat release. In propulsion systems, the triggering of combustion instabilities can lead to catastrophic engine failures, and the underlying physical mechanisms are sensitive to many parameters including the geometry and fuel-to-oxidizer ratio. Therefore, a significant number of geometric configurations and operating conditions have to be evaluated and analyzed in rocket engine design. Though the advancement of modern computational technology has enabled routine modeling of laboratory-scale rockets  [1] and other experimental configurations such as atmospheric combustors [2], direct numerical simulations (DNS) [3] and large eddy simulations (LES) [4, 5, 1, 6, 7, 8] continue to be prohibitively expensive for the simulation of full-scale engines [9, 10, 11]. This motivates research on alternative approaches such as reduced-fidelity modeling [12, 13, 14, 15, 16] and reduced order modeling [17, 18].

Reduced-fidelity modeling approaches usually adapt certain levels of simplifications in physical and numerical models to achieve higher computational efficiency. Sirignano and Popov [14] developed a two-dimensional model for transverse-mode combustion instability in a cylindrical rocket motor and further extended it to rectangular configurations [15]. Though reduced-fidelity modeling approaches can be efficiently used for engine design with satisfying accuracy for critical quantities of interest, many challenges are encountered in accurately predicting crucial combustion dynamics that can trigger instability.

In addition to approaches that model combustion dynamics by solving transport equations, empirical approaches have been pursued. These approaches formulate the combustion response as a function of well-characterized flow variables such as pressure and velocity. This class of techniques include the so-called flame transfer (FTF) [19] and flame describing functions (FDF) [20]. In Ref. [16], more factors such as variations in geometry are taken into account and the response function is successfully applied to a swirl-stabilized combustor. Though FTF/FDF approaches impose a simple model form for the flame dynamics that can be easily derived from both experimental and computational data, their applications are usually restricted to combustion system with low-amplitude perturbations with regards to the mean flow due to the fact that the derivation of the flame model is based on linear assumptions in frequency domain. Therefore, to model highly nonlinear flame dynamics which is typical in rocket engines, the validity of the FTF/FDF approach remains to be evaluated.

In contrast, reduced order modeling (ROM) approaches pursue mathematically formal reductions, offering the promise of generalization. Projection-based ROMs [21, 22] have recently been applied to combustion instability problems [23, 24, 25, 26, 27]. In projection-based ROMs, the full order model is projected onto a set of informative basis functions using proper orthogonal decomposition (POD)-Galerkin method [28], resulting in a low dimensional set of equations that can retain most of the modeling fidelity. Thus, the need for commonly used empirical modeling strategies such as combustion response functions can be eliminated. Preliminary explorations of projection-based ROMs on representative combustion model problems can be found in Refs. [23, 25, 26, 27]. It should be pointed out that POD-based ROMs have been more widely used in non-reacting flow problems such as aerodynamics [29, 28], aeroelasticity [30, 31] and flow control [32, 33].

However, applications of ROMs to rocket engines involve the following challenges:

  • A single high fidelity simulation or Full Order Model (FOM) of a full engine may be prohibitively expensive.

  • Even if full engine FOM simulations are affordable, conventional ROM construction is based on FOM simulations for different injector configurations. If there is any change to the geometry, the ROM has to be re-generated or re-parameterized, which restricts the scope of ROM applications.

  • POD-ROMs may be incapable of providing predictions beyond the training interval and parameters.

  • Since POD basis construction involves global minimization, the resulting basis will be dominated by high-energy regions. The dynamics in low-magnitude regions can be under-resolved and consequently affect stability and robustness.

The present work extends the multi-fidelity framework methodology developed by Huang et al. [26, 27, 34] and assesses ROM performance on the major aspects listed above based on a quasi-1D model combustor problem, which was originally developed by Smith et al. [35] to model the experimental work by Yu et al. on the single-element combustor [36, 37]. It should be mentioned that additional challenges exist in ROMs of more realistic problems, where complex physics such as three dimensional turbulent combustion, vortex shedding, and flame-vortex interactions are present. While progress is being made in addressing these issues [25, 38, 39], the present work focuses on the issue of ROM development in problems where the full order solution is not available for the entire domain. The approach itself is generic for problems of different levels of complexity. The framework pursues a multi-domain approach where the computational domain is divided into two components: the upstream one covering the critical physics in the heat addition/acoustic-flame interaction region and the downstream part dominated by flow and acoustics dynamics. Reduced order models are used to model the critical dynamics in the upstream component. The ROM is developed and trained based on FOM solutions of the truncated domain corresponding to the upstream component with a perturbed characteristic downstream boundary condition. This training strategy is designed to incorporate rich information in the resulting ROM and proves to be effective in enhancing predictive capabilities. The truncated domain ROM is then integrated within a multi-fidelity solver 111For the current demonstration, we replace the reduced-fidelity part with the full order model for consistency of evaluations. The multi-fidelity model is then compared to ROMs based on the FOM solutions of the entire computational domain, which will be referred as conventional ROMs in the current paper.

The rest of this paper is organized as follows: The quasi-1D combustor setup and formulations for the FOM and ROM are described in Sec. II. In Sec. III, the present framework and ROM training with characteristic boundary are introduced. Numerical experiments are presented in Sec. IV. Concluding remarks are given in Sec. V

Ii Formulation

ii.a Full order model

In this work, a quasi-1D version of the single element Continuously Variable Resonance Combustor (CVRC) [36, 37] is used as the full order model. The geometry of the CVRC is sketched in Fig. 1. Numerical studies are performed with multiple chamber lengths, which exhibit various instability behaviors. The geometric parameters for other sections are given in Table 2. To avoid invalidating the quasi-1D assumption, the back-step and the converging part of the nozzle are sinusoidally contoured in this study, in contrast to a discontinuous change in the cross-sectional area in the experimental setup. The same gas properties and operating conditions as in Ref. [23] are adopted and listed in Table 3.

Section Injector Back-step Nozzle converging part Nozzle diverging part
Length (m) 0.1397 0.0064 0.0127 0.034
Radius (m) 0.0102 0.0102 to 0.0225 0.0225 to 0.0104 0.0104 to 0.0195
Table 2: Geometry parameters
Parameter Value
Fuel mass flow rate, kg/s 0.027
Fuel temperature, K 300
Oxidizer mass composition
Oxidizer mass flow rate, kg/s 0.32
Oxidizer temperature, K 1030
Fuel composition
Equivalence ratio 0.8
Table 3: CVRC operating conditions

The governing equations in this problem are the quasi-1D unsteady Euler equations with species transport, represented as




In the conservative variable vector , is the density, is the velocity, is the total internal energy, is the oxidizer mass fraction, and is the cross sectional area. The corresponding convective fluxes are represented by the vector , where is the static pressure.

In the three source terms, is due to area variations and the latter two terms represent combustion. In the experiment [37], the fuel is injected through an annular ring located at the back-step and reacts at a finite rate with the oxidizer injected. In this work, we follow the choice of Frezzotti et al. [40, 41] and assume an infinitely-fast one-step combustion model. An important assumption behind the model is that when fuel is injected, it will react with the oxidizer instantaneously to form products and no intermediate species are produced, thus only one species transport equation is needed. This process is accounted for in . As suggested by Frezzotti et al. [40, 41], to reproduce the combustion region of a finite width and avoid discontinuities, the fuel injection process is modeled in a sinusoidal shape at the rate of , yielding


where is the total mass flow rate of fuel injection, and and are the starting and ending location of the flame, respectively. Due to the infinitely-fast model, the consumption rate of the oxidizer is


where is the stoichiometric fuel-to-oxidizer ratio. The shape of the fuel injection and resulting flame given by Eq. (3) is also shown in Fig. 1 along with the CVRC geometry.

Fig. 1: Computational domain and flame shape

The final source term models the unsteady heat release and represents the response function introduced by Crocco et al. [42], given in Eq. (5). The function takes into account the coupling between acoustics and combustion by expressing the unsteady part of the heat release as a function of pressure with an amplification parameter and a time lag .


In the simulation, a steady state is first achieved with the unsteady source term turned off. From the steady state, a low-amplitude perturbation is applied to the inlet boundary to trigger the instability at the beginning of predictive unsteady simulations (either FOM or ROM). After the perturbation, the unstable response will become self-excited and start to grow under certain parameters. The instability behavior is described in Sec. IV.

Numerical solutions of Eq. (1) will be used as a surrogate for the high-fidelity solution in the current work, and the corresponding simulation is referred as FOM.

ii.b Projection-based ROM

To derive the projection-based ROM equations, we first rewrite the full order model given by Eq. (1) for each variable as


where the variable index corresponds to the four conserved variables respectively. The residual vector for variable consists of the corresponding flux and source terms


For each variable, defining an individual POD projection basis spanning a subspace , and a complementary basis spanning , such that , the following decomposition can be derived:


The two bases and are obtained by performing singular value decomposition (SVD) on the snapshots of the -th conserved variable. Defining to be the designed size of the reduced order system for variable , is constructed column-wise from the first left-singular vectors from the SVD 222The remaining left-singular vectors are . The conserved variables are then approximated as


where is the -th element in the reduced order variable vector , and is the -th column in the POD basis . Since each left-singular vector represents a spatial mode of the FOM solution, the approximation can be viewed as a superimposition of the first leading modes, with being the temporal coefficient in the evolution of the -th mode.

The Galerkin-projected ROM for coefficients is then


For clarity of presentation, is used to denote the approximated conserved variables, including .

In practice, to achieve further reduction in computational cost and to improve robustness, additional constructs such as sparse sampling techniques and stabilization will be required. While we are developing such techniques [23, 24, 25, 43], they will not be considered in the present work, as the focus is on multi-domain modeling.

Iii Framework Overview

In this section, we explore characteristic-based ROM training on a truncated domain and contrast it with conventional ROM development. A schematic is given in Fig. 2. The major steps include

  1. Truncate the domain into two sub-domains, the first containing the physics-complex area that covers the fuel injection and the flame, and the second containing the variable-length chamber dominated by acoustics. The ROM will be trained for the first sub-domain.

  2. Perform a FOM training simulation on the first sub-domain with a broadband perturbation added on the truncation interface, which is treated as a characteristic boundary.

  3. Generate the POD bases for different variables using SVD of the solution from the training simulation.

  4. Use the POD bases in a ROM solver for the first sub-domain and couple it with a FOM solver for the variable-length, acoustics-dominated chamber for predictions.

It can be noted that in the framework ROM is only developed for the first sub-domain, whereas the second sub-domain is solved in FOM. There are two reasons for this choice:

  1. The complex, computation-intensive area is fully covered in the first sub-domain. The accurate modeling of this domain requires high-resolution simulations. The FOM computation in the second sub-domain is expected to be much less challenging than the first domain and affordable with coarser resolution modeling approaches (e.g. coarse-mesh LES and URANS), which makes a ROM replacement unnecessary.

  2. In design evaluations, the chamber length in the second sub-domain is variable, thus a new ROM training is required for each chamber length, which violates the goal of reducing the computing cost using ROM.

Fig. 2: Schematic of framework

iii.a Characteristic ROM training on reduced-domain

In traditional ROM formulations, POD bases are generated using FOM solutions of the target full domain. In the proposed approach, the domain is split into two sub-domains with the interface at m. The first sub-domain includes the injector, back-step and the leading part of the chamber where the flame is located, and the second includes the rest of the geometry, which is variable in design evaluations.

To obtain a basis that is representative of the physics in the first sub-domain, a FOM simulation is performed in a reduced domain, treating the truncation interface as a characteristic outlet boundary, which helps to eliminate the undesirable resonant acoustic modes corresponding to the reduced-domain geometry. As demonstrated by Huang et al. [27], the presence of such resonant acoustic modes can significantly affect the predictive capabilities of the ROMs. A similar training procedure as in Refs.  [26, 27] is used. The training FOM simulation is also started from the steady state as in the self-excited simulations on the full domain (Sec. II). The difference is that instead of a single frequency inlet perturbation, a broadband perturbation is imposed on the incoming characteristic wave at the truncation boundary to cover the range of responses created by the resonant frequencies corresponding to different chamber lengths. This perturbation is imposed over the entire duration of the FOM simulation.

To maintain consistency, the rest of the specified properties at the boundary are as follows:


where is the speed of sound, is one-dimensional approximation of the characteristic invariant for the in-coming acoustic wave, is obtained from the steady state solution, and the other three primitive variables are extrapolated from the values of the interior cells.

iii.b Multi-domain coupling

Following the aforementioned simulations, the bases for the conserved variables are obtained from the reduced domain solution snapshots, and used in POD-Galerkin projection (Eq. (10)) to derive the ROM for the first sub-domain. The second sub-domain is computed using the full order equations (1). Due to the hybrid integration, chamber lengths can be changed without any further modifications to the ROM. The two sub-domains communicate at every time-step by exchanging interface conditions to be the value at the neighboring cell of the adjacent domain. Similar to the interior of the computing domain, Roe’s upwind flux [44] is used at the interface. Using subscript I and II for the first and second sub-domain, and superscript and for the first and last cell in the corresponding sub-domain respectively, the flux at the interface is given by


where is the Roe-averaged state calculated from and . is computed from the last row of the basis and the reduced variable, .

iii.c Control groups for comparison

The methodology detailed above is compared to the following control groups using the conventional ROM approach:

Control group A

also uses a hybrid multi-domain solver, i.e. the first sub-domain is solved using the ROM and the second solved using the FOM. The difference from the proposed framework is in the training data generation and collection stage. In control group A, for each combination of chamber length and , a full-domain FOM simulation is performed instead of the proposed characteristic training on the truncated-domain. Then the solution is restricted to the first sub-domain and collected to generate the POD bases for the ROM of the sub-domain.

Control group B

uses a traditional full-domain ROM. The same FOM training simulations as in control group A are used and the POD bases are directly generated on the full-domain solution.

To summarize, let and be the number of chamber lengths and amplification factors studied, respectively. In the proposed framework, FOMs simulated on the reduced domain are used in total 333although as shown in the next section, this can be reduced using characteristic perturbations. Then the first sub-domain is simulated using the ROM, second sub-domain using the FOM. In control groups A&B, self-excited FOM simulations are conducted. In group A, the first sub-domain is solved using a ROM, and the second is solved using the FOM. In group B, the whole domain is solved using a ROM.

For conciseness, the proposed framework will be referred to as “reduced-domain training and multi-domain solver (RD-MD)", control group A as “full-domain training and multi-domain solver (FD-MD)", and control group B as “full-domain training and full-domain solver (FD-FD)". A schematic of the different methods is given in Fig. 3.

Fig. 3: Schematic of training (outlined in black) and prediction (outlined in blue) stages, curves indicating perturbation added at the boundary

Iv Results

Numerical tests are conducted for chamber lengths ranging from 0.254 m to 0.762 m with an interval of m and at and . In the proposed training method, a FOM simulation is performed in reduced domain for each . From a priori [45] analysis, the longitudinal frequency of the chamber spans approximately between 700 and 2000 Hz, and thus the broadband perturbation signal is set to


which is a superimposition of frequencies, starting from with an incremental interval .

To train the control groups, full-domain self-excited FOM simulations are conducted for each combination of and . In all three methods, snapshots of the training solution are collected every 100 time-steps over a period s, with the initial steady state set at . Once the basis is obtained using the SVD of the corresponding snapshot, ROM simulations are conducted over s.

iv.a FOM results

Results from the full-domain FOM simulations are presented to characterize the instability behavior. The predicted response is visualized using pressure signals obtained 0.0127 meters upstream of the converging part of the nozzle, which is a typical location selected to probe combustion instabilities in the work conducted in Ref. [40].

With different combinations of parameters, two general categories of responses are identifiable: one with positive growth rate, in which the pressure oscillation grows and settles into a limit cycle oscillation (LCO), and one with negative growth rate in which the instability starts to decay after the perturbation ends and the flow converges to a steady state. The definition of the growth rate, , is based on the peak-to-peak amplitude of the unsteady part of the pressure signal, whose growth is fitted to a exponential function plotted in Fig. 4, given by


To better distinguish the categories, two representative examples (at m for the growing category and m for the decaying category, respectively) are presented in Fig. 5 and 6, where the difference in growth rate can be clearly observed via the monitored pressure signals and the spatio-temporal diagrams of the pressure evolution. Similar responses can be found in previous studies [23, 46]. To provide more insight into the dynamics and the coupling of pressure and heat release described by Eq. (5), the pressure distribution and the unsteady heat release term are also plotted in the same figures. Each figure includes two snapshots at the beginning of the simulation, and two snapshots towards the end. It can be seen that in the case showing pressure amplitude growth (typically referred to as unstable case), the unsteady heat release has grown significantly with the pressure oscillation at the end of the simulation compared with the beginning due to the fact that the oscillations of pressure and heat release rate are in-phase, which is considered favorable to drive combustion instability. While for the other case with pressure amplitude decay (usually referred to as stable case), for some time instances, the pressure and heat release rate oscillations are out-of-phase, which is recognized as a combustion instability damping mechanism.


Fig. 4: Definition of growth rate
(a) Pressure signal at the monitored location
(b) Spatio-temporal diagram of the pressure evolution
(d) Pressure and unsteady heat release profiles. Dotted line: steady state pressure, dashed line: pressure signal monitor location
Fig. 5: Example response with positive growth rate, m.
(a) Pressure signal at the monitored location
(b) Spatio-temporal diagram of the pressure evolution
(d) Pressure and unsteady heat release profiles. Dotted line: steady state pressure, dashed line: pressure signal monitor location
Fig. 6: Example response with negative growth rate, m.

iv.b Modal decomposition

iv.b.1 Singular values

The singular values from a POD of snapshots of density, from different training methods are compared at m. There are 472 cells in the first sub-domain, 897 in the full domain at m, 1405 at 0.508 m and 1913 at 0.762 m. Since the number of snapshots is larger than the mesh size, the maximum number of modes and singular values from SVD is limited to the mesh size.

The result is visualized by the complementary part of the cumulative energy in Fig. 7. For modes, the cumulative energy is defined as


It should be noted that there is only one curve for RD-MD as it uses the same training simulation and POD bases for all chamber lengths. It can be observed that for both FD-MD and FD-FD, the decay in the complementary part increases as the chamber length decreases. This is due to the fact that the dynamics have a higher frequency when the chamber length is shorter, which is easier to be captured by fewer modes in SVD. In all cases, the decay in FD-FD is slower than FD-MD, because it covers a larger domain with more spatial variations in geometry and physics, and more modes are therefore needed to contain the same portion of information stored in the training data. For medium-to-high , the proposed framework gives the best decay since its training data contains more high-frequency contents. At small , the decay in the leading modes for the two full-domain training methods become slightly better than RD-MD as the response frequency in their training data is comparable to the highest one included in the reduced-domain training, whereas there are also lower frequencies contained in the latter. Results for other variables and follow a similar pattern.

Fig. 7: Complementary part of cumulative energy in the singular values for at

iv.b.2 Mode shapes

Comparison of just the singular values does not provide enough information on the compatibility of the basis and physical features, since each training method is different. To provide further insight, the first five most energetic modes for at m, are shown in Fig. 8. It should be noted that the basis for FD-FD is on the full domain and only truncated in the plot to the first sub-domain for comparison purposes. The first mode for all methods appears to be identical. For lower energy modes, the basis for the RD-MD is quite different from the other two methods. The basis for the FD-MD and FD-FD are also different for , which suggests that performing POD on different domains can result in different mode shapes, even when the training data is collected from the same simulation. This is because that the SVD process for POD bases generation solves a global minimization problem that will yield a basis with minimal least square projection error for the input snapshots. Thus, the scope of the minimization problem is changed with the size of the domain. There is hence no guarantee that the same local mode shapes will be obtained, though the input snapshots are the same locally. Consistent observation has been made in Ref. [47], where the mode shape in full domain and different sub-domains of a 3D simulation of the CVRC is compared (equivalent to the FD-MD and FD-FD comparison).

To better demonstrate the evolution of the modes, the FOM solution at m is projected onto the POD bases. Due to orthogonality, the projection can be performed independently and the resulting reduced variables are given by


where is the -th POD mode.

The coefficients for the five leading modes present in Fig. 8 are shown in Fig. 9. For better visualization, only a period of s is used. It can be seen that the first mode evolves at a single frequency corresponding to the large-scale instability characterized by the pressure oscillation described in the previous section. This frequency is exhibited by the latter modes, while there are also additional frequencies corresponding to the high-frequency dynamics contained in the modes.

It can be noticed from Fig. 8 that the spatial modes that exhibit high-frequency temporal oscillations also have a higher spatial frequency. To make this trend clearer, the dominant spatial frequencies of the 50 leading modes of are plotted in Fig. 10. While showing a similar trend, the growth rates are different between the methods. The reasons are two-fold: 1) using the same single-frequency training, FD-MD shows a faster decay in the singular values, which indicates that the low-frequency high-energy dynamics is better captured in the leading modes, and thus the high-frequency modes are realized at smaller mode numbers than in FD-FD. 2) There are multiple frequencies in the reduced-domain training. These frequencies are still significantly higher than those of the corresponding POD mode numbers in FD-FD and FD-MD. The corresponding modes representing these frequencies cause the RD-MD to have the slowest increase in the dominant spatial frequencies. Also presented in Fig. 10 are the dominant temporal frequencies exhibited in the evolution of the POD modes. For RD-MD, the frequency increases with mode number steadily before the 35-th mode. This demonstrate a good correlation between the temporal and spatial frequency captured in the modes. In contrast, many high-spatial-frequency modes from the two full-domain training methods actually evolve at a low temporal frequency. This inconsistency indicates a less efficient basis compared to RD-MD.

(a) RD-MD
(b) FD-MD
(c) FD-FD
Fig. 8: Leading spatial modes for at m
(a) RD-MD
(b) FD-MD
(c) FD-FD
Fig. 9: Evolution of coefficients at m
(a) Spatial
(b) Temporal
Fig. 10: Dominant spatial/temporal frequency of POD modes for at m.
Solid line: first mode, shaded: band for the first five frequencies

iv.b.3 Projection error

For a given set of POD modes, the projection error is defined as


where is the FOM solution snapshot matrix, is the first columns of the left singular vector from the SVD of training data. Again, at m, the projection error of up to is shown in Fig. 11. It can be observed that both RD-MD and FD-MD provide a monotonic decrease in the projection error. However for FD-FD, the error increases with at a few cases, which is not desirable for ROM development as it makes the choice of basis size more uncertain. The increase can be again attributed to the fact that the POD method solves a global minimization problem and the optimal bases for the full domain is not necessarily optimal for the reduced domain. The improved projection error property with the multi-domain method implies a higher accuracy and stability of ROM, which is confirmed in the following sections.

The above analysis are focused on and results for other variables follow the same trend.

(a) m
(b) m
(c) m
Fig. 11: Projection error in for

iv.c ROM evaluations

Before integrating the ROM into the Multi-fidelity Framework, the predictive capability of the ROM obtained from the reduced-domain simulation is evaluated. In this test, the ROM trained with broadband characteristic perturbation is tested by feeding single frequency perturbations on the characteristic boundary. This is to ensure that the ROM is able to predict the essential dynamics when the downstream section is dominated by single frequency acoustics. Aside from the perturbation frequency, there is no difference between the conditions used in the test and the truncated domain training (Sec. III).

The tests are conducted over s at five frequencies between 700 and 2000 Hz. It should be noted that some of the frequencies are not included in the broadband training. The number of ROM modes and amplification factor are , respectively. The results are compared with the reference FOMs using the same conditions, and the L2 error over the simulated period as well as at the last time step is plotted in Fig. 12. It can be seen that at all perturbation frequencies, the ROM is able to accurately predict the FOM behavior with an L2 error below . These evaluations confirm the fact that the ROM trained with broadband perturbation is capable of predicting dynamics over a wide range of frequencies.

Fig. 12: Sanity test for L2 error

iv.d Multi-fidelity framework evaluations

The solution from the three different techniques are compared against the FOM over the same parameters and time period. From Fig. 11 it can be seen that the projection error for RD-MD at is comparable to that for FD-FD at . The three methods are compared at the two basis sizes, in this section. Comparisons of the relative L2 error of the ROM solution are shown in Fig. 13. It should be mentioned that in several cases, the ROMs in the control groups are numerically unstable. The error is set to 1 under these circumstances, which is higher than the error in all stable cases.

When the basis size is sufficiently large, i.e. at , RD-MD performs similarly to FD-MD, except for one case at m, where the latter is unstable. At high , where the instability is stronger, the advantage of RD-MD becomes more significant. This can be observed from the medium cases at , where the gap between RD-MD and FD-MD is apparent. Also, FD-FD performs worse as increases, resulting in higher errors and more unstable cases. It should be noted that FD-FD sometimes outperforms the other two methods at low , which is consistent with an improved decay in the singular values. However at m, it becomes unstable.

When the basis size is reduced to , all the predictions deteriorate considerably. The two full-domain training methods become unstable in most conditions. In contrast, RD-MD remains stable at all combinations of parameters, although the ROM solution error is approximately one order higher than that of the cases. The advantage in stability confirms the conclusion that the reduced-domain training results in a better set of basis functions. By virtue of using a smaller basis set, the computational cost, e.g. the projection computation in Eq. (10), decreases linearly when an explicit time-marching scheme is used. In contrast to the traditional full-domain training methods, for which the stability is less predictable, the proposed framework provides a broader stability envelope, and consequentially a more flexible balance between ROM accuracy and efficiency.

Fig. 13: Relative L2 error of ROM solutions, numerically unstable cases set to 1

Major quantities of interest in rocket combustor design include the dominant acoustic frequencies, the growth rates and LCO peak-to-peak amplitudes of the pressure oscillations. It should be noted that the dominant acoustic frequencies have been well-predicted by all three methods. Therefore only the comparisons of the other two quantities, growth rates and LCO amplitudes, are shown in the current study, as can be seen in Figs. 14 and 15 respectively for all the methods. The relative performance between different methods follows a similar relation as in the L2 error analysis. It can be seen that the proposed RD-MD framework is able to predict the relation between the growth rate and accurately with an error below at and below at in most cases, which illustrates its effectiveness.

Fig. 14: Growth rate of pressure oscillation 0.0127 m upstream of nozzle
Fig. 15: LCO peak-to-peak amplitude of pressure oscillation 0.0127 m upstream of nozzle

Finally, to provide a more direct comparison between the ROM solution from the framework and the FOM solution, spatial pressure profiles for m, are provided at several time instances in Fig. 16. These plots are focused on the back-step area for better visualization, and are selected by the end of the simulation, when the error is maximum. It is observed that the RD-MD solution correlates well with the FOM, and connects smoothly across the interface.

(a) s
(b) s
(c) s
Fig. 16: Spatial profiles of pressure

iv.e Off-design condition performance

To further assess the RD-MD framework, the ROM trained using is evaluated at several off-design conditions. The evaluations include two cases at chamber lengths m and 1.016 m. These cases are characterized by dominant acoustic frequencies 2550 and 500 Hz, respectively, which are outside the range of the training frequencies given by Eq. (13). The other conditions have amplification factors that deviate significantly from the training simulation, including . The results for the L2 error (), frequency (), growth rate (), LCO peak-to-peak amplitude () and their relative errors are summarized in Table 4. It should be noted that due to the characteristic training method, the training case does not have a growing response as in the predictions, thus no growth rate or LCO amplitude is reported for the training set. Moreover, the results are plotted along with the designed conditions in Fig. 17 and 18 for a better illustration of their relation.

It is observed that, at all the off-design conditions, the relative L2 errors are below . For the case with shorter and higher instability frequency (OD ), and the four cases with deviated (OD to OD ), the ROM performance is comparable to that in the designed conditions. However for the case with longer and lower instability frequency (OD ), the error in LCO amplitude is , which is more than 5 times higher than the other cases. The result indicates that for this problem, when operating within the frequency range of the training perturbation, the RD-MD method is not only independent of the chamber lengths, but also insensitive to changes in the unsteady heat release term. When beyond the training range, the instability frequency influences the predictive capabilities of the framework, and demonstrates the importance of the multi-frequency perturbation in the characteristic training.

Case (m) (Hz) () (Mpa)
Training N/A 3.25 N/A 700 to 2100 N/A N/A N/A N/A N/A
OD 0.1905 3.25 2516/2516 0 -139.98/-139.79 0/0 0
OD 1.016 3.25 490/490 0 12.91/12.84 0.4576/0.4688
OD 0.508 2.85 1000/1000 0 28.87/28.87 0.0522/0.0521
OD 0.508 3.05 1000/1000 0 61.99/62.00 0.1545/0.1543
OD 0.508 3.45 980/980 0 122.38/122.41 0.6366/0.6366
OD 0.508 3.65 980/980 0 130.19/130.27 0.8368/0.8368
Table 4: Off-design condition results ( are listed in a “FOM/ROM" style, all errors are relative)
Fig. 17: Results at off-design
Fig. 18: Results at off-design

V Conclusions

In this work, a multi-fidelity ROM development framework is investigated on a quasi-1D version of the CVRC model rocket combustor. An model that couples the pressure oscillation and heat release is used to represent the unstable behavior of the combustor. In the proposed framework, the domain is split into two sub-domains with the first sub-domain containing the main reaction dynamics and the second covering the adjustable chamber of the CVRC. Characteristic training is conducted in the first sub-domain, resulting in a ROM that can be directly integrated with a FOM solver of the second sub-domain with any length.

Numerical tests are conducted at different chamber lengths and amplification factors . The framework is compared against the FOM solution, and other ROM approaches where traditional training approaches are taken. A major advantage of the current framework is that it significantly reduces the number of FOM simulations required to train ROMs while the traditional method requires a separate FOM simulation for each individual target chamber lengths, which does not fit in the needs for more efficient rocket engine design. Moreover, the proposed method shows faster decay of the singular value spectrum at medium-to-high chamber lengths, which implies better basis quality and ROM reliability at a low number of modes. The advantage is especially distinct in the tests with , where the proposed method shows significantly improved stability.

In predictive tests at conditions outside the training set, the framework showed significant improvement in both stability and accuracy over the traditional methods, especially when the instability is more pronounced at high , and when number of modes is low. Comparison of L2 error, growth rate and LCO peak-to-peak amplitude shows that the framework is able to predict the these quantities accurately at all combinations of and at .

In summary, the multi-fidelity framework proves to be promising approach for modeling rocket combustion instability. The modularity of the framework will be useful in efficient training and integration of components such as injector elements. The flexibility of the multi-domain approach also enables potential developments such as the application of different time-marching schemes, basis sizes and acceleration methods in different parts of the computational domain to achieve a better balance between efficiency and reliability in the future ROM development.

While the results are encouraging, it should be recognized that the flow and combustion models used in this study are highly simplified and further studies are required on more complex flow problems to evaluate the capabilities of the multi-fidelity framework.

Vi Acknowledgments

The authors acknowledge support from the Air Force under the Center of Excellence grant FA9550-17-1-0195, titled “Multi-Fidelity Modeling of Rocket Combustor Dynamics.”


  • Harvazinski et al. [2015] Harvazinski, M. E., Huang, C., Sankaran, V., Feldman, T. W., Anderson, W. E., Merkle, C. L., and Talley, D. G., “Coupling between hydrodynamics, acoustics, and heat release in a self-excited unstable combustor,” Physics of Fluids (1994-present), Vol. 27, No. 4, 2015, p. 045102.
  • Matsuyama et al. [2016] Matsuyama, S., Hori, D., Shimizu, T., Tachibana, S., Yoshida, S., and Mizobuchi, Y., “Large-Eddy Simulation of High-Frequency Combustion Instability in a Single-Element Atmospheric Combustor,” Journal of Propulsion and Power, 2016, pp. 628–645.
  • Domingo et al. [2005] Domingo, P., Vervisch, L., and Réveillon, J., “DNS analysis of partially premixed combustion in spray and gaseous turbulent flame-bases stabilized in hot air,” Combustion and Flame, Vol. 140, No. 3, 2005, pp. 172–195.
  • Ihme et al. [2012] Ihme, M., Zhang, J., He, G., and Dally, B., “Large-eddy simulation of a jet-in-hot-coflow burner operating in the oxygen-diluted combustion regime,” Flow, turbulence and combustion, 2012, pp. 1–16.
  • Huang et al. [2014] Huang, C., Yoon, C., Gejji, R. M., Anderson, W., and Sankaran, V., “Computational study of combustion dynamics in a single-element lean direct injection gas turbine combustor,” 52nd Aerospace Sciences Meeting, 2014, p. 0620.
  • Lacaze et al. [2009] Lacaze, G., Cuenot, B., Poinsot, T., and Oschwald, M., “Large eddy simulation of laser ignition and compressible reacting flow in a rocket-like configuration,” Combustion and Flame, Vol. 156, No. 6, 2009, pp. 1166–1180.
  • Hernández-Pérez et al. [2011] Hernández-Pérez, F., Yuen, F., Groth, C., and Gülder, Ö., “LES of a laboratory-scale turbulent premixed Bunsen flame using FSD, PCM-FPI and thickened flame models,” Proceedings of the Combustion Institute, Vol. 33, No. 1, 2011, pp. 1365–1371.
  • Srinivasan et al. [2015] Srinivasan, S., Ranjan, R., and Menon, S., “Flame dynamics during combustion instability in a high-pressure, shear-coaxial injector combustor,” Flow, Turbulence and Combustion, Vol. 94, No. 1, 2015, pp. 237–262.
  • Urbano et al. [2016] Urbano, A., Selle, L., Staffelbach, G., Cuenot, B., Schmitt, T., Ducruix, S., and Candel, S., “Exploration of combustion instability triggering using large eddy simulation of a multiple injector liquid rocket engine,” Combustion and Flame, Vol. 169, 2016, pp. 129–140.
  • Staffelbach et al. [2009] Staffelbach, G., Gicquel, L., Boudier, G., and Poinsot, T., “Large Eddy Simulation of self excited azimuthal modes in annular combustors,” Proceedings of the Combustion Institute, Vol. 32, No. 2, 2009, pp. 2909–2916.
  • Wolf et al. [2012] Wolf, P., Balakrishnan, R., Staffelbach, G., Gicquel, L. Y., and Poinsot, T., “Using LES to study reacting flows and instabilities in annular combustion chambers,” Flow, turbulence and combustion, Vol. 88, No. 1-2, 2012, pp. 191–206.
  • Baukal Jr et al. [2000] Baukal Jr, C. E., Gershtein, V., and Li, X. J., Computational fluid dynamics in industrial combustion, CRC press, 2000.
  • Frezzotti et al. [2014] Frezzotti, M. L., Terracciano, A., Nasuti, F., Hester, S., and Anderson, W. E., “Low-order model studies of combustion instabilities in a DVRC combustor,” 50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference, 2014, p. 3485.
  • Sirignano and Popov [2013] Sirignano, W. A., and Popov, P. P., “Two-dimensional model for liquid-rocket transverse combustion instability,” AIAA Journal, Vol. 51, No. 12, 2013, pp. 2919–2934.
  • Popov and Sirignano [2016] Popov, P. P., and Sirignano, W. A., “Transverse combustion instability in a rectangular rocket motor,” Journal of Propulsion and Power, Vol. 32, No. 1, 2016, pp. 620–627.
  • You et al. [2005] You, D., Huang, Y., and Yang, V., “A generalized model of acoustic response of turbulent premixed flame and its application to gas-turbine combustion instability analysis,” Combustion Science and Technology, Vol. 177, No. 5-6, 2005, pp. 1109–1150.
  • Lieu et al. [2006] Lieu, T., Farhat, C., and Lesoinne, M., “Reduced-order fluid/structure modeling of a complete aircraft configuration,” Computer methods in applied mechanics and engineering, Vol. 195, No. 41-43, 2006, pp. 5730–5742.
  • Lucia et al. [2004] Lucia, D. J., Beran, P. S., and Silva, W. A., “Reduced-order modeling: new approaches for computational physics,” Progress in Aerospace Sciences, Vol. 40, No. 1-2, 2004, pp. 51–117.
  • Durox et al. [2009] Durox, D., Schuller, T., Noiray, N., and Candel, S., “Experimental analysis of nonlinear flame transfer functions for different flame geometries,” Proceedings of the Combustion Institute, Vol. 32, No. 1, 2009, pp. 1391–1398.
  • Noiray et al. [2008] Noiray, N., Durox, D., Schuller, T., and Candel, S., “A unified framework for nonlinear combustion instability analysis based on the flame describing function,” Journal of Fluid Mechanics, Vol. 615, 2008, pp. 139–167.
  • Chatterjee [2000] Chatterjee, A., “An introduction to the proper orthogonal decomposition,” Current science, 2000, pp. 808–817.
  • Lucia and Beran [2003] Lucia, D. J., and Beran, P. S., “Projection methods for reduced order models of compressible flows,” Journal of Computational Physics, Vol. 188, No. 1, 2003, pp. 252–280.
  • Xu and Duraisamy [2017] Xu, J., and Duraisamy, K., “Reduced-Order Modeling of Model Rocket Combustors,” 53rd AIAA/SAE/ASEE Joint Propulsion Conference, 2017, p. 4918.
  • Xu et al. [2018] Xu, J., Huang, C., and Duraisamy, K., “Multi-Domain Reduced-Order Modeling with Sparse Acceleration of Combustion Instability,” 2018 Joint Propulsion Conference, 2018, p. 4680.
  • Huang et al. [2018a] Huang, C., Xu, J., Duraisamy, K., and Merkle, C., “Exploration of Reduced-Order Models for Rocket Combustion Applications,” 2018 AIAA Aerospace Sciences Meeting, 2018a, p. 1183.
  • Huang et al. [2016a] Huang, C., Anderson, W. E., Merkle, C. L., and Sankaran, V., “Multi-Fidelity Framework for Modeling Combustion Instability,” Tech. rep., AFRL/RQR Edwards AFB United States, 2016a.
  • Huang et al. [2017] Huang, C., Anderson, W. E., and Merkle, C., “Multi-fidelity Framework Explorations for Nonlinear Euler Equations,” 53rd AIAA/SAE/ASEE Joint Propulsion Conference, 2017.
  • Rowley et al. [2004] Rowley, C. W., Colonius, T., and Murray, R. M., “Model reduction for compressible flows using POD and Galerkin projection,” Physica D: Nonlinear Phenomena, Vol. 189, No. 1-2, 2004, pp. 115–129.
  • Ravindran [2000] Ravindran, S. S., “A reduced-order approach for optimal control of fluids using proper orthogonal decomposition,” International journal for numerical methods in fluids, Vol. 34, No. 5, 2000, pp. 425–448.
  • Amsallem and Farhat [2008] Amsallem, D., and Farhat, C., “Interpolation method for adapting reduced-order models and application to aeroelasticity,” AIAA Journal, Vol. 46, No. 7, 2008, pp. 1803–1813.
  • Lieu and Farhat [2007] Lieu, T., and Farhat, C., “Adaptation of aeroelastic reduced-order models and application to an F-16 configuration,” AIAA Journal, Vol. 45, No. 6, 2007, pp. 1244–1257.
  • Barbagallo et al. [2012] Barbagallo, A., Dergham, G., Sipp, D., Schmid, P. J., and Robinet, J.-C., “Closed-loop control of unsteadiness over a rounded backward-facing step,” Journal of Fluid Mechanics, Vol. 703, 2012, pp. 326–362.
  • Barbagallo et al. [2011] Barbagallo, A., Sipp, D., and Schmid, P. J., “Input–output measures for model reduction and closed-loop control: application to global modes,” Journal of Fluid Mechanics, Vol. 685, 2011, pp. 23–53.
  • Huang et al. [2019a] Huang, C., Anderson, W. E., Merkle, C., and Sankaran, V., “Multifidelity Framework for Modeling Combustion Dynamics,” AIAA Journal, 2019a.
  • Smith et al. [2008] Smith, R., Ellis, M., Xia, G., Sankaran, V., Anderson, W., and Merkle, C., “Computational investigation of acoustics and instabilities in a longitudinal-mode rocket combustor,” AIAA Journal, Vol. 46, No. 11, 2008, pp. 2659–2673.
  • Yu et al. [2012] Yu, Y., Sisco, J. C., Rosen, S., Madhav, A., and Anderson, W. E., “Spontaneous longitudinal combustion instability in a continuously-variable resonance combustor,” Journal of Propulsion and Power, Vol. 28, No. 5, 2012, pp. 876–887.
  • Yu [2009] Yu, Y. C., “Experimental and analytical investigations of longitudinal combustion instability in a continuously variable resonance combustor (CVRC),” Ph.D. thesis, Purdue University, 2009.
  • Huang et al. [2019b] Huang, C., Duraisamy, K., and Merkle, C., “Investigations and Improvement of Robustness of Reduced-Order Models of Reacting Flow,” AIAA Scitech 2019 Forum, 2019b, p. 2012.
  • Huang et al. [2018b] Huang, C., Duraisamy, K., and Merkle, C., “Challenges in Reduced Order Modeling of Reacting Flows,” 2018 Joint Propulsion Conference, 2018b, p. 4675.
  • Frezzotti et al. [2015a] Frezzotti, M. L., Nasuti, F., Huang, C., Merkle, C., and Anderson, W. E., “Response Function Modeling in the Study of Longitudinal Combustion Instability by a Quasi-1D Eulerian Solver,” 51st AIAA/SAE/ASEE Joint Propulsion Conference, 2015a, p. 3840.
  • Frezzotti et al. [2015b] Frezzotti, M. L., Nasuti, F., Huang, C., Merkle, C., and Anderson, W. E., “Parametric Analysis of Response Function in Modeling Combustion Instability by a Quasi-1D Solver,” 6th EUROPEAN CONFERENCE FOR AEROSPACE SCIENCES, 2015b.
  • Crocco et al. [1958] Crocco, L., Grey, J., and Harrje, D., “On the importance of the sensitive time lag in longitudinal high-frequency rocket combustion instability,” Jet Propulsion, Vol. 28, No. 12, 1958, pp. 841–843.
  • Parish et al. [2018] Parish, E., Wentland, C., and Duraisamy, K., “A Residual-Based Petrov-Galerkin Reduced Order Model with Memory Effects,” Submitted, Journal of Computational Physics, 2018.
  • Roe [1986] Roe, P. L., “Characteristic-based schemes for the Euler equations,” Annual review of fluid mechanics, Vol. 18, No. 1, 1986, pp. 337–365.
  • Grenda et al. [1995] Grenda, J., Venkateswaran, S., and Merkle, C., “Application of computational fluid dynamics techniques to engine instability studies,” PROGRESS IN ASTRONAUTICS AND AERONAUTICS, Vol. 169, 1995, pp. 503–528.
  • Wang et al. [2018] Wang, Q., Hesthaven, J. S., and Ray, D., “Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem,” Tech. rep., 2018.
  • Huang et al. [2016b] Huang, C., Anderson, W. E., Harvazinski, M. E., and Sankaran, V., “Analysis of self-excited combustion instabilities using decomposition techniques,” AIAA Journal, 2016b, pp. 2791–2807.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description