Abstract
To acquire the erosive characteristics of two-channel combustion chambers, a quasi-one-dimensional internal ballistics model is proposed. Combining the model with two different erosive burning models, the modified L-R expression, and the Mukunda-Paul expression, the flow parameters and the internal ballistics performance of a test solid rocket motor are computed. The results show good agreement with experimental data. According to the results, the more severe the erosion is, the earlier, longer, and gentler the tail-off stage becomes. During the tail-off stage, the dramatic drop of pressure leads to a low normal burning rate and makes it easier for erosion burning to occur. For this reason, notable erosive burning might appear during tail-off if the case-to-throat area ratio is extremely low. The results also show that flows in the inner channel and the outer channel are similar but not identical. This leads to different erosive burning behaviours in the two channels. Also, the two erosive burning rate models involved in this paper are compared. It seems that the M-P expression provides better results than the modified L-R expression does, since it reveals the threshold phenomena. Besides, the M-P expression has great advantages for its universality for most propellants and different SRM geometries.
1. Introduction
In a solid rocket motor (SRM) with a high load fraction, the high flow velocity often induces an augmentation of the burning rate, which is called “erosive burning.” Erosive burning has been the subject of many studies, yet a satisfactory interpretation has not been proposed because of the complexity of its nature. Nevertheless, some empirical models have been derived by studies involving experimentation and theoretical analysis.
Among these theories, the expression that gives a linear relationship between erosive burning rate and velocity might be the simplest one. Vandenkerckhove [1] proposed that the erosive burning rate is proportional to the difference between the crossover velocity and the “threshold velocity.” Based on the same erosive expression, Rout et al. [2] carried out one-dimensional (1D) numerical computation for rectangular and cylinder grains and discussed the influences of different factors on erosive burning. Many studies [3, 4] were analysed by Landsbaum [5], and the erosive burning rate was processed as a function of mass velocity (i.e., or ) beyond a threshold value.
Lenoir and Robillard [6] discussed the sources of heat to the burning surface and classified them into two categories, which are dependent on pressure and velocity, respectively. Correspondingly, the total burning rate of the propellant was written as the expression , where is the normal burning rate and is the erosive burning rate. Afterward, an implicit expression for the erosive burning rate was proposed, in which the distance to the head end of the grain is used as the characteristic length. Though the L-R expression is a rough theory which does not properly interpret the mechanism of propellant burning, it is widely accepted and used due to its good prediction.
Due to the L-R expression’s overestimation of erosive burning in large SRMs, some researchers have come up with a revised version of the L-R expression, where the hydraulic diameter is chosen to be the characteristic length instead of the length to the head of the grain [7–10]. Yet, to the authors’ knowledge, a rigorous derivation for the modification has not been given.
To validate the models referred to for a specific SRM, there are always constants to be confirmed. This makes it difficult to obtain a reliable internal ballistics prediction of an SRM before firing tests. Mukunda and Paul [11] proposed a dimensionless model in 1997. The ratio of total burning rate to normal burning rate was expressed as a function of the nondimensional flux. After vast comparisons with experimental data, parameters in the model were determined. This makes the model valid for most practical propellants. The model provides a reliable way to predict the internal ballistics performance for SRMs, with good accuracy even without a firing test.
In this paper, a two-channel combustion chamber with a noninhibitor grain is studied. A quasi-1D model is proposed, and, by combining this model with two different burning rate models (modified L-R expression and Mukunda-Paul expression), the authors calculate and discuss the internal ballistics performance of the combustion chamber and draw some interesting conclusions.
2. Mathematical Models
2.1. Quasi-1D Model for One-Channel Combustion Chambers
First, a steady, one-dimensional situation is considered for a combustion chamber with a single channel. As shown in Figure 1, a 1D control volume of the chamber has a length of , and the main flow goes from left to right. The transpiration flow has a velocity denoted as and a mass flow rate expressed as where is the density of the propellant, is the total burning rate of the grain, and is the area of the burning surface of this control volume. At the left section of the control volume, the mass flow rate, port area, pressure, velocity, and temperature of the flow are denoted as , , , , and , respectively. Due to the port area difference and mass flow from transpiration, these variables at the right section change into , , , , and .
The continuity equation can be written as where is the density of the flow.
The momentum equation for the control volume is where is the axial component of the transpiration velocity. By introducing the variable , Equation (3) can be written as where is the Mach number of the main flow at the left section of the control volume, is the specific heat ratio, and is gas constant of the combustion gas.
Assuming that the transpiration gas and the gas in the main flow have the same value of specific stagnation enthalpy, the energy equation can be written as where is the specific heat at constant pressure. Taking the derivative on both sides of Equation (5) and considering the expression , Equation (5) is modified as
From Equations (2), (4), and (6), the derivatives of , , and can be deduced as functions of and as follows:
Derivatives of other parameters such as and can be derived from Equations (7), (8), and (9).
The variable in Equations (7), (8), and (9) depends on the geometry of the grain. In most SRMs, the angle between the generating line of the grain and the axis is small, even under the condition of the erosive burning. This makes the value of to be relatively small compared to one. The values of in References [6, 7] are about 0.4% and 0.8%, for instance. Thus, is assumed to be zero, and the expression is considered as 1 in this paper.
Equations (7), (8), and (9) (with the consumption) offer a relationship of flow parameters at two adjacent sections. According to these equations, for a combustion chamber whose geometric parameters (i.e., areas of different sections) and burning rate model (which will determine the transpiration mass flow rate according to Equation (1)) are given, if there is some section where the flow parameters are known, the 1D flow field of the whole combustion chamber can be computed by an iterative procedure from one section to another. In this paper, we choose the aft-end section of the combustion chamber as the “known section.”
During motor operation, the velocity at the throat of the nozzle would be the sound velocity, and at a specific section in the nozzle, the flux function equals the area ratio and can be expressed as [12] where and represent the area of the specific section and throat, respectively. With the expansion loss neglected, the aft end of the chamber can be seen as the inlet of the nozzle. In such situation, Equation (10) can be used to compute the Mach number at the aft end of the chamber. Equation (10) yields two solutions for . Since we focus on the flow in the combustion chamber, only the subsonic value is taken. With the Mach number computed, the temperature at the aft end of the combustion chamber can be obtained from the relation where the combustion temperature is used as the total temperature of the combustion gas.
The density or pressure at the aft-end section is still need to be ascertained to obtain all the flow parameters at the aft-end section. To extend the 1D approach to SRMs with two channels, an estimation of mass flux at aft-end section is given at each time step to determine the density () and then pressure () at the aft-end section.
After an iterative procedure among sections, the parameters including mass flux at all sections in the chamber can be computed. At the head end of the chamber, the mass flux should equal to zero. This provides a criterion to determine whether the estimation mass flux is accurate. According to the procedure described above, the mass flux at the head end can be considered to be a function of the estimated value of the of aft-end mass flux, . The criterion can be expressed as an equation
The secant method can be used to solve Equation (12) to get a good estimation of and then the flow state of the whole combustion chamber.
Now that the flow parameters have been computed for a time step, we can compute the total burning rate and the geometry of the grain for next time step using a proper burning rate model. Then, the iteration goes to the next time step, until the chamber pressure drops below the ambient pressure.
2.2. Quasi-1D Model for Two-Channel Combustion Chambers
For a two-channel combustion chamber with a free-standing noninhibitor grain (which means all the surfaces of the grain burn during the SRM operation), a new situation emerges. In a combustion chamber of this kind, the flow parameters and the erosive conditions are different between the inner and outer channels. The exchange of gas between the two channels can often appear. Thus, for two-channel combustion chambers, the traditional 1D assumption cannot hold anymore.
We assume that the flow variables at the aft ends of the two channels are identical. The validity of this assumption can be proved by a steady two-dimensional computational fluid dynamics (CFD) procedure. Figure 2 shows the contour of axial velocity near the aft end of the combustion chamber of an SRM. The parameters used are from the SRM which will be described in Section 3. The Mukunda-Paul expression is used to obtain the burning rate of the grain. The values of the average velocities of the inner channel and outer channel at the aft-end section are 160 m/s and 155 m/s, respectively. The average pressures in the both channels are all 10.6 MPa.
Under the assumption, the Mach number in the two channels is computed from Equation (10), only the area is now the sum of the areas of the two channels (represented by and , respectively). As in one-channel combustion chamber conditions, at the aft-end section, an estimation of the total mass flow rate of both the channels is given. The mass flow rate at each aft-end section can be easily computed:
Two iterative procedures among the sections from aft ends to head ends can be separately carried on in the two channels. For each channel, the iteration will be just the same as that under the one-channel condition. The sum of the mass flow rate at the head-end section is now considered to be a function of . For a chamber with a noninhibitor grain, the sum of the mass flow rate at the head-end sections of the channels equals to the gas generated from the head-end burning surface . This is now the criterion to ascertain the true value of and is expressed as
The secant method is used to solve Equation (15), and then the true value of the total mass flow rate at aft-end section of the chamber is obtained. This yields the flow field variables of the whole combustion chamber at the current time step.
The burning rate and geometry of the grain are updated according to the flow field at the current time step, and the iteration goes on to the next time step until the pressure drops below the ambient pressure.
In this paper, two different erosive burning rate models are used. They are the modified L-R expression and Mukunda expression. For comparison, an iterative procedure using the de Saint Robert law, i.e., , is also carried out in this paper. The internal ballistics results obtained by the different computations are compared.
2.3. Modified L-R Expression
In Reference [6], the heat transfer coefficient for flow over a flat plate, written as is used to represent the heat transfer coefficient on the burning surface without transpiration, while the heat transfer coefficient under the transpiration condition is written as . The proportional relationship is postulated between the erosive burning rate and the heat transfer coefficient. Considering the expression of the normal burning rate, the total burning rate can be written as where is the proportional coefficient and is expressed as where is the specific heat of the propellant, is the temperature of the burning surface, and is the ambient temperature. Then, the L-R expression is deduced: where the parameter is
As indicated in Reference [8], in a regular solid rocket motor, the gas velocity will decrease as the burning surface regresses and the port area becomes larger, which leads to a decrease in erosive burning. Yet, the L-R expression cannot reveal the decrease of the erosive burning under this condition.
In this paper, a slight modification is realized by using the Dittus-Boelter equation as the heat transfer coefficient : where is the Reynolds number of tube flow and has the form of
As Equation (21) is suitable for tube flow, it should be more appropriate to use in SRM rather than Equation (16).
A procedure similar to the one in Reference [6] is carried out, and the total burning rate now can be expressed as
Equation (23) has a very similar format with Equation (19) proposed by Lenoir and Robillard except for the characteristic length and the expression of parameters. It should be noted that, while the parameter has the same value as in Equation (19), the first parameter in the modified expression is written as
It has a similar but not identical format with given in Equation (20).
A similar modified expression has been mentioned in References [7–10], but it seems that the researchers just emphasized the replacement of the characteristic length, and the difference between and was not discussed, so the expression for parameter has not been presented explicitly.
Like in Reference [6], the parameters and can be obtained simultaneously by a parameter identification procedure, or can be computed from Equation (24) if the relevant variables are known, followed by a trial-and-error procedure to determine .
2.4. Mukunda-Paul Expression
As we can see, the constants in the different forms of the L-R expressions are unknown for a regular motor. The firing test and trial-and-error procedure are needed to determine these parameters. Mukunda and Paul [11] proposed a dimensionless relationship to obtain the ratio of the total burning rate to the normal burning rate (). Strictly, parameters in the expression also need to be confirmed. Yet they found that for various propellants, the parameters can be close. With the recommended parameters, the expression can be considered universal with good accuracy [13, 14].
In their study, Mukunda and Paul raise a parameter , where is the normal burning rate of the grain, and is the port diameter. A dimensionless parameter is expressed as a function of , where , and is the mass flux through the port. The burning rate ratio is then expressed as
In Equation (26), is the Heaviside step function, and are parameters to be confirmed, and can be seen as a critical dimensionless mass flux. Mukunda and Paul recommend that , .
3. Experimental Apparatus and Analysis Procedure
In this paper, a large-aspect-ratio SRM with a free-standing grain with no inhibitor is developed. As the schema in Figure 3 shows, the gas in the motor will flow through two channels, one inside of the grain and the other outside of the grain between the grain and motor case. A conical nozzle is set up at the aft end of the motor. The throat radius is 7 mm. The ablation of the nozzle is neglected because of the short operation time. The grain burns out during the experiment with no interruption. The main geometric parameters of the motor are shown in Table 1.
A double-base propellant is adopted in the motor. The primary constituents are nitrocellulose, nitroglycerin, and dinitrotoluene, with the proportions 59.5%, 25%, and 8.8%, respectively. The density of the propellant is 1605 kg/m3. The normal burning rate of the propellant can be expressed as , where the unit of is MPa and the unit of is m/s. A chemical equilibrium calculation is carried out with the Chemical Equilibrium with Applications (CEA). The properties of the combustion gas obtained are shown in Table 2. These properties are assumed to be constant.
The pressure transducer installed at the head of the chamber collects data every millisecond, and the pressure-time trace at the head of the motor is obtained. As the computational model in this paper is quasi-steady and no ignition data is computed, the experimental data from 0.08 s after the ignition to the beginning of the tail-off is selected to fit with the calculation results. In the fitting procedures, the sum of squared errors (SSE) is used as the criterion.
4. Results and Discussion
The 1D model described in Section 2 is implemented to obtain the 1D flow field in the combustion chamber and the internal ballistics performance of the test motor. The two different burning rate models mentioned earlier are used. The geometric and physical parameters employed are from the SRM mentioned in Section 3.
To ensure that sufficient cells are used for the computation, three solutions obtained on different numbers of cells are compared. In each computation, the same time step 0.001 s is used. The Mukunda-Paul expression is used as the burning rate model and and are chosen, as Mukunda and Paul recommend. The finest mesh has 1600 cells. For each solution, pressure values at different times are compared to the values from the finest-mesh case at the same times. The largest relative differences are shown in Table 3. The relative differences obtained on different grid levels shows that the numerical error decreases as the number of cells increases. Since the difference between the solutions with the cell numbers 800 and 1600 is small, the cell number 800 is chosen for the subsequent computations. Similarly, the time step is chosen to be 0.001 s.
4.1. Identification of and in the Modified L-R Expression
The modified L-R expression can only be used to get the ideal data with appropriate and values. As previously discussed, a parameter identification procedure is needed to find out the values of and .
Since both and have an influence on the internal ballistics, a two-parameter optimization is required to determine the appropriate values. The sum of squared errors between experimental and computational curves, , is used as the criterion to evaluate the accuracy of the result. To avoid falling into local optimal traps, a relatively large range for and is chosen, and the distribution of SSE over the range is computed. Considering the values in References [3, 8, 10], the range of is set as to , and the range of is 20 to 200. The contour of SSE in this range is given in Figure 4. It can be seen that for a constant value of , the SSE decreases and then increases as the value of rises. Only one minimum value exists for a specific value. Figure 5 depicts the curve of minimum SSE values versus constant values. The figure shows that there is an overall minimum value among these minimum SSE values; this overall minimum SSE value in the range of and is found at the point where and . Considering the large range selected at first, this should be the true point where the SSE is minimum, not just a local minimum value. Using the values of and obtained, the curve is computed and shown in Figure 6. The computational results show good agreement with the experimental data.
As mentioned before, the value of can be computed theoretically by Equation (24). Limited by conditions, not all the parameters in Equation (18) are obtained in this paper, but a rough value, , is computed by combining the relevant parameters from Reference [6] and the propellant flame temperature . This value is fairly close to the value obtained by the optimization procedure.
4.2. Identification of and in the Mukunda-Paul Model
Using the value of and recommended by Mukunda and Paul [11] (, ), a numerical procedure is carried out to obtain the pressure-time trace of the combustion chamber. The results are shown in Figure 6. As we can see, the results do not fit the experimental data very well. However, before the tail-off stage, the relative error at each time is within 10%; the result is coincident with Mukunda and Paul’s work [11]. Considering that no parameter identification is needed to implement, the Mukunda and Paul model can still be a promising universal model for a coarse prediction of the internal ballistics performance of an SRM. After a trial-and-error procedure in which is maintained, a good fit is achieved when for the SRM in this paper.
The procedure to confirm the constants in the M-P expression is much simpler than that for the L-R expression. This is because the M-P expression is a universal model and the constants recommended are not far to their ideal values for most SRMs. Moreover, since the M-P expression is an explicit expression, it takes much less computation time to compute the total burning rate for each channel section at a time step using the M-P expression, than using the L-R expression.
4.3. Erosive Characteristics of the Motor
The plots in Figure 6 show that with proper erosion burning model, a good agreement with experimental results can be achieved by computation with the 1D model proposed in this paper. With the parameters from parameter identification procedures, both the L-R expression and the M-P expression can lead to good results. The pressure-time traces from the two models are similar to each other.
The deviations from experimental data at the tail-off stage are mainly due to the quasi-steady model used. In the calculation, the pressure changes immediately as the burning surface area changes and then influences the burning rate in turn. At the beginning of the tail-off stage, this makes the pressure to decrease more abruptly than that in the experimental data, and the regression of the burning surface goes slower. Thus, more propellant remains near the end of the operation, causing a slower pressure decrease at last.
The curve with erosion ignored is computed, with the normal burning rate expression, , as the burning rate model. With the erosive effect ignored, the pressure decreases approximately linearly and is lower than the pressure under the erosive condition. The deviation between these two curves is large at first and gets smaller. This indicates that the erosive burning is severe at the beginning of the motor operation and abates over time. The tail-off stage starts much later when the erosive burning is ignored. Once the tail-off begins, the pressure drops very quickly. This is because under this condition, the burning rate is dependent only on pressure and just a small difference exists along the grain. Thus, the propellant at different axial positions burns out almost at the same time.
The local burning rate ratios () in the inner channel at five different times are shown in Figure 7. The two sets of plots are from solutions using the L-R expression, and the Mukunda-Paul expression, respectively. The velocity at the aft-end section of the combustion chamber is also plotted.
(a)
(b)
As Figure 7 shows, the aft-end velocity decreases due to the increasing of the port area before the tail-off stage. During the tail-off stage, the area of the aft end of the grain equals to the combustion case section area, so the aft-end velocity remains. The burning rate ratio results from the two different burning rate models are similar to each other. Generally, the burning rate ratio increases with the axial coordinate. Before the tail-off stage, as time goes on, the augmentation of the burning rate decreases. This is the result of the velocity decrease due to the increment of the port areas. When the propellant near the aft end burns out and the motor runs into the tail-off stage, the burning rate ratio of the grain begins to increase again (the plot for the time at 0.45 s). This is because the normal burning rate decreases as the pressure drops. As many studies have revealed [2, 15], the erosive burning tends to be more severe at the lower normal burning rate. The aft-end velocity is nearly 80 m/s at the tail-off stage. This is large enough to cause erosive burning for some propellants, especially at low pressure.
There also some differences between the two solutions using the modified L-R expression and the M-P expression. Because the modified L-R expression does not consider the threshold velocity, erosive burning can emerge when the mass flux or the velocity of cross flow is low. This will lead to overestimation of the erosive burning. As we can see in Figure 7, at each time, the erosion burning from the L-R expression results appears earlier in the channel than it does in the M-P expression results. Besides, in the plots from the L-R expression solution, erosive burning appears near the head end of the grain. Since the maximum value of the reverse velocity at the head-end section is only about 28 m/s, this should be unlikely. The Mukunda-Paul expression, on the other hand, gives more reliable results. The effect of the threshold mass flux is clearly reproduced by the M-P expression. And the erosive burning does not show up near the head-end section. The two curves of the time 0.45 s and 0.5 s shows that even if the velocity of the aft end of the grain is the same, the erosive burning differs. This is because the erosive is controlled by the dimensionless mass flux .
The burning surfaces of the inner and outer channels at different times are shown in Figure 8. The results are from the solution using the M-P expression. The “reference face” in the figure is actually the mid-face of the grain. Burning surfaces of the inner channel at different moments are under the reference face, while outer channel burning surfaces are basically above the reference face. It can be clearly seen that, due to the erosive burning, the propellant near the aft end of the grain has a greater burning rate and burns out earlier. The erosive effects of the two channels are not identical. The erosion is more severe in the outer channel, leading to a larger total burning rate. That is the reason why the burning surface of the outer channel crosses the reference face.
Figure 8 also shows that near the head-end section of the grain, the burning surface is almost parallel to the axis of the motor. This indicates that no erosive burning exists, and the burning rate in this area is controlled by the pressure only. Since the pressure descend along the axial direction is not very obvious, the burning rate in this region is nearly uniform.
5. Conclusions
A quasi-1D model for internal ballistics prediction of two-channel combustion chambers is developed under the assumption that the flow variables are identical at the aft end of the grain. Combining it with two different erosive burning rate models, the L-R expression, and the M-P expression, the internal ballistics of the SRM is computed.
Results show that the more severe the erosion is, the sooner the tail-off begins and the gentler it becomes. For the two-channel combustion chamber studied in this paper, the inner and outer channels have similar erosive characteristics, yet the erosive effect in the outer channel is a little more serious.
Before tail-off, the erosive burning becomes less serious as the burning surface regresses. However, during the tail-off stage, as the normal burning rate decreases due to the obvious pressure drop, the erosive burning tends to be more likely to occur.
The two different erosive burning models are compared. Results show that with proper constants, both the models can lead to pressure-time trace with good accuracy. However, since the threshold effect is not considered in the L-R expression, the erosive burning will be overestimated in the low-mass-flux region. On the contrary, the M-P expression takes account of the threshold effect and can provide a more accurate result. Besides, as the M-P model is dimensionless and universal, it shows great advantage when dealing with different propellants and SRMs. The M-P expression is also a low time-consumption expression that provides good computing performance.
The good agreement between the results from the computational and experimental results shows that with a proper erosive burning model (the M-P expression is recommended), the 1D model proposed in this paper is valid for flow field and internal ballistics prediction for the two-channel combustion chambers.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.