Unsteady Simulation of Ignition of Turbulent Reactive Swirling Flow of Novel Design of Solid-Fuel Ramjet Motor

This paper presents a numerical investigation of ignition and combustion stabilization of a novel design of a solid-fuel ramjet (SFRJ) motor with and without swirl flow. The proposed design includes two solid fuels, retaining the simple design of the classic SFRJ. Numerical simulations of unsteady, turbulent, reactive, and swirling flow coupled with solid-fuel pyrolysis have been performed using an in-house CFD solver. Experiments on SFRJ were conducted via a connected-pipe test facility to validate the developed code. Furthermore, the code was validated for chemical reactions, heat diffusion, and swirl flow by using benchmark test cases of shock-induced, semi-infinite plate, and dump combustor with swirl flow, respectively. Then, the proposed and classic designs were simulated for the same inflow conditions and configurations, and the results were analyzed and discussed. It is found that the mixing degree, reactant residence time, mass flux, ignition delay time, and regression rate improve when using the proposed design. Moreover, the proposed design reveals interesting observations of a new flame being created and merged with the main flame.


Introduction
It is well known that a solid-fuel ramjet (SFRJ) uses ambient air as an oxidizer during its flight. This oxidizer reacts inside the combustor with the solid fuel at high flight speed, since the SFRJ starts working at supersonic conditions. Because of this, the initiation of the reaction is extremely complex and difficult to attain. Normally, SFRJ does not start its operation from zero speed, but by using an auxiliary device to reach a certain speed, such as booster rocket or military airplane. This makes the ignition process challenging for designers and researchers, and a hotspot area of research. Several parameters affect the ignition process, including igniter type, solid-fuel ignitability, combustor design, and operating conditions. Solid-fuel ignition can be defined as a complex phenomenon that involves many physicochemical processes with a sequence of events that starts with energy transfer to the solid fuel by an external stimulus, then heating and subsequent decomposition of the solid phase, diffusion of vaporized gases near the surface, and, finally, subsurface, heterogeneous, and/or gas-phase reactions. Sustained ignition is achieved once the heat released from the chemical reactions overcomes the heat losses. Once the external stimulus is removed, a sustained combustion will follow the ignition, otherwise the ignition is incomplete.
The ignition delay can be expressed as the period from starting the external stimulus to achieving sustained ignition. In general, three characteristic time intervals control the ignition delay, the time required for heating, the mixing (diffusion plus convection), and the reaction. In the study of ignition, ignition delay is one of the most important parameters. Nevertheless, the instant of sustained ignition is not easily recognized precisely.
There are well-documented studies about SFRJ motors and other air-breathing engines [1][2][3][4][5][6]. Most of these studies focus on the internal design and geometry effects in the performance of SFRJ. Although SFRJ contains no moving parts, the internal design has great influence on its performance. Gany et al. [7] experimentally investigated how geometry effects the combustion characteristics of SFRJs. It was found that the local regression rate is closely related to the local convective heat flux, downstream conditions do not influence the regression pattern, and mean regression rate decreased with the increase of port diameter. Moreover, Gong et al. [8,9] pointed out that the star solid fuel demonstrated higher space-averaged regression rate than that of the traditional single-port cylindrical solid fuel. They conducted an experimental and numerical investigation on the combustion characteristics of the SFRJ with star solid fuel.
Swirl flow has a major impact on flame stabilization and performance of SFRJ, as reported by many studies [10][11][12][13][14][15][16][17][18][19]. William et al. [20] experimentally examined the SFRJ at low air mass flux with swirl intensity varying from 0 to 0.649. They reported that the regression rate of HTPB (Hydroxyl Terminated Polybutadiene) increased significantly when swirl flow was used. The study of solid-fuel ignition is complicated by the fact that the contribution of individual unit processes may shift with changes in solid-fuel composition, pressure, heating mode, hydrodynamic factors, and gaseous environment. The ignition of solid fuels and propellants has received considerable attention in recent years. However, it should be noted that the SFRJ ignition is different from that in rockets' solid propellant, since a combination of fuel and oxidizer is used by rockets' solid propellant while ambient air is used by SFRJ as oxidizer. Tahsini et al. [21][22][23] numerically studied the ignition in SFRJ and back-step flows, as well as the igniter jet dynamics in SFRJs. Furthermore, Tahsini [24] studied supersonic hydrogen/air cylindrical mixing layer in scramjets to explore the influence of inlet swirl on ignition delay time.
To determine the ignitability of solid fuel, a simple experimental method may be used with an arc image furnace to subject the solid-fuel surface to a radiant flux and then detecting the ignition delay. However, this method neglects the real working conditions and motor design, thus, it is not enough to fully describe the ignition process. Furthermore, using a laboratory to adequately simulate the ignition process or static test is the major difficulty surrounding SFRJ studies. This is due to the fact that for the entire ignition event, the time period is very short, and the region of major activity is extremely small. Accordingly, no broad experimental database is as yet available, because it is hard to probe and observe the ignition region. Therefore, the current paper studies the ignition process in SFRJ motor thoroughly, using numerical simulation.
This study aims to numerically explore the ignition process and combustion stabilization of SFRJ motor for different combustor's designs, with and without swirl flow. First, a novel design is proposed and tested for both swirl and non-swirl flow. Then, both classic and proposed designs are simulated using an in-house CFD code developed via FORTRAN language and parallel computing.

Fluid Domain
To include all velocity and Reynolds stress components in axisymmetric solvers for three-dimensional (3-D) flow field with angular symmetry (∂/∂θ = 0), the simplifying condition of angular symmetry is used. Thus, the conservation of energy, momentum, mass, and species equations that are represented by the 3-D Reynolds-average Navier-Stokes (RANS) equations with angular symmetry condition in integral form are given by [25][26][27]: where Ω is the control volume, t is the time, n = n x i + n y j is the unit normal vector in outward direction of the boundary surface s. The viscous flux vectors F v , G v , convective flux vectors F, G, conservative vector U, chemical reactions source term S, and axisymmetric source terms H, H ν , are: • ε = 1 for axisymmetric swirl flow in which x and y denote z and r in cylindrical coordinates system, respectively. • ε = 1 and w = 0 for axisymmetric flow in which x and y denote z and r in cylindrical coordinates system, respectively. • ε = 0 and w = 0 for two-dimensional plane flow in which x and y denote x and y in Cartesian coordinates system, respectively.
By using this technique, three solvers are combined in one code, which makes the developed code more flexible. In addition: where p, T, u, v, w, E and ρ represent the pressure, temperature, axial, radial, tangential velocities, total energy per unit mass, and total density, respectively; s ρ i , s ρE , s ρv , s ρu , s ρw , and s ρ , are the species source terms due to solid-fuel pyrolysis and chemical reaction at the first cell in the fluid domain attached with the solid fuel, energy, momentum, and mass source terms, respectively. And τ is the shear stress, N is the number of species, k is the thermal conductivity coefficient of the gases which equals to the summation of its laminar k l and turbulent k t = C p (µ t /Pr t ) parts, Pr t is the turbulent Prandtl number, ω i is the mass production rate of species i due to chemical reactions, c i = ρ i /ρ is the mass fraction, and µ is the total effective viscosity which equals to the summation of laminar µ l and turbulent µ t viscosities. The mixture's laminar viscosity is determined by: where X i , M i , and µ li are the mole fraction of species, molecular weight, and laminar viscosity, respectively. In the mixture and for species i, the mass diffusivity D i can be determined via and the turbulent and laminar Schmidt numbers are adopted as constants, S ct = Pr t = 0.8, S c = 0.5, respectively. The Equation (4) is used to calculate the total energy per unit volume and then the temperature is determined by the Newton iteration method.
where h 298 i is the heat of formation at reference temperature (298 K), and the specific heat at constant pressure C pi at high temperatures is determined by: where R u is the universal gas constant and the chemical kinetics package [28] is used to get the coefficients a ji (j = 1, .., 5). The equation of state for thermally perfect gases is: where the ratio of total specific

Finite Rate Model
In the present work, chemical reaction processes are estimated using finite rate model neglecting the turbulence-chemistry interaction. The reaction mechanism for I reactions can be expressed by and the rate of mass production of species N in I reactions is given bẏ where v ik and v ik are stoichiometric coefficients and for each reaction,ṙ b k andṙ f k are the backward and forward reaction rate constants, respectively; L M = 1 when there is the third-body (M) and L M = 0 elsewhere, C j is the third-body efficiency. To calculate reaction rate constants, the Arrhenius formula (Equation (9)) is used.
where T w is the wall temperature at solid-fluid interface, n k is the temperature exponent, E a is the activation energy, and A k is the pre-exponential factor.

Solid Domain
In the solid domain, to model the heat diffusion, Fourier's equation (Equation (10)) is used in which RANS equations (Equation (1)) are reduced by considering only one of the species (solid) and evaluating all velocities to zeros. Then, the Fourier equation is given by and where k sol , c sol and ρ sol denote the thermal conductivity of the solid material, specific heat capacity, and density, respectively. S sol is the source term due to the chemical reaction of the burning surface attached with the fluid domain, T sol is the solid domain temperature, and h p is the reaction energy per unit mass.

Solid-Fluid Interaction
In this work, the high-density polyethylene (HDPE) has been employed as solid fuel for both the proposed and the classic designs with the assumption that C 2 H 4 is the only pyrolysis product of the solid fuel due to difficulties of having a precise chemical kinetics model besides the very complex reactions inside the SFRJ. The reaction model of a 3-step gas-phase [29] is used to reduce the computing uncertainties, since the 10-step model gives nearly similar results as a 3-step model [29,30]. Table 1 reports the gas-phase chemical reactions kinetic of HDPE [31]. The zeroth-order (n k = 0) Arrhenius equation is employed to model the solid-fuel pyrolysis in the solid-gas interface, then the fuel freed by pyrolysis mass rate is calculated by: The subscript "sol" denotes the solid fuel and the characteristics of HDPE are as follows: white color, density ρ sol of 940 kg/m 3 , thermal conductivity of 0.38 W/(m.K), A sol of 8750 m/s , E a sol of 130 KJ/mol [32,33]. The pyrolysis rate of HDPE due to heat transfer from fluid domain into the solid surface before combustion (gasification step) and during combustion are treated as a source term affecting the fluid-domain conservation equations.

Boundary Conditions
The boundary conditions are reported in Table 2; more information can be found elsewhere [34]. To calculate the wall temperature T w , direct coupling method is employed, which uses the energy balance as:

Numerical Simulation Procedure
The aforementioned RANS equations have been solved by means of finite-volume method, cell-centered, multi-block with structured grids, and density-based approach. Whereas AUSMPW+ [39] is selected to discretize the inviscid terms and is enhanced by Van Albada [40] limiter function and the third-order upwind-biased MUSCL (Monotone Upstream centered Scheme for Conservation Laws) scheme [41], the second-order central differences are selected for the viscous terms, Menter's shear-stress transport (SST) turbulence model [42] for turbulence closure, and the LU-SGS implicit dual time-stepping algorithm [43] for temporal discretization. In order to calculate chemical source terms, the finite rate model is selected. However, the effect of turbulent fluctuating (turbulence-chemistry interaction) is neglected, and the flame is considered to be laminar flame, at which the Arrhenius expression is used to calculate the mean reaction rate. More details are available elsewhere [44]. The solid domain is modeled for heat diffusion by using the second-order central differences to discretize Fourier's equation and explicit backward Euler approach for temporal discretization [43]. The direct coupling method is implemented to couple the solid and fluid domains, besides adding the source terms into the fluid governing equations.

Code Structure
To model the multi-physics coupling of chemical kinetics, heat transfer, fluid mechanics, thermodynamics, and solid pyrolysis, FORTRAN language is used to develop a multi-physics coupling code along with using OpenMP (Open Multi-Processing) programming interface, which makes it capable of being used in parallel computers. The flowchart of general calculation procedure and code structure is illustrated in Figures 1 and 2, and can be summarized as: • Pre-processing module: The physical structure is generated, the mesh is created, then computational domain is built by the code.
To select the type of simulation, such as swirling/non-swirling, reacting/non-reacting, steady-state/unsteady, inviscid/turbulent, or symmetric/axisymmetric, and other parameters, control parameters can be adjusted. Material property parameters, inlet, and exit parameters can be defined as well.
• Numerical simulations module: In this step, the numerical calculations are performed with chemical reactions coupled with fluid mechanics simulations, steady/unsteady simulation of compressible flow field, fluid-solid coupling simulation, conduction heat transfer simulation in solids, and coupled heat transfer simulation.

•
Post-processing module: Finally, the results data are exported to be visualized and processed in other software.

The New Design of SFRJ
The classic design of the SFRJ motor is presented in Figure 3a in which the central oxidizer flows from the inlet to the outlet without contributing to the combustion process. In the proposed design, however, the central oxidizer flow region is occupied by a new solid fuel while keeping the simplicity of the overall design as shown in Figure 3b. The newly added solid fuel is in rod shape configuration with a circular cross-section area. Therefore, a solid-fuel surface area contacted with the fluid flow is enlarged, and then the heat transfer and regression rate can be increased as well. Furthermore, the combustion will take place on the inner surface of tubular solid fuel, alongside the outer surface of the rod solid fuel. Thus, during flight, the tubular grain inner diameter will increase and the rod diameter will decrease.

Calculations Procedure
In this work, numerical simulations are conducted on both classic and proposed designs with and without swirling flow to explore the ignition and combustion stability of the SFRJ motor. Hence, four different cases are simulated in which the proposed design is simulated with and without swirl (Cases 4 and 3, respectively) beside two cases of the classic design with and without swirl flow (Cases 2 and 1, respectively). The inlet flow conditions and geometries are the same for all cases except the swirl intensity. However, the swirl number equals 0.6 for swirl cases and equals zero for non-swirl cases. The mass flow rate is set to be 0.6 kg/s and the total temperature is fixed at 540 K. As presented in Figure 4, the aft-mixing chamber diameter is 0.074 m, grain length is 0.300 m, grain port inner diameter is 0.070 m, and the ramjet inlet inner diameter is 0.040 m. The exit nozzle diameter is 0.053 m, the nozzle throat diameter is 36 mm, and the diaphragm diameter is 0.062 m. These dimensions are based on a typical experimental test engine. The meshes are 2-D structured grids and multi-block with total mesh number 104,982 and 82,746 cells for proposed and classic designs, respectively. The height of the first cell from the wall is 10 −6 m, which gives y + ≈ 1. The rod is 0.020 m diameter and 0.303 m in length.
Near the shear layer regions and along the walls, the meshes are clustered to precisely resolve the flame stabilization point and the steep gradients (see Figure 4). In addition, the first cell height from the interface wall is set to be fine enough to ensure grid independence of the heat flux estimation and to accurately resolve the temperature gradient in the solid domain. Alongside the adopted meshes and to check the grid sensitivity, simulations were conducted on the coarsest and finest meshes. No significant differences between the finest and adopted meshes were observed for all the flow parameters.

Experiments Setup
The connected-pipe facility is used to conduct firing tests on the SFRJ motor as shown in Figures 5 and 6. The facility consists of five basic sections: ramjet (nozzle, aft-chamber, and combustor); ignition unit; the air supply system; acquisition and measurements equipment. Initially, the size is set to the required values and the solid fuel is positioned with appropriate dimensions (as described in Section 6.1). Then, the incoming air is passed through the ramjet with mass flow rate (∼1.5 kg/s), temperature (∼567 • C), and pressure (∼1.7 MPa). Flight conditions up to Ma = 2.5 at sea level can be simulated. To start the experiment, the combustion chamber is provided with high-temperature gases from reacting the double-base (DB) propellant. Subsequently, during burning time of about 20 s, combustion of HDPE is initiated, and sustained combustion can be achieved. Lastly, nitrogen is employed to quench the combustion. Moreover, the pressure and temperature measurement sensors are placed at four and two different positions, respectively, as depicted in Figure 5. The 3-D scanner apparatus is used to measure the local regression rates after burning ( Figure 6); detailed information can be found in [10].

Model Validation
The developed code was verified and validated in our previous work [45], and in this section, more cases are tested to check its accuracy and robustness for solid-fuel decomposition, swirling flow, combustion, and heat diffusion.

Experiments on SFRJ
To assess the developed solver for predicting the solid-fuel decomposition and regression rate, experiments were conducted on a SFRJ motor using a connected-pipe facility. Two cases were tested with and without swirl flow and compared with numerical simulation. For the non-swirl case, the geometry configuration is adjusted as 300 mm grain length, 35 mm inlet inner diameter, and 70 mm grain inner diameter. For swirl case, swirl number was 0.6, 40 mm inlet inner diameter, and grain length was 500 mm. The inflow conditions that simulate the flight condition of Ma = 2.0 at sea level were 540 K of total temperature and 0.3 kg/s of mass flow rate.
As shown in Figure 7, the static pressure of the inlet air equals the static pressures of the combustor and aft-chamber before ignition (<2 s). Thereafter, all pressures raised suddenly at the start of the ignition process and then declined for steady pressures during sustained combustion period to about 20 s and after quenching the combustion; pressures were dropped to the initial values (>22 s). After experiment, the 3-D scanner was used to detect the regression level (surface morphology) and the results are displayed in Figures 8b and 9b for non-swirl case and swirl case, respectively.
Moreover, local regression rates for experiments and numerical simulation are compared and presented in Figures 8a and 9a for non-swirl case and swirl case, respectively. It is shown that for regression rate, the results are slightly under-estimated, which might be due to ignoring the effect of the radiation heat transfer, which influences the heat transfer coefficient and then the estimated regression rate. Nonetheless, the estimated result reveals good agreement with the measurement. It could be concluded that the developed solver can accurately predict the combustion and the regression rate in the SFRJ.   Table 4 shows the average regression rate and characteristic velocities that were obtained from numerical and experimental studies. According to the deviation calculated in Table 4, it can be concluded that the accuracy of this code for engineering investigations is acceptable, and the in-house code that is used in this research has enough precision in predicting the combustion and flow characteristics in SFRJ.

Flow Field Characteristics
In the swirl flow condition, the combustion characteristics of SFRJ is mainly determined by the flow field characteristics. In this section, case1 is used for analyzing the flow field characteristic; the combustion and regression rate become steady at the time of 1.352 s, thus, all of the results that are presented in following figures are the state in 1.352 s. According to Figure 16, the temperature reaches 2800 K around the recirculation zone and the reattachment point. In the redevelopment zone, the reaction of gaseous fuel and incoming air occurs, where the temperature reaches 2400 K. It could be seen from Figures 16 and 17 that the interface of the region with high CO2 ratio and the region with high O2 ratio has the maximum temperature. As described in Figure 17b, the region with higher mass fraction of CO coincides with the high-temperature region. Thus, it is obvious that the (a) (b) Figure 9. Measured and simulated local regression rate along HDPE surface (a) 3D scanner measurements of HDPE grain after burning (b) for swirl case.

The Swirl Case
Since the developed code can predict swirl flow, this feature has to be validated with the benchmark case, which was done with the test case of Nejad et al. [46]. The flow variables (u, v, w, p, k, ...etc.) with respect to the circumferential coordinate are negligible, then the flow is axisymmetric in their experiments. As a result, we considered only 2-D axisymmetric coordinate in the present work to simulate the swirl flow.
The measurement data are available starting from step downstream at x/H = 0.38 (the step height is H = 25.3 mm). The Re = 1.25 × 10 5 based on the inlet tube diameter (D = 51 mm) and the centerline axial velocity (U 0 = 19.2 m/s) at upstream of sudden expansion is considered to be a reference velocity. The simulations for swirl and non-swirl flows have been done for S = 0 and S = 0.5 (see Figure 10). The swirl number S can be given by: where R is the inlet radius. The normalized axial velocity is compared with the measurement data, as depicted in Figure 11. It is shown that for non-swirling flow, the developed solver agrees very well with the experimental data. For S = 0.5, the normalized axial and tangential velocities profiles are presented in Figure 12a,b, respectively, and good agreement is achieved for the stream-wise velocity close to the inlet. However, the agreement decreased for the axial velocity from x/H = 6 on, and for the tangential velocity the solver agrees very well with the experimental data.

The Chemical Reaction Case
Furthermore, the code is developed to predict the reacting flow and its predictability is evaluated using a shock-induced combustion benchmark case of Lehr [47]. The 600 × 300 structured grids were used for computational domain. Eight reactions with seven species are employed in the chemistry model. The comparison of steady-state results with measurement data are presented in Figure 13a for temperature and pressure distributions along the stagnation line, and in Figure 13b for mass fractions of the species (O 2 , H 2 , and H 2 O) at which good agreement was achieved. It can be concluded that the combustion solver is capable of predicting the combustion phenomena accurately.

Heat Conduction
The heat diffusion along the solid fuel needs to be validated as well; to do so, a semi-infinite plate is tested. Therefore, a uniform mesh (20 × 150) for rectangular plate (0.002 × 0.15 m) is employed. At x = 0, two different conditions are used: a constant heat flux condition q w = 70 kW (case b) and a constant temperature T w = 900 K (case a); and other boundaries are adiabatic. The initial temperature is fixed at 300 K. The plate has density of 7840 kg/m 3 , a thermal conductivity of 49.8 W/(m·K), and specific heat of 465 J/(kg·K). The analytical solutions are calculated using Equations (14) and (15) for transient temperature distributions. Excellent agreements are observed when comparing numerical and analytical results (see Figure 14).
T(x, t) where α = k/ρc is the thermal diffusivity of the material and er f () is the error function.

Results and Discussion
The developed solver has been employed to simulate the cases presented in the previous sections and the results are demonstrated in this section. Four different cases were simulated to examine the ignition and combustion characteristics of SFRJ motor. The cases represent the classic and proposed designs of SFRJ with and without swirl flow. This is to explore the influence of swirl flow and combustor design on the ignition and combustion characteristics of SFRJ.

Ignition Mode
To initiate the combustion in the simulations, hot gases were added to the incoming air flow with a certain mass flow rate (0.4 kg/s) and both flows are directed into the combustion chamber. The DB propellant is used to generate the hot gases, at which only the combustion products of the DB propellant are considered in the simulation. Therefore, the gases contain N 2 , H 2 O, and CO 2 with 2500 K total temperature. This ignition procedure is adopted according to real conditions of the experimental facility, as described in Section 5. The incoming flow (air and gases) is continued for a period of time and then the DB propellant gases have stopped, and only incoming air continues flowing into the combustor. This period of time is used in order to make sure that sustained combustion is achieved.
For the classic design cases, the hot gases were continued for one second and 0.60 s and 0.35 s for the proposed design cases with and without swirl, respectively.

Ignition Delay Time
The ignition delay time starts from the time spent on heating (gasification), mixing, and finally chemical processes. These three periods of time represent the ignition delay time; however, mixing time delay is often ignored, since it is very small compared to the others. In this study, we clearly observed both heating and reaction time delays, as will be discussed in this section.
Starting from the fact that the reactants and product concentrations somehow represent the ignition and combustion in the combustion chamber, this paper adopts mass fractions and temperature of the reactants to discuss the ignition and combustion especially near the solid-fuel surface. Figure 15a,b illustrates the mass fractions of the reactants (oxygen and C 2 H 4 ) for the proposed design with swirl flow at the surface of the solid fuel. Herein, the total combustion process is divided into three time zones-before ignition (Zone A), after ignition (Zone B), and sustained combustion time (Zone C). As soon as the oxide thickness reaches zero, the ignition will start, and sustained combustion may be achieved. Therefore, as presented in Figure 15a, the heating starts from the beginning and continues until around 0.14 s, at which the concentration of O 2 decreases gradually with time from 0.01 s, 0.1 s, 0.11 s, 0.12 s, and 0.13 s because of the HDPE pyrolysis with time, and thus the gas-phase of HDPE was raised; hence, this time zone is named the gasification step. Furthermore, oxygen concentration increases suddenly between 0.13 s and 0.14 s. This may indicate that the ignition may start during this period of time (0.13-0.14 s) and then continue until it reaches sustained combustion as the concentration of oxygen is zero and a fuel-rich region is achieved at the region close to the solid-fuel surface from 0.15 s until the end of combustion.
It is well known that in such configurations, the ignition will start at the recirculation zone, at which the pyrolyzed fuel gas and oxidant have higher residence time and good mixing degree. This supports our observation about the ignition being started between 0.13 s and 0.14 s; the mass fraction of O 2 at 0.14 s equals zero in the recirculation zone (<75 mm), as illustrated by the red curve of 0.14 s. This clearly confirms that the ignition has already started before 0.14 s. Note that due to higher heat flux at the reattachment point, which is located at the recirculation region's end, the temperatures of the gas close to the solid fuel are higher than that in the recirculation region, thus the vaporization rate of the solid fuel is high.
In contrast to oxygen, however, C 2 H 4 concentration is zero at the beginning at the grain surface and raises depending on the temperature, since it is a product of solid-fuel decomposition (see Figure 15b). The same trend as oxygen is observed: at first, C 2 H 4 increased gradually from the beginning until 0.13 s, and then suddenly increased at 0.14 s. This is due to a higher increase in temperature near the surface because of combustion initiation, as displayed in Figure 16. In this figure, the temperature at the grain surface is presented at different time intervals during the numerical simulations. The grain surface is heated first and its temperature is raised gradually as time progresses, until it reaches the gasification temperature of HDPE in which HDPE is pyrolyzed (Zone A), with a narrow peak temperature zone at the reattachment point. Thus, in this case, the time before ignition (Zone A) ends at 0.14 s, at which point the combustion is initiated. In this case, after starting the ignition, the hot gases continue to flow with the incoming air to avoid the combustion blow-out until 0.6 s, then only incoming air continues after stopping the hot gases. This time zone starts from 0.15 s until 0.60 s, which is herein named Zone B. In this zone, the oxygen concentration and temperature profiles are constant, and C 2 H 4 concentration showed steady behavior. As a result of stopping the ignition gases, the combustion fluctuates at 0.60 s and then sustained combustion is observed with almost constant oxygen concentration and temperature; the same trend is observed for C 2 H 4 concentration. Based on the fact that the ignition is complete if sustained combustion follows the ignition process, this case presents sustained combustion after ignition and then ignition is complete. After sustained combustion was achieved, the solver stopped at 0.67 s, and this time period is named time Zone C, which starts by ending the ignition gases, and ends at the stopping the simulation.
The proposed design contains two solid fuels-one is a tubular grain and the other one is a rod grain located at the center of the combustor, as mentioned in the previous sections. The ignition process is also discussed along the surface of the rod solid fuel. Figure 17a shows mass fractions of oxygen along the surface of the rod grain during the whole simulation process. At the first interval time (0.01 s), the mass fraction is about 0.15%, which is less than the oxygen mass fraction in the air. This is due to the ignition gases, which reduce oxygen concentration at this region. Then, the oxygen mass fraction decreases gradually as time progresses until 0.13 s at the end of the rod solid fuel (<300 mm). The concentration drops suddenly between 0.13 s and 0.14 s especially at 20 mm from the inlet, which means the combustion generated in the recirculation zone mentioned above has contributed to the ignition process near the rod solid fuel by adding high heat flux to the rod solid fuel, as presented in Figure 18. This figure clearly shows a temperature increase between 0.13 s and 0.14 s (red curve), and after that, the temperature is suddenly increased to the maximum at 0.15 s and near the fuel end (<300 mm). This indicates that the combustion is initiated starting from the end of the rod solid fuel and the affected area is increased towards the inlet as demonstrated in Figure 17a at which the mass fraction of oxygen is zero at 0.15 s between 300-280 mm, 0.16 s between 300-180 mm and so on until it reaches 10 mm. The mass fraction of C 2 H 4 illustrates similar behavior as presented in Figure 17b, in which the rod solid-fuel pyrolysis starts from the combustor's end and grows towards the inlet as time progresses. This is due to the recirculation zone formed at end of the rod solid fuel. A rich-fuel region is created after stopping the ignition gases, and sustained combustion is achieved. The same trend is also observed for the proposed design without swirl flow, as displayed in Figures 19-21 for the tubular and rod grains, respectively. In this case, the ignition gases are continued until 0.35 s and the simulation was stopped at 0.44 s after obtaining steady-state combustion. The time zones for this case are clearly demonstrated in Figure 22a.

The Regression Rate
The instantaneous regression rates of the solid fuel for the proposed and classic designs with and without swirl at all time zones are displayed in Figures 22a and 23a, respectively. To describe the impact of the proposed design, the behavior of sustained combustion and ignition of different designs is discussed. The major concern is whether the proposed design can positively affect ignition and combustion of SFRJ for non-swirl and swirl flows. For all cases presented in Figures 22 and 23, the regression rates were small but afterwards increased gradually, and the amount of fuel released converts to pyrolyzed fuel vapor. This period of heating is named Zone A, at which the behavior is different from each case depending on the flow type. Swirl cases, for both designs, revealed a higher regression rate, especially in the recirculation region, which indicates the swirl enhances the solid-fuel pyrolysis. Moreover, the heating process (Zone A) of the proposed design without swirl ends at 0.19 s, since the regression rate increased suddenly, at which point the ignition starts, between 0.18 s and 0.19 s, as shown in Figure 22a. This means the swirl decreases the ignition delay time for the proposed design, since swirl cases presented low time delay. However, the presence of swirl enhances the degree of mixing and the residence time of the reactants due to the decrease of axial velocity especially close to the centerline, then enhances the swirl velocity near the wall. This significantly increases the heat transferred into the grain surface, and thus decreases the delay time of the ignition.
For time Zone B, the ignition gases are continued after combustion is initiated and completed, to avoid combustion blow-out. The time period for the classic design cases was 1 s and for proposed designs are 0.35 s and 0.60 s for non-swirl and swirl flow, respectively. All cases demonstrated a steady tendency of regression rates in Zone B with a higher rate for the proposed design with swirl (see Figure 22b). The narrow peak appeared in swirl cases caused by the central recirculation zone, which is located at the combustor's inlet, then increased the radial velocity gradient. Therefore, the corner recirculation length decreases with swirl flow then regression rate, and local heat flux is increased gradually in the recirculation region until it reaches the maximum at the reattachment point; thereafter, it decreases progressively alongside the reattachment point.
As soon as ignition gases were stopped, the regression rates dropped suddenly, and then steady profiles were achieved due to the loss of incoming mass flow rates. Therefore, from these points onwards, the combustion is steady, and this region is named time Zone C, as can be seen in Figures 22 and 23. However, an additional maximum peak of regression rate close to the combustor's end is observed for non-swirl case of the proposed design (see Figure 22a-Zone C). This is due to the appearance of a second flame at this region and then a high rate of heat transferred into a tubular solid-fuel surface. In addition, combustion stabilization can be observed from this figure as well by noticing the fluctuation of regression rates at the time of stopping ignition gases. Swirl cases for both designs showed less fluctuation than non-swirl cases; alongside this, regression rates at Zone C for swirl cases have shown a steady manner, unlike non-swirl cases, which reveals that the swirl enhances the combustion stabilization. Furthermore, insignificant differences are observed between the proposed and classic designs regarding combustion stabilization.
The fuel released from the rod grain created the second diffusion flame near to its surface, which pushes the first flame close to the tubular grain. For the proposed design with swirl flow, however, this manner was not observed, since the two flames were merged inside the combustion chamber, unlike the non-swirl case, which merged in the aft-burning chamber. This could be obviously seen in temperature contours visualized in Figure 24.
Additionally, since the inlet and port cross-section areas decreased in the proposed design, the mass flux is increased; then, the corner recirculation length is decreased, and local regression rate is increased. This can be noticed if one compares Figure 22b with Figure 23b and Figure 22a with Figure 23a. Hence, the proposed design gives higher mass flux and then improves the regression rate for the same configuration as the classic SFRJ. Finally, one can conclude that the proposed design reveals higher regression rates than the classic design for both swirl and non-swirl flows. Furthermore, near the outer surface of rod solid fuels, the regression rate propagation is examined as well. It is observable from Figure 25a that the solid fuel starts the pyrolysis process near the combustor's end (between 0.15-0.2 m) for non-swirl cases, and from the combustor's inlet (between 0-0.05 m) and outlet for swirl cases, Figure 25b, in Zone A, due to the central recirculation zone. For both cases, regression rate exhibits much higher values than tubular grain for proposed design, and then higher average regression rate than the classic SFRJ cases.
Furthermore, the proposed design contour visualization of temperature is demonstrated in Figure 26, at which the reacting, turbulent, unsteady, and swirl flow in the combustion chamber is illustrated for the whole simulation process. It can be seen that ignition starts between 0.13 s and 0.14 s, sustained combustion using ignition gases happens in 0.15-0.60 s, and finally sustained combustion is reached at the simulation's end.
Moreover, in sudden expansion combustors, such as an SFRJ combustor, the recirculation zone is used to provide a stabilization mechanism to the diffusion flame. Therefore, the feature of recirculation flow in the combustor has a major influence on the reaction process, starting from the ignition process until flame generation. The reattachment point follows the recirculation zone, at which point the redevelopment zone starts and then within the redeveloped boundary layer a turbulent diffusion flame is formed (diffusion-controlled flame) (see Figure 27). For this reason, in this paper, the axial velocities close to the tubular solid-fuel inner surface for both proposed and classic designs are displayed in Figure 28a,b, respectively. Comparing the two results, it can be seen that there was a significant difference between the two designs, in which the proposed design reduces the recirculation zone length to about 28 mm (close-up view in Figure 28a) from 70 mm (close-up view in Figure 28b) for classic design in time Zone C. Finally, Figure 29 shows the central recirculation region generated due to the presence of swirl flow for the proposed design.

Conclusions
In this paper, the ignition and combustion stability of a novel design of a SFRJ motor were investigated by means of a developed in-house CFD code. Two solid fuels were employed by the proposed design with the same simplicity of the classic design. Moreover, a connected-pipe facility was used to conduct firing tests on a SFRJ without and with swirling flow. To perform reactive, unsteady, turbulent, and swirling flow simulations coupled with solid domain, an in-house CFD code was developed. The solver has been assessed for predicting combustion by shock-induced benchmark case, swirl flow by dump combustor, benchmark case, heat diffusion by semi-infinite plate, and solid-fuel decomposition by using experimental data of the SFRJ. Good agreement was observed. In addition, the code has been used to perform numerical simulations on the proposed and classic designs, with and without swirl flow. It is shown that the proposed design increases the mixing degree, reactant residence time, mass flux, and regression rate, and decreases the ignition delay time. For both designs, the presence of swirl increased the regression rate, reactant residence time, and mixing degree, and decreased the ignition delay time. Furthermore, in the proposed design, a new flame was generated near the rod solid-fuel surface and merged with the first flame inside the combustor for the swirl case, and in the aft-burning chamber for non-swirl case.

Conflicts of Interest:
The authors declare no conflict of interest.