1. Introduction
With the development of modern propulsion technology, both the energy efficiency and the energy density reach a considerably high level. Combustion instability has been an unsolved problem since the early application of both liquid rocket engines and solid rocket motors (SRMs) [
1,
2].
The sources of the combustion instabilities in SRMs are believed to be the interaction of the flow field, the acoustic fields and the unsteady combustion inside the chamber [
1]. The coupling between the unsteady combustion and the acoustic fields, or the thermoacoustic mechanism, plays a central role in various SRMs, especially in small and medium SRMs [
3]. While for the large scale SRMs with high length-to-diameter ratio, the hydrodynamical instability governs.
The thermoacoustic instability in SRMs has been attracting most of the researchers in the community since the 1960s. According to Culick [
4], a reliable modeling of the thermoacoustic instability should properly include the non-linear effects in the acoustic system of the chamber and the combustion system of the propellants. Actually, the three combinations of linear acoustics with linear combustion response, the linear acoustics with non-linear combustion response and the non-linear acoustics with linear combustion response, are not able to model the important phenomenon of bi-stability and triggering. Following Culick’s approach [
5,
6], the non-linear acoustics in the chamber can be well-modeled. Nevertheless, on the other side, the modeling of the unsteady combustion remains difficult until now. Traditional approaches are based on different overall response functions [
4]. The major limit of this approach is the historical effects of temperature distribution inside the solid propellants. On account of this, Mariappan and Sujith [
7] improved the model by replacing the overall response function by a spatial position specified unsteady combustion model, yielding more detailed non-normal and non-linear features. In their model, the spatial temperature distribution in the solid propellant is taken into consideration, which leads to the historical effects that are included. Indeed, under the same pressure disturbance, the solid phase will respond differently with different temperature distribution, which is determined by the combustion history of the propellant. The unsteady combustion model used by Mariappan and Sujith [
7] is the quasi-steady homogeneous one-dimensional (QSHOD) model equipped with a flame model, named the Krier-T’ien-Sirignano-Summerfield (KTSS) model [
8], where the flame is in the gas phase, which is assumed to be a uniform one in a region above the burning surface. The author of the present paper finds that the Zel’dovich-Novozhilov(ZN) model serves better, and investigates the triggering instability and bi-stable region thoroughly in a previous work [
9,
10].
Unlike the steady combustion models, the unsteady combustion model is designed to be used in the combustion instability study of SRMs. The detailed modeling of the various physical and chemical processes in the steady model is no longer appropriate due to its mathematical complexity. Thus, the unsteady model should not be too simplified so as to include the major dynamics related to the combustion instability, and also not too complicated to be mathematically intractable [
1]. The unsteady combustion models of solid propellants were developed separately by the researchers in the USA and the former USSR/Russia. Researchers in the USA adopt the QSHOD approach, where the processes in the condensed phase, the gas phase and the interface between them are detailed, modeled and then combined together [
1]. Meanwhile, the Russian researchers follow the phenomenological approach, originally proposed by Zel’dovich, where the burning rates and surface temperature are related to the ambient pressure and temperature directly [
11]. The comparison of the two classes of models have been discussed thoroughly since the 1990s, and this is summarized in [
12], where the linear and non-linear responses of the two models under pressure oscillations with relative high amplitude are compared. It is well known that the QSHOD model and ZN model are equivalent in the linear region mathematically. In [
12], the non-linear responses are compared based on this linear equivalence, which leads to the inevitable drawback that the equivalence is derived under linear assumption while the non-linear responses are investigated.
On account of this, the present paper returns to the original phenomenological construction of the models. The steady combustion features of solid propellants based on the QSHOD model are thoroughly calculated in a wide range of ambient pressure and temperature. The results are then used to determine the phenomenological parameters of the ZN model, yielding a set of ZN parameters for the wide range of ambient conditions. The responses of the ZN model with the parameters calculated from traditionally linear equivalence and the new phenomenological approach are compared, from the linear to non-linear regimes.
After the present introductory part, the mathematical formulations of the QSHOD model and the ZN model, as well as the relations between them, are discussed in
Section 2. Based on the mathematical characteristics of the two models,
Section 3 discusses the numerical algorithm used to investigate the equations. Using this algorithm, a phenomenological modeling of a QSHOD model is constructed, and the responses of both models are illustrated and discussed in
Section 4. Finally, the concluding remarks are given in
Section 5.
2. Mathematical Formulations
The mathematical formulations of the QSHOD model and ZN model share most parts in common, which has been discussed thoroughly in the literature [
12]. In the present section, the governing equations and the boundary conditions are discussed, as well as their physical significance. However, not all of the derivations and assumptions are included here.
2.1. QSHOD Model with KTSS Flame
To start with, the non-dimensionalized governing equation of the model is [
8]
where
are the non-dimensional temporal and spatial variables,
are the non-dimensional temperature and pressure, and
R is the non-dimensional burning rate. It should be noted that the exponential dependence of the burning rate on the surface temperature is an approximation of the exponential Arrhenius’ Law. Theses variables are defined as follows.
where
H represents a non-dimensional energy comparing the total heat release of the combustion
and the heat capacity of the propellant. The condition under investigation is determined by the ambient pressure
and temperature
. Under this condition, the steady burning rate and the burning surface temperature are
,
, respectively. The spatial and temporal variables are non-dimensionalized by
and the heat diffusivity of the propellant
.
The QSHOD model is in a class of similar models with the same governing equations [
12]. The QSHOD-KTSS model is a typical one of the class, where the above governing equation for the condensed phase is equipped with a flame in the gas phase, assumed to be uniform in a region. Under this assumption, the corresponding boundary condition can be derived as [
8]
where the dependence on pressure is introduced. The boundary condition is derived from the mass and energy balance at the phase interface [
8]. It should be noticed that the governing Equation (
1) is a partial differential equation (PDE) with strong coupling with its boundary condition (
2). The equation is derived from a heat conduction problem with the regressing boundary, leading to a term with a coefficient
R determined by the boundary value of temperature. Thus, the equation governs the temperature field on a moving domain, and the velocity of the moving boundary, which is also the burning rate
R, which is determined by the temperature field. This complicated interaction situation precludes the possibility of an analytical solution. The previous study of the model focus on the analytical solution of its linearization system or on the numerical solution of the full system.
From the equations above, the three major parameters of the QSHOD-KTSS model are n, m, and H. The former two determine the dependence of the burning rate on the ambient pressure and burning surface temperature. H describes the energy released by the combustion of the propellant to the heat capacity.
2.2. ZN Model
The ZN model is a phenomenological approach without detailed modeling of the flame following the notion of Zel’dovich [
11,
13,
14]. The governing equation for the solid phase is the same as the QSHOD model,
where the non-dimensional variables are defined as the same as the QSHOD model.
According to Zel’dovich’s notion, there is no need to model the processes in the gas phase in detail. Alternatively, a phenomenological approach in the macroscopic perspective is adopted. To start with, there are two major features of steady combustion, the steady burning rate
and the steady burning surface temperature
, and two major ambient conditions, ambient pressure
and temperature
. In steady combustion conditions, there exists dependence relations between them.
The specific forms of the functions
are not known, and they are not necessary to be known. Usually, a derived form of the above relations is used, based on the steady distribution of temperature fields in the solid propellant [
14]. It should be noticed that the major feature of the above relation can be described by the following partial derivatives.
The central assumption of the ZN model is Zel’dovich’s notion that the chemical and gas dynamical processes in the gas phase are much faster than the heat conduction processes inside the condensed phase. Under this condition, the relaxation times of the gas phase are much smaller than those in the condensed phase, and the gas phase reaches equilibrium immediately, as the temperature in the condensed phase varies. Mathematically, the Equation (
4) can be extended to the unsteady case [
11]. For example, the burning rate
under a pressure disturbance
can be expressed as
. In this manner, the burning rate and the burning surface temperature of the unsteady combustion can be expressed by the parameters
, defined in Equation (
5). The burning rate and boundary condition of the ZN model can be derived as
where the heat flux at the burning surface is
. The readers can refer to the literature [
11,
14] for detailed derivation processes.
Thus, the ZN model, consisting of Equations (
3) and (
6) and the QSHOD model, consisting of Equations (
1) and (
2) have the same governing equation as the two models share the same solid phase modeling. The burning rate expressions and boundary conditions in the two models are completely different, as the modeling methods for the gas phase and phase interface are different. Despite the obvious difference, it is well known that the linearization of the two models yields the equivalent results [
1,
14].
In the previous work [
12], the two models are also compared based on a linear equivalence. The model variables of the QSHOD model are
, while the model variables of the ZN model are
. Except for the pressure exponent
n, the equivalent relation of the other two variables is
which can be derived in the sense that the derivative of the burning rate on the surface temperature is the same [
12]. Obviously, this linear equivalence holds near the reference state, where the derivatives are taken, which may not valid for a wide parameter range. Thus, a new approach is proposed in the present paper.
2.3. The Phenomenological Reconstruction of the QSHOD Model
From the discussion in the previous subsection, the linear equivalence holds accurately only on one state point. In the present paper, the calculations are taken in a wide state range, and the ZN parameters are evaluated based on the information throughout the range comprehensively.
The method adopted in the present paper is based on the phenomenological notion of the ZN model. The steady state combustion states are evaluated by the QSHOD model with varying ambient pressures and temperatures, and the ZN parameters are then calculated according to the definition (
5). In the QSHOD-KTSS model, the two major coefficients determining the main feature of the model are
H and
m. Thus, with the fixed
, the steady combustion can be solved by evolving the following equation to the steady state,
where only the temperature is reverted to the dimensional form to take the ambient temperature variation into consideration.
For each ambient state,
, one can solve a steady state from the above equations, characterized by the burning rate and surface temperature
and
. Based on the two matrices of
and
, the parameters can be evaluated by the rows and columns. For example,
k can be evaluated by each row of the matrix
, where the information in row
will yield a
k at fixed ambient pressure
with varying ambient temperature
by a regression algorithm. In the present paper, an
mesh of ambient states are adopted, which can be adjusted for specific ranges as needed. Solving the above system for
different
P and
different
yields an
matrix of
and an
matrix of
, the four coefficients in Equation (
5) can then be calculated by regression in each row and column (the total number for each parameter is
). In this manner, a ZN model based on the results of the QSHOD-KTSS model is derived. The detailed modeling of the flame model is replaced by the phenomenological relation Equation (
4) or Equation (
5).
4. Results and Discussion
In this section, the unsteady combustion under varying pressure oscillations of the models will be discussed. The models in the comparison are the QSHOD model and the two equivalent models, in the sense of linear equivalence and the phenomenological approach, respectively.
4.1. Model Parameters
The basic parameters for the solid propellant under investigation are collected in
Table 1. The calculations are conducted in the non-dimensional form, according to (
8), and the results are reverted back and summarized in the dimensional form.
4.2. The Equivalent ZN Models
Traditionally, the equivalent counterpart of a QSHOD model in the ZN framework should be based on the linear equivalence [
12,
15]. In the present model, the corresponding ZN parameters can be calculated by Equation (
7)
where the superscript denotes the linear equivalence.
In the present paper, a new phenomenological approach is shown in
Figure 1, where the logarithmic coordinates are used when the logarithm of the quantity is used in Equation (
5). In each subfigure, five typical conditions are selected to draw the variation. It can be seen from the figure that the lines are of good linearity, and the slope of each line is almost the same. This indicates that, in the selected parameter range, the ZN parameters are close to constant. However, this does not imply that the above linear equivalence method is valid.
In the present paper, the ZN parameters are calculated according to the algorithm in
Section 3.2. Under the
mesh, there are
results for each parameter, and their mean value is used for the final results,
It is noticeable that the parameters generated from the new approach are quite different from the traditional linear equivalence results
, and the pressure exponent is a bit larger.
Comparing the two methods, the linear equivalence is derived from the derivatives calculated at the reference state, while in the phenomenological method, the parameters are calculated from the comprehensive effects throughout the parameter range. The difference originates from two aspects. The first one is the effects of different parameters, which are coupled together, and not separated. The pressure exponent
n is an example. Although the pressure exponent is introduced into the QSHOD-KTSS model explicitly in the boundary condition (
8), its influence is not straightforward, and the pressure exponent of the full system differs from the original set value. Another important aspect is the wide state range. The four parameters are not constant throughout the state range and the dependences of
on
are not strictly linear. The results illustrate the discrepancies between the linear equivalence at a single state point and the overall variations in the state range.
4.3. The Unsteady Combustion of the QSHOD-KTSS Model
The comparison of the unsteady combustion of the QSHOD model and the ZN model has been discussed in [
12] under various conditions. In the present paper, the harmonic excitation with increasing amplitude is discussed.
In the harmonic excitation problem, the unsteady combustion is driven by a periodically varying ambient pressure, which is the case in the chamber of an SRM under combustion instability. The driving signal is
where the non-dimensional form is used. According to Equation (
11), the pressure is oscillating around the reference pressure at an amplitude
.
It has been shown in the literature that the non-linear dynamics of the unsteady combustion under the increasing driving amplitude is a series of period-doubling bifurcations to chaos [
11,
14]. In the present paper, the comparison is also made in this problem. In each calculation, time is evolved to
and the range
is plotted to cut off the transient processes. The reason for cutting off the transient processes is the non-linear dynamics of the system, which varies mainly on the final states, which is the attractor of the system. Indeed, it will be shown that the combustion system will be attracted to a series of geometrically different attractors with increasing amplitude.
To start with, at
, which corresponds to a steady constant pressure condition, the burning rate response is also steady at unity. With an increasing
, there will be a burning rate oscillation, correspondingly. When the amplitude is small, the linear simplification applies, and a harmonic pressure perturbation will lead to a harmonic burning rate response with the same frequency. The burning rate responds to the pressure oscillation in a “linear” manner when
, as shown in
Figure 2a. In the three subfigures, the top one is the time series of the burning rate, the middle one is the phase portrait by phase reconstruction and the bottom one is the FFT of the time series. Actually, the response in this region has begun to deviate from the linear response, but the period of the burning rate is the same as the driving pressure. Mathematically, the attractor of the system is qualitatively the same for
.
At the point
, a period-doubling bifurcation occurs, where the limit cycle becomes self-intersected and the period of the cycle doubles, as shown in
Figure 2b. As the attractor becomes geometrically different, the final state of the system is qualitatively different from the
case, and a bifurcation occurs. With the further increasing of the
, there will be a series of period-doubling to period-4 (
), period-8, period-
and finally to chaos (
), as shown in
Figure 2c,d.
It is noticeable that the burning rate response becomes extremely high as , which is unphysical. This is due to the over-simplification of the flame structure, which may lead to wrong results when applied in the combustion instability study of the SRMs.
4.4. The Unsteady Combustion of the Two ZN Models
Based on the linear equivalency, the ZN parameter should be
,
,
. The corresponding unsteady combustion responses are collected in
Figure 3. The overall bifurcation scenario is the same as the QSHOD model, qualitatively. Nevertheless, from the quantitative perspective, there are considerable differences. The bifurcation points are smaller, or each period-doubling becomes earlier. Additionally, the burning rate responses are confined in a realistic range. Thus, the ZN model should be more appropriate for the combustion instability study than the QSHOD model, as discussed in [
12]. However, although the ZN model by linear equivalence is equivalent to the QSHOD model according to non-linear dynamics, in the sense of the bifurcation scenario, the quantitative features are quite different; especially the exact position of the bifurcation points.
In the new phenomenological approach proposed in the present paper, the ZN parameter should be
. The corresponding unsteady combustion responses are shown in
Figure 4. Again, the overall bifurcation scenario is the same as the previous two cases. But, the bifurcation points of the new model become more close to the original QSHOD model than the ZN model based on the linear equivalency.
Indeed, the typical bifurcation points of the three systems are collected in
Table 2, and the ZN model from linear equivalence has much smaller bifurcation points. It should be noted that the unphysical large values of both bifurcation points and response amplitudes in the QSHOD-KTSS model are due to the simplifications in the modeling processes. For example, the burning rate
is a simplification of the exponential function in Arrhenius’ Law, which is valid near the reference state or
. Admittedly, the accuracy of the QSHOD model is not the focus of the present paper. However, the results show that the unphysical late bifurcation effects are better captured by the phenomenological approach than the linear equivalence approach.
By the comparison of the two ZN models, there also exists the non-linear equivalency in the sense of the same bifurcation scenario. The ZN model from the phenomenological approach reconstructs the system better, as each bifurcation point is closer to the original QSHOD model than the linear equivalency. Indeed, the phenomenological approach makes full use of the information in a wide range of parameters, while the linear equivalency only use the information near the reference state.
5. Conclusions
A new unsteady combustion modeling framework is investigated based on the Zel’dovich phenomenological approach. The original unsteady combustion model based on flame modeling is reconstructed in the phenomenological form, with both good qualitative and quantitative agreements.
The spectral method is modified to investigate the combustion features of the QSHOD-KTSS model in a wide state range, yielding combustion characteristics under different ambient conditions. The results are then analyzed by linear regressions to calculate the four ZN parameters, in a phenomenological manner.
The unsteady burning rate responses of the the new constructed model are then compared with the traditional ZN model based on linear equivalence. The results show that the new phenomenological approach is able to reconstruct the original model better than the traditional linear equivalence model, in the sense that the bifurcation scenario is closer. This indicated that the phenomenological approach, in which information in a widely varying state range is taken into consideration, is able to capture dynamical properties of the original model better than the linear equivalence approach, in which information in only a single state is considered.
The proposed methodology can be used for more complicated combustion modeling for solid propellants, or even liquid propellants. Indeed, the modeling of unsteady combustion responses of different solid or liquid propellants focuses on the burning rate responses of various perturbations of ambient conditions. The phenomenological approach can be applied in these modeling problems.