Reduced Order Modeling Framework for Combustor Instabilities Using Truncated Domain Training
Abstract
A multifidelity 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 projectionbased 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. Geometryspecific 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 quasi1D model of a singleelement rocket combustor and the present framework is compared to traditional ROM development approaches.
Reduced Order Modeling Framework for Combustor Instabilities Using Truncated Domain Training
Nomenclature
=  amplification factor 
=  cross sectional area 
=  stoichiometric fueltooxidizer 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 fueltooxidizer 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 laboratoryscale 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 fullscale engines [9, 10, 11]. This motivates research on alternative approaches such as reducedfidelity modeling [12, 13, 14, 15, 16] and reduced order modeling [17, 18].
Reducedfidelity modeling approaches usually adapt certain levels of simplifications in physical and numerical models to achieve higher computational efficiency. Sirignano and Popov [14] developed a twodimensional model for transversemode combustion instability in a cylindrical rocket motor and further extended it to rectangular configurations [15]. Though reducedfidelity 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 wellcharacterized flow variables such as pressure and velocity. This class of techniques include the socalled 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 swirlstabilized 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 lowamplitude 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. Projectionbased ROMs [21, 22] have recently been applied to combustion instability problems [23, 24, 25, 26, 27]. In projectionbased 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 projectionbased ROMs on representative combustion model problems can be found in Refs. [23, 25, 26, 27]. It should be pointed out that PODbased ROMs have been more widely used in nonreacting 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 regenerated or reparameterized, which restricts the scope of ROM applications.

PODROMs 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 highenergy regions. The dynamics in lowmagnitude regions can be underresolved and consequently affect stability and robustness.
The present work extends the multifidelity framework methodology developed by Huang et al. [26, 27, 34] and assesses ROM performance on the major aspects listed above based on a quasi1D model combustor problem, which was originally developed by Smith et al. [35] to model the experimental work by Yu et al. on the singleelement 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 flamevortex 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 multidomain approach where the computational domain is divided into two components: the upstream one covering the critical physics in the heat addition/acousticflame 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 multifidelity solver ^{1}^{1}1For the current demonstration, we replace the reducedfidelity part with the full order model for consistency of evaluations. The multifidelity 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 quasi1D 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 quasi1D 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 quasi1D assumption, the backstep and the converging part of the nozzle are sinusoidally contoured in this study, in contrast to a discontinuous change in the crosssectional 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  Backstep  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 
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 
The governing equations in this problem are the quasi1D unsteady Euler equations with species transport, represented as
(1) 
where
(2) 
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 backstep 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 infinitelyfast onestep 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
(3) 
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 infinitelyfast model, the consumption rate of the oxidizer is
(4) 
where is the stoichiometric fueltooxidizer 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.
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 .
(5) 
In the simulation, a steady state is first achieved with the unsteady source term turned off. From the steady state, a lowamplitude 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 selfexcited 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 highfidelity solution in the current work, and the corresponding simulation is referred as FOM.
ii.b Projectionbased ROM
To derive the projectionbased ROM equations, we first rewrite the full order model given by Eq. (1) for each variable as
(6) 
where the variable index corresponds to the four conserved variables respectively. The residual vector for variable consists of the corresponding flux and source terms
(7) 
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:
(8) 
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 columnwise from the first leftsingular vectors from the SVD ^{2}^{2}2The remaining leftsingular vectors are . The conserved variables are then approximated as
(9) 
where is the th element in the reduced order variable vector , and is the th column in the POD basis . Since each leftsingular 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 Galerkinprojected ROM for coefficients is then
(10) 
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 multidomain modeling.
Iii Framework Overview
In this section, we explore characteristicbased ROM training on a truncated domain and contrast it with conventional ROM development. A schematic is given in Fig. 2. The major steps include

Truncate the domain into two subdomains, the first containing the physicscomplex area that covers the fuel injection and the flame, and the second containing the variablelength chamber dominated by acoustics. The ROM will be trained for the first subdomain.

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

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

Use the POD bases in a ROM solver for the first subdomain and couple it with a FOM solver for the variablelength, acousticsdominated chamber for predictions.
It can be noted that in the framework ROM is only developed for the first subdomain, whereas the second subdomain is solved in FOM. There are two reasons for this choice:

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

In design evaluations, the chamber length in the second subdomain is variable, thus a new ROM training is required for each chamber length, which violates the goal of reducing the computing cost using ROM.
iii.a Characteristic ROM training on reduceddomain
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 subdomains with the interface at m. The first subdomain includes the injector, backstep 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 subdomain, 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 reduceddomain 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 selfexcited 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:
(11) 
where is the speed of sound, is onedimensional approximation of the characteristic invariant for the incoming 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 Multidomain coupling
Following the aforementioned simulations, the bases for the conserved variables are obtained from the reduced domain solution snapshots, and used in PODGalerkin projection (Eq. (10)) to derive the ROM for the first subdomain. The second subdomain 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 subdomains communicate at every timestep 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 subdomain, and superscript and for the first and last cell in the corresponding subdomain respectively, the flux at the interface is given by
(12) 
where is the Roeaveraged 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 multidomain solver, i.e. the first subdomain 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 fulldomain FOM simulation is performed instead of the proposed characteristic training on the truncateddomain. Then the solution is restricted to the first subdomain and collected to generate the POD bases for the ROM of the subdomain.
 Control group B

uses a traditional fulldomain ROM. The same FOM training simulations as in control group A are used and the POD bases are directly generated on the fulldomain 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 ^{3}^{3}3although as shown in the next section, this can be reduced using characteristic perturbations. Then the first subdomain is simulated using the ROM, second subdomain using the FOM. In control groups A&B, selfexcited FOM simulations are conducted. In group A, the first subdomain 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 “reduceddomain training and multidomain solver (RDMD)", control group A as “fulldomain training and multidomain solver (FDMD)", and control group B as “fulldomain training and fulldomain solver (FDFD)". A schematic of the different methods is given in Fig. 3.
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
(13) 
which is a superimposition of frequencies, starting from with an incremental interval .
To train the control groups, fulldomain selfexcited FOM simulations are conducted for each combination of and . In all three methods, snapshots of the training solution are collected every 100 timesteps 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 fulldomain 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 peaktopeak amplitude of the unsteady part of the pressure signal, whose growth is fitted to a exponential function plotted in Fig. 4, given by
(14) 
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 spatiotemporal 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 inphase, 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 outofphase, which is recognized as a combustion instability damping mechanism.








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 subdomain, 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
(15) 
It should be noted that there is only one curve for RDMD as it uses the same training simulation and POD bases for all chamber lengths. It can be observed that for both FDMD and FDFD, 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 FDFD is slower than FDMD, 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 mediumtohigh , the proposed framework gives the best decay since its training data contains more highfrequency contents. At small , the decay in the leading modes for the two fulldomain training methods become slightly better than RDMD as the response frequency in their training data is comparable to the highest one included in the reduceddomain training, whereas there are also lower frequencies contained in the latter. Results for other variables and follow a similar pattern.
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 FDFD is on the full domain and only truncated in the plot to the first subdomain for comparison purposes. The first mode for all methods appears to be identical. For lower energy modes, the basis for the RDMD is quite different from the other two methods. The basis for the FDMD and FDFD 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 subdomains of a 3D simulation of the CVRC is compared (equivalent to the FDMD and FDFD 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
(16) 
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 largescale 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 highfrequency dynamics contained in the modes.
It can be noticed from Fig. 8 that the spatial modes that exhibit highfrequency 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 twofold: 1) using the same singlefrequency training, FDMD shows a faster decay in the singular values, which indicates that the lowfrequency highenergy dynamics is better captured in the leading modes, and thus the highfrequency modes are realized at smaller mode numbers than in FDFD. 2) There are multiple frequencies in the reduceddomain training. These frequencies are still significantly higher than those of the corresponding POD mode numbers in FDFD and FDMD. The corresponding modes representing these frequencies cause the RDMD 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 RDMD, the frequency increases with mode number steadily before the 35th mode. This demonstrate a good correlation between the temporal and spatial frequency captured in the modes. In contrast, many highspatialfrequency modes from the two fulldomain training methods actually evolve at a low temporal frequency. This inconsistency indicates a less efficient basis compared to RDMD.








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
(17) 
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 RDMD and FDMD provide a monotonic decrease in the projection error. However for FDFD, 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 multidomain 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.



iv.c ROM evaluations
Before integrating the ROM into the Multifidelity Framework, the predictive capability of the ROM obtained from the reduceddomain 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.
iv.d Multifidelity 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 RDMD at is comparable to that for FDFD 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 , RDMD performs similarly to FDMD, except for one case at m, where the latter is unstable. At high , where the instability is stronger, the advantage of RDMD becomes more significant. This can be observed from the medium cases at , where the gap between RDMD and FDMD is apparent. Also, FDFD performs worse as increases, resulting in higher errors and more unstable cases. It should be noted that FDFD 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 fulldomain training methods become unstable in most conditions. In contrast, RDMD 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 reduceddomain 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 timemarching scheme is used. In contrast to the traditional fulldomain 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.



Major quantities of interest in rocket combustor design include the dominant acoustic frequencies, the growth rates and LCO peaktopeak amplitudes of the pressure oscillations. It should be noted that the dominant acoustic frequencies have been wellpredicted 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 RDMD 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.












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 backstep area for better visualization, and are selected by the end of the simulation, when the error is maximum. It is observed that the RDMD solution correlates well with the FOM, and connects smoothly across the interface.



iv.e Offdesign condition performance
To further assess the RDMD framework, the ROM trained using is evaluated at several offdesign 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 peaktopeak 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 offdesign 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 RDMD 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 multifrequency 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 




V Conclusions
In this work, a multifidelity ROM development framework is investigated on a quasi1D 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 subdomains with the first subdomain containing the main reaction dynamics and the second covering the adjustable chamber of the CVRC. Characteristic training is conducted in the first subdomain, resulting in a ROM that can be directly integrated with a FOM solver of the second subdomain 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 mediumtohigh 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 peaktopeak amplitude shows that the framework is able to predict the these quantities accurately at all combinations of and at .
In summary, the multifidelity 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 multidomain approach also enables potential developments such as the application of different timemarching 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 multifidelity framework.
Vi Acknowledgments
The authors acknowledge support from the Air Force under the Center of Excellence grant FA95501710195, titled “MultiFidelity Modeling of Rocket Combustor Dynamics.”
References
 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 selfexcited unstable combustor,” Physics of Fluids (1994present), Vol. 27, No. 4, 2015, p. 045102.
 Matsuyama et al. [2016] Matsuyama, S., Hori, D., Shimizu, T., Tachibana, S., Yoshida, S., and Mizobuchi, Y., “LargeEddy Simulation of HighFrequency Combustion Instability in a SingleElement 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 flamebases 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., “Largeeddy simulation of a jetinhotcoflow burner operating in the oxygendiluted 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 singleelement 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 rocketlike configuration,” Combustion and Flame, Vol. 156, No. 6, 2009, pp. 1166–1180.
 HernándezPérez et al. [2011] HernándezPérez, F., Yuen, F., Groth, C., and Gülder, Ö., “LES of a laboratoryscale turbulent premixed Bunsen flame using FSD, PCMFPI 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 highpressure, shearcoaxial 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. 12, 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., “Loworder 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., “Twodimensional model for liquidrocket 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 gasturbine combustion instability analysis,” Combustion Science and Technology, Vol. 177, No. 56, 2005, pp. 1109–1150.
 Lieu et al. [2006] Lieu, T., Farhat, C., and Lesoinne, M., “Reducedorder fluid/structure modeling of a complete aircraft configuration,” Computer methods in applied mechanics and engineering, Vol. 195, No. 4143, 2006, pp. 5730–5742.
 Lucia et al. [2004] Lucia, D. J., Beran, P. S., and Silva, W. A., “Reducedorder modeling: new approaches for computational physics,” Progress in Aerospace Sciences, Vol. 40, No. 12, 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., “ReducedOrder 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., “MultiDomain ReducedOrder 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 ReducedOrder 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., “MultiFidelity 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., “Multifidelity 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. 12, 2004, pp. 115–129.
 Ravindran [2000] Ravindran, S. S., “A reducedorder 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 reducedorder 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 reducedorder models and application to an F16 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., “Closedloop control of unsteadiness over a rounded backwardfacing 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 closedloop 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 longitudinalmode 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 continuouslyvariable 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 ReducedOrder 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 Quasi1D 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 Quasi1D 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 highfrequency 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 ResidualBased PetrovGalerkin Reduced Order Model with Memory Effects,” Submitted, Journal of Computational Physics, 2018.
 Roe [1986] Roe, P. L., “Characteristicbased 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., “Nonintrusive 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 selfexcited combustion instabilities using decomposition techniques,” AIAA Journal, 2016b, pp. 2791–2807.