Based on the aforementioned research works, a typical combustion chamber model of the powder-fuel ramjet is created (Figure 1). The oxidizer, air, enters the combustion chamber through two separated air inlets mounted on the combustion chamber along the streamwise direction, shown as air inlet 1 and air inlet 2 in the figure. The air supplied from air inlet 1 makes the air-to-powder ratio close to the equivalent ratio and then meets the requirements of engine startup and ignition (Auerswald et al., 2020). Air from air inlet 2 increases the air-to-powder ratio to improve the specific impulse of the engine and combustion efficiency. There are some exhaust inlets arranged on the head of the combustion chamber. The high-temperature and high-speed exhaust is injected into the combustion through these exhaust inlets. The powder fuel is injected from the central inlet on the head of the chamber by the fluidization gas. The working medium of the fluidization process is air. The fluidizing powder mixes with the exhaust from the exhaust inlet, air from air inlet 1 and air inlet 2 in the combustion chamber.

From the published articles, some parameters of a solid rocket ramjet are given in Table 1. In the table, the geometric configuration of the experimental combustor and test conditions are available. Figure 2 shows the related geometric parameters of the combustion chamber. The diameter of the combustion chamber is 160 mm (d) with a length of 1,560 mm (l). The simplified engine has six main exhaust nozzles, with a diameter of 16 mm (d ₂). They are distributed circumferentially at an angle of 60 degrees. The diameter of the powder inlet is 21 mm (d ₁). The lengths of the two rectangular air inlets on the side wall are 82 mm (l ₁) and 52 mm (l ₂), respectively. From the top view, the angle between the two air inlets is 90 degrees. To study the influence of the structural shape in the head and exhaust inlet of the combustion chamber on the mixing and engine performance, dome height 1, dome height 2, and the distance between the high-temperature exhaust inlet center and the fuel inlet center are fixed at 110 mm (h ₁), 330 mm (h ₂), and 37 mm (l ₃), respectively.

Based on the traditional design with a cylindrical head and six exhaust inlets type (case A) Chuang et al. (1989), two modified designs varied in head shapes and three exhaust inlet patterns were built aiming at improving the mixing efficiency. The head shapes of case B and case C are designed as a round head type and a coned head type, respectively, as shown in Figure 3 (case B) and (case C). In Figure 4, the high-temperature exhaust inlet is designed as four-hole type in case A1 and eight-hole type in case A2. Case A3 is designed as six inclined holes, with an inclination angle of 30 degrees. The details of the six configurations are described as follows:

Case A: basic model, cylindrical head, six exhaust inlets. Case B: round head with spherical diameter of 243 mm, six exhaust inlets. Case C: coned head with the inclined contracted angle of 22 degrees, six exhaust inlets

Case A1: cylindrical head, four exhaust inlets

Case A2: cylindrical head, eight exhaust inlets

Case A3: cylindrical head, six inclined exhaust inlets with an inclination angle of 30 degrees

In the calculations, three turbulence models are tested and compared. Commercial CFD software ANSYS Fluent 19.0 Liu et al. (2020b) is used for simulation. It uses parallel algorithm to solve the three-dimensional Navier–Stokes equation in a structured multi-block grid system. In order to simulate the compressible flow regime, a density-based solver scheme Roe (2003) is applied.

Murty and Chakraborty (2012) used the k-ε standard model in the study of mixing effects and obtained pretty accurate results. The transport equations of k and ε are shown as following two equations: ∂ ∂ t ( ρ k ) + ∂ ∂ x i ( ρ k u i ) = ∂ ∂ x j [ ( μ + μ t σ k ) ∂ k ∂ x j ] + G k + G b − ρ ε − Y M + S k and ∂ ∂ t ( ρ ε ) + ∂ ∂ x i ( ρ ε u i ) = ∂ ∂ x j [ ( μ + μ t σ ε ) ∂ ε ∂ x j ] + C 1 ε ε k ( G k + C 3 ε G b ) − C 2 ε ρ ε 2 k + S ε .

In the two equations, μ _t is written as follows: μ t = ρ C μ k 2 ε .

The species transport equation is also used in this article, and this conservation equation takes the following general form: ∂ ∂ t ( ρ Y i ) + ∇ ⋅ ( ρ v ¯ Y i ) = − ∇ ⋅ J ¯ i + R i + S i , where R _i is the net rate of production of species i by chemical reaction and S _i is the rate of creation by addition from the dispersed phase plus any user-defined sources.

The k-ε standard turbulence model has been used in this work. The pressure values predicated by the turbulence model and the pressure measured by the experiments in the combustion chamber of the ramjet are compared in Figure 5. Obviously, the predicted results by the k-ε standard model have good agreement with the trend with the tested result and the ratio of pressure to total pressure at each pressure measurement locations along the axis in different operation conditions is nearly overlapped. Since the three cross-sectional areas show a trend of increasing and then decreasing, the pressure ratio shows a trend of decrease and then increase. In addition, the results show that the model has good stability and the obtained value has acceptable deviations in a certain range of fluctuations. It should be pointed out that there are some inevitable differences between simulation results and experimental results. For example, the pressure test region has relatively large disturbance by the instruments. However, the fluctuation can be ignored in the simulation process, and the pressure change is not obvious. Therefore, the k-ε standard model has shown enough validations for later calculations.

In the present work, the mixing effect of high-temperature exhaust and the powder is replaced by the mixing of high-temperature exhaust (H₂O and CO₂) with air ejecting from the central inlet. Table 2 gives the related parameters of air and high-temperature exhaust and the boundary conditions of the work. The assumptions in the calculations are listed as follows.

1 The fluidized gas of powder is assumed as air.

In order to reduce the calculation effort, the computational domain is half of the entire domain with the symmetry boundary conditions. The structured meshes are generated to improve the calculation accuracy. The meshes in the head region are dense because the flow structures and mixing effect in this region are complex. The structured meshes in the whole domain and some typical regions are displayed in Figure 6.

The computational domain and grid are generated using commercial software ANSYS ICEM 19.0. A multiblock grid approach is used (as shown in Figure 6). To ensure the accuracy of the numerical simulation, cells are densely close to the wall surface, the high-temperature exhaust inlet, and the powder inlet. The height of the first cell row is made a value of 0.1 mm near the wall, and the total number of cells is 0.9 million.

In order to eliminate the influence of the grid sizes on the calculation results, three different mesh systems are built for the grid independence study, the grids are coarse, medium, and fine, and the number of grids is 0.26, 0.91, and 1.48 million, respectively. For different grid regimes, the concentration distribution of carbon dioxide on the section of the combustion chamber is compared. Figure 7 shows the concentration distribution of carbon dioxide at each interface of the combustor displayed by the three different grid regimes. The y-z cross section at x/l = 0.25 is selected from the location of the side wall air inlet mounted the combustion chamber, and the concentration value shows great fluctuation. Therefore, the carbon dioxide concentration in the cross section of the combustion chamber is significantly reduced after x/l = 0.25. From Figure 7, the calculated carbon dioxide concentration on the cross section by three grid systems is basically overlapping. The medium grid and fine grid are closer to each other at the separation point. In order to balance the calculation efficiency and accuracy, medium grid is used for subsequent calculation. The average value of y plus calculated on the grid wall is 22, which is set based on the requirement of the standard k-ε turbulence model.

In order to figure out the flow field structures inside the combustion chamber in each design, the streamline distributions of the symmetry plane under various conditions are presented. Figure 8 shows the contours of the mole fraction (CO₂) and streamline distributions on the x-z central sections of the combustion chamber in each case. The high-temperature exhaust at the inlet of the combustion chamber is water vapor and carbon dioxide, and the central fluidizing gas is air. It can be seen that the high-temperature exhaust mixes with the supplied air at air inlet 1 and air inlet 2 in each case, which forms a recirculation zone inside the head of the combustion chamber. In different cases, the scales of the recirculation zone are different. The small recirculation zones cause the powder to gather together and have a negative effect on the operation of the combustion chamber. Moreover, the powder cannot burn totally in the combustion chamber because of the short residence time so that the production of the powder combustion gathers inside the head of the combustion chamber. Compared with case A, case B and case C have similar symmetrical recirculation zones and the scale of the recirculation zones is relatively large, which can cover the powder injection region. In case B and case C, high-temperature exhaust and air are mixed to form a recirculation zone with a symmetrical structure in the head, which increases the stay time of the powder and can easily achieve fully burned. In case A, some metal powder is wrapped in the recirculation zone and cannot be transported to the core combustion region for full combustion. Compared with other cases, the high-temperature exhaust in case B and case C surrounds the powder particles closely to peel off the surface oxide layer and liquid film of the powder particles continuously, which helps the powder particles to burn completely. In case A2 and case A3, it can be seen that increasing the number of high-temperature exhaust inlets is helpful to the fluidization effect of the powder and the formation and extension of the recirculation zone. In case A3, it can be clearly seen that the heightened effect of the recirculation zone is obvious, and the scouring effect of high-temperature exhaust on the core reaction region is improved. When the exhaust and powder fuel are injected at high temperature with a certain angle, the mixing and combustion efficiency of powder and exhaust are stronger than other cases.

Contours of relative density of the mixing exhaust/air on two y-z sections (x/l = 0.125 and x = 0.25) for all the cases are displayed in Figure 9. The relative density is the ratio of the total density to the air density. relative density = total   density densit y air .

The high relative density is found on the middle of the cross section x/l = 0.125, and it also distributed two sides in the section of x/l = 0.25 for all cases. When the relative density of the certain area is close to 1, it indicates that the density of the area is close to the air. It also means that the mixing degree is better. And the relative density gradient describes the uniformity of the mixing field. From Figure 9A, the high relative density is found in small area inside cylindrical head and round head at the section of x/l = 0.25. It has disadvantages for the powder mixing with the gases around in the region. On the other side, the relative density distribution of the coned head is uniform with a good mixing efficiency. In addition, a better relative density is obtained in the case with the change of exhaust inlet patterns to the inclined six holes type. In the section of x/l = 0.125, a large region with relative high density is found in the coned head shape, which indicates that this design provides a better mixing performance than the other cases.

Pressure distributions on the x-y central sections and six y-z sections along the streamwise direction for all the cases are shown in Figure 10. A region with high static pressure is formed at the region of the two air side inlets. The high-pressure regions in case A and case B are larger than those in case C, which prevents the powder injection and mixing in the combustion chamber. The pressure loss is also shown in this figure, and the coned head type obtains the smallest pressure loss penalty. In case A1 and case A2, there is an obvious decrease in pressure inside the head region, but it is improved in case A3. Therefore, the inclined six holes type, that is, case A3, provides the best powder injection and mixing. Compared with case A1 and case A2, the pressure distribution in case A3 is more uniform.

In 1992, Lomkov and Kopchenov (Kopchenov and Lomkov, 1992) proposed the definition of a mixing degree in a certain section is defined as follows: D = ∬ ρ u ( c − c ¯ ) 2 d A c ¯ 2 ∬ ρ u d A .

In the formula, D , u , c , and ρ represent component mixing degree, flow velocity, mass concentration, and local density, respectively. A represents area, c ¯ represents the weighted average mass concentration on the calculated section, and D   =   1 represents no fuel injection. The mixing efficiency is very low. D =   0   means that the components are fully and uniformly mixed. In order to make the mixing efficiency proportional to the mixing degree, the next equation is defined as follows: η = 1 − D .

In this expression, η is the mixing efficiency and D is the mixing degree.

Figure 11 shows the comparisons of mixing efficiency distributions for all cases. Each datum is obtained by averaging all the results on the y-z sections along the streamwise direction. The side wall air inlet 1 and air inlet 2 cause additional flow mixing in the combustion chamber. The mixing effect is not obvious when the flow is approaching in the combustion chamber. Therefore, the starting point of the distribution in the combustion chamber is set at air inlet 2, and all the six distributions curves obtain the same trend. From Figure 11, the mixing efficiency of the process increases continuously in the section of the combustion chamber from air inlet 2. The final value is equal to 1, which means that the mixing efficiency in the combustion chamber gradually reaches the maximum value. However, in the entire flow field, the mixing efficiency of the coned type (case C) combustor shows an advantage at the starting region of the inlet 2 section. At the starting point, the mixing efficiency of the round type (case B) combustor is always lower than the mixing degree of other two cases. The curve of the round head type has presented a rapid increase in the considered range of x direction. This shows that the mixing efficiency of the round head is better than that of the cylindrical type. It can be seen that change in the number and angle of high-temperature exhaust inlets can have great influence on the mixing efficiency, when cases A1, A2, and A3 are compared. It is not obvious that by changing the head structures and the exhaust inlet patterns in the combustion chamber can directly affect the mixing efficiency. In addition, the relationship between pressure loss and mixing efficiency needs to be considered together to determine the mixing performance.

Figure 12 shows correlation between pressure recovery coefficients and mixing efficiency on the y-z sections along the streamwise direction. It can be seen from Figure 12 1) that the traditional cylindrical head of the combustion chamber gives more total pressure loss under the same mixing efficiency. At this time, the coned head type has improved this problem, and it always has high mixing efficiency than other types. This also indicates that the coned head type has the lowest pressure loss under the same mixing efficiency. It can be seen from Figure 12 2) that the three curves of mixing efficiency of case A1, case A2, and case A1 are nearly overlapping. When the head shapes are the round head or coned head, the total pressure loss and mixing efficiency are improved compared with the cylindrical head. In addition, from Figure 12A, it can be found that the round head is slightly inferior to the coned head. It can be seen from Figure 12 1) and Figure 12 2) that as long as the head shape is improved, there is obvious advantages in total pressure loss and mixing efficiency.

The pressure loss in the mixing process is an important evaluation factor for the structure design in a combustion chamber. In order to further compare the pressure loss of several cases and evaluate the engine performance, the pressure recovery coefficient in the combustion chamber is calculated.

The total pressure recovery coefficient is expressed as follows: σ = P ¯ 0 o u t P ¯ 0 u n d i s t u r b e d .

In the formula, σ represents the total pressure recovery coefficient, P ¯ 0 u n d i s t u r b e d represents the average total pressure of the undisturbed air flow, and P ¯ 0 o u t represents the average total pressure of the outlet.

The total pressure recovery coefficient of the inlet is defined as: the ratio of the average total pressure of the inlet to the average total pressure of the undisturbed part. It is used to measure the pressure loss of the airflow inside the engine.

The average total pressure (Li et al. 2017) is defined as follows: P ¯ 0 = ∬ P 0 ρ u d A ∬ ρ u d A .

In the formula, P ¯ 0 represents the average total pressure; P 0 represents the total pressure, ρ represents the density, u represents the speed along the x axis, and A represents the cross-sectional area.

According to the total inlet pressure recovery coefficient, the pressure recovery coefficient of the combustion chamber is defined. The pressure recovery coefficient is defined as follows: λ = P ¯ 0 s e c t i o n P ¯ 0 i n , where λ represents the pressure recovery coefficient of each section of the combustion chamber, P ¯ 0 s e c t i o n represents the average total pressure of each section, and P ¯ 0 i n represents the average total pressure at the inlet.

Figure 13 shows the pressure recovery coefficients distributions for all cases along the axial section. The six curves have the same trend, and the section pressure recovery coefficient shows a gradually decreased trend along the axial direction. It can be seen that the pressure recovery coefficient of case A is always lower than that of other cases in Figure 13A. Even compared with case A1, case A2, and case A3, the pressure recovery coefficient of case A is the smallest. It also indicates that the pressure loss of the combustion chamber with the modified head shapes and exhaust inlet patterns is smaller than that of case A. Among them, case C has the best performance with the largest pressure recovery coefficient compared with other configurations. Therefore, it is more recommended to use a round head or coned head type to improve mixing effect of powder and high temperature exhaust. When the number of high temperature exhaust inlets is changed, the effect of improving the pressure recovery coefficient is not too obvious. It can be clearly seen from Figure 13 that the cross-sectional pressure recovery coefficient in the final stage has increased. Because the cross section is located in front of the engine throat, the pressure recovery coefficient in this region is increased with the pressure increasing.