- Email: [email protected]
Effects of thermal wall conditions on rotating detonation
Computers and Fluids 140 (2016) 59–71
Contents lists available at ScienceDirect
Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid
Effects of thermal wall conditions on rotating detonation Yuhui Wang a,∗, Jianping Wang b, Wenyou Qiao a a b
Research Center of Combustion Aerodynamics, Southwest University of Science and Technology, Mianyang 621010, China College of Engineering, Peking University, Beijing 100871, China
a r t i c l e
i n f o
Article history: Received 19 March 2016 Revised 25 June 2016 Accepted 8 September 2016 Available online 8 September 2016 Keywords: Rotating detonation Boundary oblique detonation Numerical simulation Thermal wall conditions
a b s t r a c t Rotating detonation in an annulus full of stoichiometric hydrogen and air is numerically studied to figure out the effects of thermal wall conditions on detonation. A transient density-based solver with implicit formulation, a laminar finite-rate model with one step reaction, and a standard k-epsilon model are employed. It is found that after some time from ignition, boundary oblique detonation near the walls occurs and is included into the Y-shaped rotating detonation formed from the detonation front and transverse wave. The transverse wave is comprised of detonation attached to the rotating detonation front and shock. Boundary oblique detonation occurs with high wall temperature conditions, but it does not occur with adiabatic wall conditions. It is the rapid reaction caused by the high wall temperature that induces boundary oblique detonation. The length of boundary oblique detonation is increased with increasing wall temperature and time. The average velocity of rotating detonation for higher wall temperature is higher. Why continuous propagation of rotating detonation in engines without cooling is not broken by deflagration due to hot walls is explained. Numerical results show that transverse waves produce sub-peaks of pressure traces in the experiments. © 2016 Elsevier Ltd. All rights reserved.
1. Introduction Detonation has lower entropy production than deflagration, which improves fuel efficiency. Rotating detonation engines (RDEs), which have 20% higher thermal efficiency than deflagration engines [1], have been studied theoretically, numerically and experimentally. An RDE combustor is usually annular and detonation in it is non-premixed combustion in order to avoid backfire. Though the high-pressure region behind the rotating detonation wave (RDW) blocks a part of the inlet, the reactants can still flow into the combustor through the other parts of the inlet to maintain continuous detonation. There may be one or more rotating detonation waves rotating around the combustor channel. Since the detonation wave is almost hypersonic, when the inflow has a lower velocity than the detonation wave it always stays near the combustor inlet where the mixing of fuel and oxidant is insufficient. When the inflow is hypersonic and faster than the detonation wave, a convergent inlet can be used to stabilize the rotating detonation wave axially. Thus, the RDEs apply to air-breathing engines especially. More detailed operating principle can be found in the references [2–5]. OH∗ chemiluminescence imaging was used to capture the RDW structures of different mass flow rates, equivalence ratios, air injection areas and so on [6]. The radial RDW front and the tra-
∗
Corresponding author. Fax: +008608162419256. E-mail address: [email protected] (Y. Wang).
http://dx.doi.org/10.1016/j.compfluid.2016.09.008 0045-7930/© 2016 Elsevier Ltd. All rights reserved.
jectories of triple points in the two-dimensional curved channels were observed using multi-frame short-time open-shutter photography [7]. In the experiments the specific impulse for the RDE of a hydrogen and air mixture with a choked converging area aerospike nozzle was about 430 0–560 0 s for different equivalence ratios of the mass flow rate 1.14 kg/s [8]. Since detonation has higher temperature and pressure than deflagration and the RDW is moving all the time, the cooling of RDEs is a great challenge. The mean heat fluxes to the combustor walls were almost the same for detonation and deflagration, but the peak heat flux during an unsteady operation could be two to three times higher [9]. The pressure traces of the RDW in an RDE without a cooling system showed serious zero shifts [10,11] which were harmful to the sensors. Heat fluxes were measured for both uncooled and water cooled RDEs using thin film gauge Resistance Temperature Devices [12–14]. The experiments showed that the heat flux at the base of the detonation channel could peak above 9 MW/m2 and below 1 MW/m2 , and heat flux was higher during stable detonation operation than during the unstable operation observed during the ignition time. The heat pulses matched the frequency seen by high-speed video as well as pressure waves in the detonation channel. As a lot of previous numerical simulations employing the Euler equations could not explore the effects of thermal wall conditions on the RDE flow field, viscous simulations of premixed hydrogen-air rotating detonation engines with different wall conditions were carried out [15]. Oblique shock waves were found to exist near the inner and outer walls of the combustor
60
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
annulus, on either side of the RDW front in the simulation with isothermal wall boundary conditions of 800 K. 800 K was roughly the wall temperature for a water cooled annulus made of material commonly used in the lab. The oblique shock structures not found in the previous inviscid study were considered to be caused by the influence of the hot walls on the reactants into which the RDW moved. An RDE with a long duration of use has to be integrated with a cooling system to protect the sensors and combustor. Cooling directly affects the wall temperature which influences the RDW behavior and propulsion performance. Two-dimensional numerical simulations of the RDW are conducted in the present study to explore the effects of different thermal wall conditions on the RDW flow field and to help understand the multiple peaks in one RDW cycle of pressure traces in the RDE run. Why an RDW in an RDE especially without a cooling system can keep rotating continuously and why the RDW is not interrupted by deflagration induced by hot walls will also be explained in the present study.
Fig. 1. Pressure contours with an ignition region of 1.5 MPa and 2800 K at 0 ms, inner diameter 10 mm, outer diameter 20 mm.
2. Computational setup Computational fluid dynamics (CFD) program Fluent is used to simulate RDWs. 2.1. Reaction model The laminar finite-rate model with one step reaction is employed here. The laminar finite-rate model computes the chemical source terms by Arrhenius expressions and neglects the effects of turbulent fluctuations on reactions. The model is exact for laminar flames, but is generally inaccurate for turbulent flames due to highly non-linear Arrhenius chemical kinetics. However, the laminar model may be acceptable for combustion with relatively slow chemistry and small turbulence-chemistry interaction, such as supersonic flames [16]. The forward rate constant for reaction r kf,r is computed using the Arrhenius expression
k f,r = Ar T βr exp(−Er /RT ), where Ar is the pre-exponential factor (consistent units), β r is the temperature exponent (dimensionless), Er is the activation energy for the reaction, and R is the universal gas constant. For the one-step reaction 2H2 + O2 = 2H2 O, Ar is 9.87 × 108 , β r is 0, Er is 1.255 × 108 J/kmol, R is 8314 J/(kmol•K), and r is 1. Rate exponents of H2, O2 and H2 O are respectively 1, 1and 0. Backward reaction is not included. The model of one-step reaction based on Arrhenius expressions was proven to be feasible for RDW study in the inviscid simulations [17]. Since a time step size 0.02 μs was used for computations of RDWs [15], here the time step size is 0.02 μs to keep the computations convergent. A courant number 1 is used, and is reduced automatically for computations near corners. A maximum chemical time step ratio 0.9 is employed to limit the local Courant–Friedrichs–Lewy number when the chemical time scales (eigenvalues of the chemical Jacobian) become too large to maintain a well-conditioned matrix. Inlet diffusion, full multicomponent diffusion and thermal diffusion are used for the species model and more details are found in the reference [16]. The thermal conductivity and viscosity of the species are obtained from mass-weighted mixing law. The mass diffusivity and thermal diffusion coefficient are obtained from kinetic-theory. The constant-pressure specific heat of the material is obtained according to the mixing law. The density of ideal gas is used. 2.2. Other setup Here a two-dimensional model much smaller than the experimental combustor [6,8,10,11] is used to study rotating detonation,
due to the computation speed. The combustor model on x–y plane at ignition time 0 ms is shown in Fig. 1. The inner diameter is 10 mm and the outer diameter is 20 mm. The center of the model is (0, 0). The computational domain is an annular channel between the inner and outer walls. The domain is divided by quadrilateral mesh cells. There are premixed reactants of stoichiometric hydrogen and air mixture in the combustor at the ignition time. Initial absolute pressure and temperature are respectively 0.1 MPa and 300 K. The ignition region shown in Fig. 1 is the intersection of a quad and the computational domain. The xy coordinate ranges of the quad are respectively (3, 10) and (-1, 0) mm. The pressure and temperature of the ignition region are respectively 1.5 MPa and 2800 K which approximate to the Chapman–Jouquet (CJ) detonation parameters of hydrogen-air mixture in ambient conditions. There is no slip on the walls. The adiabatic isolated wall from the point (5, 0) mm to point (10, 0) mm is used to make the detonation wave rotate clockwise. 7 cases are studied by changing thermal conditions of the inner and outer walls. The wall temperatures of 6 cases of them are respectively 40 0 K, 60 0 K, 80 0 K, 10 0 0 K, 120 0 K and 150 0 K. Another case of adiabatic walls is also studied. The inviscid model is not used here because of the heat transfer between the walls and fluids. The standard k-epsilon model has the advantages of robustness, economy, and reasonable accuracy. The RNG k-epsilon model has a higher accuracy for swirling flows and the vorticity computation is better, but it is not economic. Other viscous models such as Realizable k-epsilon, Reynolds Stress and SST k-omega models were used to compute rotating detonation waves by the authors. It was found that for these models variations of the RDW velocity or structure were very small [18]. Here, the results and discussions do not focus on vorticity. The standard k-epsilon model is enough to capture the basic RDW structure and is used in the present study. Reynolds-averaged Navier–Stokes (RANS) equations are solved. The standard wall functions are used for near-wall treatment. Wall Ystar y∗ is a nondimensional parameter defined by the equation
y∗ =
ρCμ1/4 k1P/2 yP , μ
where kP is the turbulence kinetic energy at point P, yp is the distance from point P to the wall, ρ is the fluid density, Cμ is the constant used to compute the turbulent (or eddy) viscosity μt and μ is the fluid viscosity at point P. yp is smaller than the sizes of the wall-adjacent cells. The turbulent (or eddy) viscosity is computed
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
3
by combining k and epsilon as follows:
ε
,
where Cμ is a constant. The logarithmic law for mean velocity is known to be valid for 30 < y∗ < 300. In ANSYS FLUENT, the log-law is employed when y∗ > 11.225. When y∗ is below 11.225 at the wall-adjacent cells, ANSYS FLUENT applies the laminar stress-strain relationship that can be written as U∗ = y∗ , where U∗ is the dimensionless velocity and y∗ is the dimensionless distance from the wall [16]. The lawsof-the-wall for mean velocity and temperature are based on the wall unit y∗ rather than y+. These quantities are approximately equal in equilibrium turbulent boundary layers. The temperature wall functions include the contribution from the viscous heating [19]. The transient density-based solver with implicit formulation and AUSM (Advection Upstream Splitting Method) flux type [20] is used. The density-based solver is appropriate for high speed compressible flow with combustion, hypersonic flows and shock interactions. Generalized Mach number based convection and pressure splitting functions were proposed by Liou [21] and the new scheme was termed AUSM+. The AUSM+ scheme free of oscillations at stationary and moving shocks provides exact resolution of contact and shock discontinuities. Second Order Upwind is used for the pressure, momentum, and energy equations. The advantage of the fully implicit scheme is that it is unconditionally stable with respect to time step size, as indicates that the larger time step size can be used to save the computing time. Second Order Implicit is used for Transient Formulation. Least Squares Cell-Based Gradient Evaluation is used for convection and diffusion terms in the flow conservation equations. On irregular unstructured meshes, the accuracy of the least-squares gradient method is comparable to that of the node-based gradient and they are both better than the cellbased gradient. It is more economical to calculate the least-squares gradient than the node-based gradient [16]. 3. Results and discussions 3.1. Detonation tube Detonation tube cases are employed to verify implementation of the inviscid combustion model in Ansys Fluent before the main research is carried out. The computational setup is the same as that mentioned above except that an inviscid model is employed here. Simulations are performed at an initial pressure and temperature of 10 0, 0 0 0 Pa and 30 0 K, respectively. A 2 mm driver section at 30 0 0 K and 5 MPa is used at the left hand end of the tube to initiate a detonation wave. First, one-dimensional (1D) simulations are conducted using 0.1 mm and 0.2 mm resolution grids to discretize a 0.3 m long tube. Analysis of the simulation resullts enables an evaluation of the numerical solver and chemical reaction model employed here. For this 1D case, computational predictions should be close to the ideal Chapman–Jouguet (CJ) values. The computational peak pressure, peak temperature in the detonation region and detonation velocity for grid size 0.1 mm are respectively 2.81 MPa, 3004 K and 1991 m/s, which are in good agreement with the corresponding results from Gaseq (2.78 MPa, 2956 K and 1979 m/s, respectively). Results from the 0.2 mm resolution simulations are in close agreement with those presented above, with a velocity of 1995 m/s for the detonation wave. Pressure profiles for a 1D premixed stoichiometric hydrogen-air detonation at three time instants are shown in Fig. 2. The waveform indicates that Fluent solver has little numerical dissipation to capture shock waves. Two-dimensional (2D) detonation tube simulations are conducted to demonstrate the ability of the Fluent solver to capture
2.5 Pressure, MPa
μt = ρCμ
k2
61
0.0504 ms 0.06 ms
0.036 ms
2 1.5 1 0.5 0.05
0.1
0.15 x, m
0.2
0.25
Fig. 2. Pressure profiles for a 1D premixed stoichiometric hydrogen-air detonation at three time instants, grid size 0.2 mm.
the cellular structures of a detonation wave. Here, the tube has a height of 0.03 m and a length of 0.3 m. A driver section is employed at the left hand end of the tube, and a quad spot of 1 MPa are placed downstream to produce the instabilities of the detonation wave. The x coordinate ranges of the quad is (2, 5) mm. The y coordinate range is (10, 25) mm. Fig. 3 presents instantaneous pressure and temperature contours at 0.1422 ms before the wave has reached the end of the tube. A maximum pressure time history is shown in Fig. 4, which records the maximum pressure seen at each point in the domain during the detonation propagation. After a short period of time an irregular detonation wave cellular structure is obtained in the simulation. The detonation cell size measurement method is shown in Fig. 4b. The typical cell size marked by the arrows is 8 mm, matching the cell size 8.2 mm in the previous study [22]. However, the cell sizes in Fig. 4 are not uniform and they vary from 6 mm to 10 mm. Actually, if the detonation tube is long enough, detonation cells will tend to be uniform. 3.2. Numerical accuracy for rotating detonation A mesh-convergence study shown in Table 1 and Fig. 5 is conducted for the adiabatic case and the case 1200 K to compare the computational results caused by three sets of quadrilateral mesh cells. The RDW velocity and temperature in Table 1 are on the outer wall at 24 μs for the adiabatic case. The maximum value of y∗ along the walls is situated closely behind the RDW front and on the outer wall, and it is shown in Table 1 for the adiabatic case at 24 μs. The maximum value 290 of y∗ proves the cell size 0.05 mm is feasible for numerical simulations since the logarithmic law is employed when 11.225 < y∗ < 300 and laminar stressstrain relationship is employed when y∗ < 11.225 in ANSYS FLUENT. The heat flux peaks of walls situated closely behind the RDW front and on the outer wall for the case 1200 K at 9.2 μs are also presented in Table 1. Heat flux peaks over 100 MW/m2 were obtained in the experiments [23], which may be referred to. The little difference of heat flux peaks between the cell sizes 0.02 mm and 0.05 mm also proves the cell size 0.05 mm is roughly feasible for the computations. By considering the variation for cell sizes 0.02 mm and 0.05 mm, the theoretical velocity variation is estimated as Dv = Dx/t = 0.03/0.4∗ 10 0 0 m/s = 75 m/s, where Dx is the difference 0.03 mm of mesh cell sizes and t is the time interval 0.4 μs for calculating the RDW velocity. If the numerical velocity variation is smaller than the theoretical velocity variation, it is considered that the computed results are roughly convergent. The RDW velocity difference between the cell sizes 0.05 and 0.02 mm is 28 m/s, which is lower than the theoretical velocity variation 75 m/s. Thus, the numerical results by using the cell size 0.05 mm can be used to analyze the flow field qualitatively as the results are approximate. The numerical rotating detonation velocity 2378 m/s is higher than the velocity 1240 m/s–1630 m/s [6] of the experi-
62
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Pressure, Pa
Temperature, K Fig. 3. Pressure and temperature contours at 0.1422 ms for the 2D detonation tube simulation.
(a) Global view
(b) Close-up view Fig. 4. Maximum pressure history for the 2D detonation tube simulation.
mental rotating detonation for the stoichiometric hydrogen and air mixture as the numerical rotating detonation is two-dimensional while the experimental rotating detonation is three-dimensional. The axial flow perpendicular to the rotating detonation moving direction in the experiment consumes part of the combustion heat from the rotating detonation, causing a lower rotating detonation velocity. Insufficient mixing of fuel and oxidant in the experiment reduces the combustion heat, also causing a lower rotating detonation velocity. The rotating detonation velocity 2378 m/s is also higher than the velocity 2166 m/s for the three-dimensional rotating detonation [17]. The combustion heat loss used to accelerate the axial flow for the three-dimensional RDE causes that. In ad-
dition, the velocity of Chapman–Jouguet detonation for the stoichiometric hydrogen-air mixture in the conditions of 0.1 MPa and 300 K is 1979 m/s, which is lower than the rotating detonation velocity 2378 m/s. The rotating detonation part near the outer wall is compressed and intensified by the outer wall. Thus, it usually has a higher velocity than the rotating detonation part near the inner wall and the Chapman–Jouguet detonation velocity. Pressure and temperature contours at 24 μs for the cell sizes 0.05 mm and 0.02 mm are shown in Fig. 5. The flow field behind the rotating detonation (pointed out by the arrow) is finer for the cell size 0.02 mm and more transverse waves (pointed out by arrows) are obtained. However, the rotating detonation fronts at
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
63
Table 1 Mesh convergence study. Cell size, mm
Cell amount
RDW velocity, m/s
RDW temperature, K
y∗
Heat flux peak, MW/m2
0.02 0.05 0.1
785,500 125,800 31,400
2350 2378 2475
4336 4351 4283
190 290 410
80 74 35
Fig. 5. Pressure (Pa, left) and temperature (K, right) contours at 24 μs, cell size 0.05 mm above, cell size 0.02 mm below, the adiabatic case.
almost the same position are similar to each other. The main results in the study are obtained by using the cell size 0.05 mm, as is limited by the computational conditions. 3.3. Y-shaped RDW A Y-shaped RDW is mainly comprised of an RDW front (RDWF) and a transverse wave (TW). The RDWF usually includes a normal detonation wave (ND) or an oblique detonation wave (OD) and a curved detonation wave (CD), as are divided by the TW. Actually, the Y-shaped RDW is formed from the interaction between the RDWF and TW, and the TW is caused by unsteady or uneven heat release, or pressure gradients in the flow field. Collisions between the RDW and walls result in uneven heat release and large pressure gradients, which are the main causes of the TW in this study. The development of Y-shaped RDW for the wall temperature 1200 K is shown in Fig. 6. At the beginning when the explosion wave develops, it collides with the outer wall, leading to a
reflected shock wave (i.e., TW) at 2.8 μs. While the TW is moving to the inner wall, an obvious Y-shaped RDW shown at 6 μs is formed gradually. Because the outer wall has effects of compression and collision on the RDW, the detonation part near the outer wall is usually intense enough to be normal detonation. The detonation part near the inner wall is weak curved detonation, caused by the expansion effects near the inner wall. A reflected shock wave (RS) occurs at 10 μs after a part of the TW collides with the inner wall. Because the TW moves inward and chases after the RDWF, the region between the RDWF and TW becomes narrower. That region disappears when the TW catches up with the RDWF at 15.6 μs, which causes the Y-shaped RDW to become an RDWF without the TW and intensifies the detonation part near the inner wall. Meanwhile an RS from the inner wall collides with the outer wall to produce a TW and a new Y-shaped RDW will be formed. Thus, the Y-shaped RDW appears after a shock wave collides with the outer wall and disappears when the TW catches up with the RDWF, which occurs periodically. Boundary oblique deto-
64
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Fig. 6. Development of Y-shaped RDW for the case 1200 K, pressure contours.
nation (BOD) caused by high wall temperature has been formed at 18.8 μs and the ND becomes oblique detonation at 24 μs. Close-up views of pressure contours (Pa) and contours of heat of reaction (W) in Fig. 7 show the Y-shaped RDW clearly, which was also discovered in the three-dimensional study of RDEs [24]. The Y-shaped RDW is mainly comprised of an RDWF and a TW at 6 μs. The RDWF is formed from normal detonation and curved detonation. Affected by the walls of high temperature, the RDWF becomes more complicated at 18.8 μs and is formed from two boundary oblique detonations, a normal detonation and a curved detonation. The normal detonation is intensified by the compression effects from the TW and the outer wall and it is more intense than the curved detonation and BODs, as is shown in Fig. 7. Therefore, one conclusion is obtained that a normal detonation is more intense than a curved detonation or BOD in the same reactant conditions. The BODs caused by high wall temperature were also discovered in three-dimensional RDE simulations of isothermal wall conditions [15], as shown in Fig. 8b. It demonstrates that the twodimensional simulation in this study is a valuable reference for the RDE. Contours of heat of reaction with two legends of different maximum levels at 18.8 μs are shown to figure out the transverse wave property and verify BOD. Dashed and dashdot lines respectively for the transverse wave and curved detonation or BOD are drawn where heat of reaction is so little that it cannot be displayed. It is seen that the transverse wave part and the BOD near the normal detonation have large heat of reactions and the BOD near the inner wall has a small heat of reaction, as indicates the TW part near the normal detonation and the BODs are detonation.
The BOD near the inner wall is weaker than the one near the outer wall because of the expansion and compression effects respectively near the inner and outer walls. However, the transverse wave part away from the normal detonation has no heat of reaction, indicating it is a shock wave. It is concluded that the detonation part of the transverse wave leads to the shock part, as is similar to the relation between the rotating detonation and oblique shock in an RDE. It is known in the same way that the transverse wave at 6 μs is also comprised of a detonation wave and a shock wave. The conclusion also applies to the cases of 40 0, 60 0, 80 0, 10 0 0 and 1500 K. However, the BODs for the cases below 800 K are so short that it is difficult to distinguish them from the RDWFs because the low temperatures of the walls have small effects on the kinetic rate of reaction. It is noted that with the BOD development the normal detonation near the outer wall turns gradually to be an oblique detonation wave (OD) shown at 24 μs in Fig. 6. To figure out that what has happened before 24 μs should be known. The inner peak in Fig. 9b is caused by the transverse wave and the other peak is caused by the interaction between the normal detonation and outer BOD. The transverse wave moving to the inner wall intensifies the normal detonation and makes the pressure variation along the normal detonation front except the two sides not large shown in Fig. 9b, causing the detonation velocity variation along the normal detonation to be small. Thus the angular velocity of the normal detonation part near the inner wall is higher and the ND becomes an OD at 24 μs with the inner side stretching forward into the reactants.
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Fig. 7. Close-up view of Y-shaped RDW for the case 1200 K and three-dimensional RDE [24], pressure contours (Pa), contours of heat of reaction (W).
65
66
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Fig. 8. Temperature contours with the wall temperature 800 K, (a) two-dimensional simulation in this study, (b) axial slice through the detonation zone for the threedimensional simulation [15].
5E+06 4.5E+06
Inner side
Outer side
4E+06
Pressure, Pa
3.5E+06 3E+06
2.5E+06 2E+06
1.5E+06 1E+06 500000 0
(a) Pressure contours.
0.005
x,m
0.0055
0.006
(b) Pressure distributions along the x component of ND. Fig. 9. Pressure distributions at 23.2 μs for the case 1200 K.
3.4. Effects of thermal wall conditions on rotating detonation Hot walls heating the reactants near the walls increase the reaction rates and make the detonation near the walls propagate more quickly than other parts, inducing the BOD mentioned above. As is seen in Fig. 10, obvious BODs pointed out by the arrows are formed for cases above 800 K, and the BOD length near the outer wall in Fig. 11a is increased with the increased wall temperature and time. Here the BOD length is a tangential length of the BOD front. The radial thickness h of the non-detonated boundary region heated by hot walls shown in Figs. 11b and 12 is increased with the time, promoting the BOD development. The h is not increased
in some time intervals but the average temperature in the heated region in Fig. 12 is being increased all the time. In fact chemical reactions in the boundary region heated by walls are in progress from 0 ms, which raises the average temperature in the boundary region. However, the chemical reactions are isobaric deflagration and they are much slower than that of detonation. For example, the RDW at 24 μs has rotates around the outer wall for nearly one cycle 56.8 mm but deflagration caused by the outer wall spreads for only 0.25 mm radially for the case 1200 K. Thus it is roughly estimated that detonation is 227 times faster than deflagration for the case 1200 K, which is the reason why an RDW rotating continuously in an RDE without a cooling system is not interrupted
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Fig. 10. Temperature (K) contours at 24 μs for different wall temperatures.
Fig. 11. BOD length near the outer wall and radial thickness h of the non-detonated region heated by the outer wall for the cases 10 0 0 K, 120 0 K and 1500 K.
67
68
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Fig. 12. Development of radial thickness h of the heated region by the outer wall for the case 1200 K, temperature contours (K).
Table 2 The average RDW velocities on the outer wall in different thermal wall conditions. Temperature, K Velocity, m/s
400 2022
600 2043
800 2034
by the deflagration induced by high wall temperature. Higher wall temperature causes a higher RDW velocity and a thicker deflagration region near the walls, which counteract with each other as for detonation proportion in combustion. That makes fuel consumption by detonation not decreased too much, which is helpful for the propulsion performance. Because the walls do not heat the reactants for the adiabatic case shown in Fig. 5, there are no BODs occurring. Pressure distributions along the outer wall at 24 μs are shown in Fig. 13. The angle with the vertex (0, 0) starts from the positive x axis and the positive rotation is counter-clockwise. For example, the point (0, 10) mm is corresponding to the angle π /2. The angle at which the pressure begins to rise in the close-up view (such as that pointed out by the arrow for the case 1500 K) is smaller when the wall temperature is higher, indicating that the RDW is closer to the point (10, 0) mm and the average velocity (from 0 to 24 μs) of the RDW is higher. That is also verified by Table 2. The angles of the pressure rising for the adiabatic case and the case 800 K are both nearly 0.704 rad, as is in accord with the nearly equal velocities in Table 2. The average velocity 2030 m/s of the RDW is lower than the velocity 2378 m/s at 24 μs for the adia-
10 0 0 2059
1200 2143
1500 2258
Adiabatic 2030
batic case because the developing RDW at the beginning is slower than the roughly steady one. There are main peaks and sub-peaks for the RDW of the cases 1500 K, 10 0 0 K and adiabatic shown in Fig. 13, which are also discovered in many experiments such as that in the references [25–27]. Obviously the main peak is caused by the boundary oblique detonation (a part of RDW). A transverse wave TW2 in Fig. 14 is running after the RDW along the outer wall and it causes the sub-peak. The TW2 is caused by the reflected shock wave from the inner wall, which stems from the unsteady heat release due to the collision between the RDW and the outer wall. If one transverse wave dies out due to the viscosity or pressure equilibrium, other TWs will still occur randomly due to the collision between the RDW and the outer wall. Thus, there are always TWs keeping moving back and forth between the inner and outer walls. When the region near the outer wall and behind the boundary oblique detonation is affected by a transverse wave, the sub-peak occurs. That also means if the transverse wave is near the inner wall, there will be no sub-peaks recorded by the sensor in the outer wall, as is seen in Fig. 13 for the cases 400 K, 600 K, 800 K and Fig. 14a. In this two-dimensional study the TWs move radially and tangentially while in an RDE combustor the transverse
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
69
Fig. 13. Pressure distributions along the outer wall at 24 μs for all the cases.
waves propagating along the axial flow will also result in the subpeaks, as is concluded. In the experiments the TWs moving axially have no direct relation to the RDW and therefore the main peak and sub-peak have no time sequence. Either of the two peaks may occur earlier, as shown in Fig. 14a. Because the RDW is more intense than the transverse waves, the main peak is always higher than the sub-peak. It is noted that there may be a few sub-peaks such as those pointed out by arrows in Fig. 14a for one RDW cycle in an RDE because of the multiple transverse waves resulted from the complicated flow field in the experiments. 4. Conclusions A rotating detonation wave in a two-dimensional annular channel with different thermal wall conditions is numerically studied by the standard k-epsilon model. It is found that boundary oblique detonation near the inner and outer walls occurs after some time from ignition for isothermal boundary conditions of high wall temperature above 800 K. The length of the boundary oblique detonation increases with the increasing wall temperature and time because higher wall temperature and more time cause wider heated boundary regions with higher temperature, which increase the kinetic rate of reaction when detonation occurs. Higher wall temperature results in higher average RDW velocity, which is
related to the boundary oblique detonation. As the RDW is much faster (227 times for the case 1200 K) than deflagration due to the high wall temperature, the RDW in an RDE can rotate continuously and most of the reactants will not be combusted by boundary deflagration. The Y-shaped RDW comprised of an RDW front and a transverse wave including detonation and shock is also discovered. The RDW front includes normal detonation and curved detonation. After the RDW develops for a period of time, boundary oblique detonation occurs and is included into the RDW front. The Y-shaped structure appears after the shock wave collides with the outer wall to produce a transverse wave and disappears when the transverse wave catches up with the RDW front. Transverse waves behind or before the RDW produce sub-peaks of pressure traces in the RDE experiments.
Acknowledgments This work was supported by Doctor Research Foundation of Southwest University of Science and Technology [14zx7141], National Natural Science Foundation of China [11502247] and [11602207].
70
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Fig. 14. Generation mechanism of Sub-peak of the RDW pressure traces, (b-c) for the case 1500 K.
Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.compfluid.2016.09. 008. References [1] Falempin F, Naour BL. In: R&T effort on pulsed and continuous detonation wave engines; 2009. p. 7284. [2] Bykovskii FA, Vedernikov EF. Continuous detonation of a subsonic flow of a propellant. Combust. Explo. Shock 2003;39:323–34. [3] Bykovskii FA, Zhdan SA, Vedernikov EF. Continuous spin detonation in annular combustors. Combust. Explo. Shock 2005;41:449–59. [4] Wolanski P. Detonative propulsion. Proc. Combust. Inst. 2013;34:125–58. [5] Lin W, Zhou J, Liu S, Lin Z. An experimental study on CH4/O2 continuously rotating detonation wave in a hollow combustion chamber. Exp. Therm. Fluid Sci. 2015;62:122–30. [6] Rankin BA, Richardson DR, Caswell AW, Naples A, Hoke JL, Schauer FR. Imaging of OH∗ chemiluminescence in an optically accessible nonpremixed rotating detonation engine. In: AIAA; 2015. p. 1604. [7] Nakayama H, Moriya T, Kasahara J, Matsuo A, Sasamoto Y, Funaki I. Stable detonation wave propagation in rectangular-cross-section curved channels. Combust. Flame 2012;159(2):859–69. [8] Fotia M.L., Schauery F., Kaemming T., Hoke J. Study of the experimental performance of a rotating detonation engine with nozzled exhaust flow. AIAA 2015– 0631.
[9] Bykovskii FA, Vedernikov EF. Heat fluxes to combustor walls during continuous spin detonation of fuel-air mixtures. Combust. Explo. Shock 2009;45:70–7. [10] Russo RM, King PI, Schauer FR, Thomas LM. Characterization of pressure rise across a continuous detonation engine. In: AIAA; 2011. p. 6046. [11] Wang YH, Wang JP, Li YS, Li Y. Induction for multiple rotating detonation waves in the hydrogen-oxygen mixture with tangential flow. Int. J. Hydrogen Energ. 2014;39:11792–7. [12] Theuerkauf S, King P, Schauer F, et al. Thermal management for a modular rotating detonation engine. In: AIAA; 2013. p. 1176. [13] Theuerkauf SW, Schauer FR, Anthony R, Hoke JL. Average and instantaneous heat release to the walls of an RDE. In: AIAA; 2014. p. 1503. [14] Theuerkauf SW, Schauer FR, Anthony R, Hoke JL. Experimental characterization of high-frequency heat flux in a rotating detonation engine. In: AIAA; 2015. p. 1603. [15] Cocks PAT, Holley AT, Greene CB, Haas M. Development of a high fidelity RDE simulation capability. In: AIAA; 2015. p. 1823. [16] User’s Manual, ANSYS Fluent 15.0. [17] Tang XM, Wang JP, Shao YT. Three-dimensional numerical investigations of the rotating detonation engine with a hollow combustor. Combust. Flame 2015;162:997–1008. [18] Wang YH. Rotating detonation in a combustor of trapezoidal cross section for the hydrogen-air mixture. Int. J. Hydrogen Energy 2016;41:5605–16. [19] Viegas JR, Rubesin MW, Horstman CC. On the use of wall functions as boundary conditions for two-dimensional separated compressible flows. Technical Report AIAA-85-0180, AIAA 23rd Aerospace Sciences Meeting Reno, Nevada; 1985. [20] Liou MS, Steffen CJ. A new flux splitting scheme. J. Comput. Phys. 1993;107(1):23–39.
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71 [21] Liou MS. A sequel to AUSM: AUSM+. J. Comput. Phys. 1996;129:364–82. [22] Ciccarelli G, Ginsberg T, Boccio J, Finfrock C, Gerlach L, Tagawa H, Malliakos A. Detonation cell size measurements in high-temperature hydrogen-air-stream mixtures at the BNL high-temperature combustion facility. Brookhaven National Laboratory; 1997. Technical Report NUREG/CR-6391, BNL-NUREG-52482. [23] Theuerkauf SW. Heat exchanger design and testing for a 6-inch rotating detonation engine. Graduate School of Engineering and Management, Air Force Institute of Technology (AU). Wright-Patterson AFB OH; 2013. MS Thesis, AFIT-ENY-13-M-33.
71
[24] Zhou R, Wang JP. Numerical investigation of shock wave reflections near the head ends of rotating detonation engines. Shock Waves 2013;23(5):461–72. [25] Wang YH, Wang JP. Effect of equivalence ratio on the velocity of rotating detonation. Int. J. Hydrogen Energy 2015;40:7949–55. [26] George AS, Randall S, Anand V, Driscoll R, Gutmark E. Characterization of initiator dynamics in a rotating detonation combustor. Exp. Therm. Fluid. Sci. 2016;72:171–81. [27] Peng L, Wang D, Wu XS, Ma H, Yang CL. Ignition experiment with automotive spark on rotating detonation engine. Int. J. Hydrogen Energy 2015;40(26):8465–74.
Contents lists available at ScienceDirect
Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid
Effects of thermal wall conditions on rotating detonation Yuhui Wang a,∗, Jianping Wang b, Wenyou Qiao a a b
Research Center of Combustion Aerodynamics, Southwest University of Science and Technology, Mianyang 621010, China College of Engineering, Peking University, Beijing 100871, China
a r t i c l e
i n f o
Article history: Received 19 March 2016 Revised 25 June 2016 Accepted 8 September 2016 Available online 8 September 2016 Keywords: Rotating detonation Boundary oblique detonation Numerical simulation Thermal wall conditions
a b s t r a c t Rotating detonation in an annulus full of stoichiometric hydrogen and air is numerically studied to figure out the effects of thermal wall conditions on detonation. A transient density-based solver with implicit formulation, a laminar finite-rate model with one step reaction, and a standard k-epsilon model are employed. It is found that after some time from ignition, boundary oblique detonation near the walls occurs and is included into the Y-shaped rotating detonation formed from the detonation front and transverse wave. The transverse wave is comprised of detonation attached to the rotating detonation front and shock. Boundary oblique detonation occurs with high wall temperature conditions, but it does not occur with adiabatic wall conditions. It is the rapid reaction caused by the high wall temperature that induces boundary oblique detonation. The length of boundary oblique detonation is increased with increasing wall temperature and time. The average velocity of rotating detonation for higher wall temperature is higher. Why continuous propagation of rotating detonation in engines without cooling is not broken by deflagration due to hot walls is explained. Numerical results show that transverse waves produce sub-peaks of pressure traces in the experiments. © 2016 Elsevier Ltd. All rights reserved.
1. Introduction Detonation has lower entropy production than deflagration, which improves fuel efficiency. Rotating detonation engines (RDEs), which have 20% higher thermal efficiency than deflagration engines [1], have been studied theoretically, numerically and experimentally. An RDE combustor is usually annular and detonation in it is non-premixed combustion in order to avoid backfire. Though the high-pressure region behind the rotating detonation wave (RDW) blocks a part of the inlet, the reactants can still flow into the combustor through the other parts of the inlet to maintain continuous detonation. There may be one or more rotating detonation waves rotating around the combustor channel. Since the detonation wave is almost hypersonic, when the inflow has a lower velocity than the detonation wave it always stays near the combustor inlet where the mixing of fuel and oxidant is insufficient. When the inflow is hypersonic and faster than the detonation wave, a convergent inlet can be used to stabilize the rotating detonation wave axially. Thus, the RDEs apply to air-breathing engines especially. More detailed operating principle can be found in the references [2–5]. OH∗ chemiluminescence imaging was used to capture the RDW structures of different mass flow rates, equivalence ratios, air injection areas and so on [6]. The radial RDW front and the tra-
∗
Corresponding author. Fax: +008608162419256. E-mail address: [email protected] (Y. Wang).
http://dx.doi.org/10.1016/j.compfluid.2016.09.008 0045-7930/© 2016 Elsevier Ltd. All rights reserved.
jectories of triple points in the two-dimensional curved channels were observed using multi-frame short-time open-shutter photography [7]. In the experiments the specific impulse for the RDE of a hydrogen and air mixture with a choked converging area aerospike nozzle was about 430 0–560 0 s for different equivalence ratios of the mass flow rate 1.14 kg/s [8]. Since detonation has higher temperature and pressure than deflagration and the RDW is moving all the time, the cooling of RDEs is a great challenge. The mean heat fluxes to the combustor walls were almost the same for detonation and deflagration, but the peak heat flux during an unsteady operation could be two to three times higher [9]. The pressure traces of the RDW in an RDE without a cooling system showed serious zero shifts [10,11] which were harmful to the sensors. Heat fluxes were measured for both uncooled and water cooled RDEs using thin film gauge Resistance Temperature Devices [12–14]. The experiments showed that the heat flux at the base of the detonation channel could peak above 9 MW/m2 and below 1 MW/m2 , and heat flux was higher during stable detonation operation than during the unstable operation observed during the ignition time. The heat pulses matched the frequency seen by high-speed video as well as pressure waves in the detonation channel. As a lot of previous numerical simulations employing the Euler equations could not explore the effects of thermal wall conditions on the RDE flow field, viscous simulations of premixed hydrogen-air rotating detonation engines with different wall conditions were carried out [15]. Oblique shock waves were found to exist near the inner and outer walls of the combustor
60
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
annulus, on either side of the RDW front in the simulation with isothermal wall boundary conditions of 800 K. 800 K was roughly the wall temperature for a water cooled annulus made of material commonly used in the lab. The oblique shock structures not found in the previous inviscid study were considered to be caused by the influence of the hot walls on the reactants into which the RDW moved. An RDE with a long duration of use has to be integrated with a cooling system to protect the sensors and combustor. Cooling directly affects the wall temperature which influences the RDW behavior and propulsion performance. Two-dimensional numerical simulations of the RDW are conducted in the present study to explore the effects of different thermal wall conditions on the RDW flow field and to help understand the multiple peaks in one RDW cycle of pressure traces in the RDE run. Why an RDW in an RDE especially without a cooling system can keep rotating continuously and why the RDW is not interrupted by deflagration induced by hot walls will also be explained in the present study.
Fig. 1. Pressure contours with an ignition region of 1.5 MPa and 2800 K at 0 ms, inner diameter 10 mm, outer diameter 20 mm.
2. Computational setup Computational fluid dynamics (CFD) program Fluent is used to simulate RDWs. 2.1. Reaction model The laminar finite-rate model with one step reaction is employed here. The laminar finite-rate model computes the chemical source terms by Arrhenius expressions and neglects the effects of turbulent fluctuations on reactions. The model is exact for laminar flames, but is generally inaccurate for turbulent flames due to highly non-linear Arrhenius chemical kinetics. However, the laminar model may be acceptable for combustion with relatively slow chemistry and small turbulence-chemistry interaction, such as supersonic flames [16]. The forward rate constant for reaction r kf,r is computed using the Arrhenius expression
k f,r = Ar T βr exp(−Er /RT ), where Ar is the pre-exponential factor (consistent units), β r is the temperature exponent (dimensionless), Er is the activation energy for the reaction, and R is the universal gas constant. For the one-step reaction 2H2 + O2 = 2H2 O, Ar is 9.87 × 108 , β r is 0, Er is 1.255 × 108 J/kmol, R is 8314 J/(kmol•K), and r is 1. Rate exponents of H2, O2 and H2 O are respectively 1, 1and 0. Backward reaction is not included. The model of one-step reaction based on Arrhenius expressions was proven to be feasible for RDW study in the inviscid simulations [17]. Since a time step size 0.02 μs was used for computations of RDWs [15], here the time step size is 0.02 μs to keep the computations convergent. A courant number 1 is used, and is reduced automatically for computations near corners. A maximum chemical time step ratio 0.9 is employed to limit the local Courant–Friedrichs–Lewy number when the chemical time scales (eigenvalues of the chemical Jacobian) become too large to maintain a well-conditioned matrix. Inlet diffusion, full multicomponent diffusion and thermal diffusion are used for the species model and more details are found in the reference [16]. The thermal conductivity and viscosity of the species are obtained from mass-weighted mixing law. The mass diffusivity and thermal diffusion coefficient are obtained from kinetic-theory. The constant-pressure specific heat of the material is obtained according to the mixing law. The density of ideal gas is used. 2.2. Other setup Here a two-dimensional model much smaller than the experimental combustor [6,8,10,11] is used to study rotating detonation,
due to the computation speed. The combustor model on x–y plane at ignition time 0 ms is shown in Fig. 1. The inner diameter is 10 mm and the outer diameter is 20 mm. The center of the model is (0, 0). The computational domain is an annular channel between the inner and outer walls. The domain is divided by quadrilateral mesh cells. There are premixed reactants of stoichiometric hydrogen and air mixture in the combustor at the ignition time. Initial absolute pressure and temperature are respectively 0.1 MPa and 300 K. The ignition region shown in Fig. 1 is the intersection of a quad and the computational domain. The xy coordinate ranges of the quad are respectively (3, 10) and (-1, 0) mm. The pressure and temperature of the ignition region are respectively 1.5 MPa and 2800 K which approximate to the Chapman–Jouquet (CJ) detonation parameters of hydrogen-air mixture in ambient conditions. There is no slip on the walls. The adiabatic isolated wall from the point (5, 0) mm to point (10, 0) mm is used to make the detonation wave rotate clockwise. 7 cases are studied by changing thermal conditions of the inner and outer walls. The wall temperatures of 6 cases of them are respectively 40 0 K, 60 0 K, 80 0 K, 10 0 0 K, 120 0 K and 150 0 K. Another case of adiabatic walls is also studied. The inviscid model is not used here because of the heat transfer between the walls and fluids. The standard k-epsilon model has the advantages of robustness, economy, and reasonable accuracy. The RNG k-epsilon model has a higher accuracy for swirling flows and the vorticity computation is better, but it is not economic. Other viscous models such as Realizable k-epsilon, Reynolds Stress and SST k-omega models were used to compute rotating detonation waves by the authors. It was found that for these models variations of the RDW velocity or structure were very small [18]. Here, the results and discussions do not focus on vorticity. The standard k-epsilon model is enough to capture the basic RDW structure and is used in the present study. Reynolds-averaged Navier–Stokes (RANS) equations are solved. The standard wall functions are used for near-wall treatment. Wall Ystar y∗ is a nondimensional parameter defined by the equation
y∗ =
ρCμ1/4 k1P/2 yP , μ
where kP is the turbulence kinetic energy at point P, yp is the distance from point P to the wall, ρ is the fluid density, Cμ is the constant used to compute the turbulent (or eddy) viscosity μt and μ is the fluid viscosity at point P. yp is smaller than the sizes of the wall-adjacent cells. The turbulent (or eddy) viscosity is computed
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
3
by combining k and epsilon as follows:
ε
,
where Cμ is a constant. The logarithmic law for mean velocity is known to be valid for 30 < y∗ < 300. In ANSYS FLUENT, the log-law is employed when y∗ > 11.225. When y∗ is below 11.225 at the wall-adjacent cells, ANSYS FLUENT applies the laminar stress-strain relationship that can be written as U∗ = y∗ , where U∗ is the dimensionless velocity and y∗ is the dimensionless distance from the wall [16]. The lawsof-the-wall for mean velocity and temperature are based on the wall unit y∗ rather than y+. These quantities are approximately equal in equilibrium turbulent boundary layers. The temperature wall functions include the contribution from the viscous heating [19]. The transient density-based solver with implicit formulation and AUSM (Advection Upstream Splitting Method) flux type [20] is used. The density-based solver is appropriate for high speed compressible flow with combustion, hypersonic flows and shock interactions. Generalized Mach number based convection and pressure splitting functions were proposed by Liou [21] and the new scheme was termed AUSM+. The AUSM+ scheme free of oscillations at stationary and moving shocks provides exact resolution of contact and shock discontinuities. Second Order Upwind is used for the pressure, momentum, and energy equations. The advantage of the fully implicit scheme is that it is unconditionally stable with respect to time step size, as indicates that the larger time step size can be used to save the computing time. Second Order Implicit is used for Transient Formulation. Least Squares Cell-Based Gradient Evaluation is used for convection and diffusion terms in the flow conservation equations. On irregular unstructured meshes, the accuracy of the least-squares gradient method is comparable to that of the node-based gradient and they are both better than the cellbased gradient. It is more economical to calculate the least-squares gradient than the node-based gradient [16]. 3. Results and discussions 3.1. Detonation tube Detonation tube cases are employed to verify implementation of the inviscid combustion model in Ansys Fluent before the main research is carried out. The computational setup is the same as that mentioned above except that an inviscid model is employed here. Simulations are performed at an initial pressure and temperature of 10 0, 0 0 0 Pa and 30 0 K, respectively. A 2 mm driver section at 30 0 0 K and 5 MPa is used at the left hand end of the tube to initiate a detonation wave. First, one-dimensional (1D) simulations are conducted using 0.1 mm and 0.2 mm resolution grids to discretize a 0.3 m long tube. Analysis of the simulation resullts enables an evaluation of the numerical solver and chemical reaction model employed here. For this 1D case, computational predictions should be close to the ideal Chapman–Jouguet (CJ) values. The computational peak pressure, peak temperature in the detonation region and detonation velocity for grid size 0.1 mm are respectively 2.81 MPa, 3004 K and 1991 m/s, which are in good agreement with the corresponding results from Gaseq (2.78 MPa, 2956 K and 1979 m/s, respectively). Results from the 0.2 mm resolution simulations are in close agreement with those presented above, with a velocity of 1995 m/s for the detonation wave. Pressure profiles for a 1D premixed stoichiometric hydrogen-air detonation at three time instants are shown in Fig. 2. The waveform indicates that Fluent solver has little numerical dissipation to capture shock waves. Two-dimensional (2D) detonation tube simulations are conducted to demonstrate the ability of the Fluent solver to capture
2.5 Pressure, MPa
μt = ρCμ
k2
61
0.0504 ms 0.06 ms
0.036 ms
2 1.5 1 0.5 0.05
0.1
0.15 x, m
0.2
0.25
Fig. 2. Pressure profiles for a 1D premixed stoichiometric hydrogen-air detonation at three time instants, grid size 0.2 mm.
the cellular structures of a detonation wave. Here, the tube has a height of 0.03 m and a length of 0.3 m. A driver section is employed at the left hand end of the tube, and a quad spot of 1 MPa are placed downstream to produce the instabilities of the detonation wave. The x coordinate ranges of the quad is (2, 5) mm. The y coordinate range is (10, 25) mm. Fig. 3 presents instantaneous pressure and temperature contours at 0.1422 ms before the wave has reached the end of the tube. A maximum pressure time history is shown in Fig. 4, which records the maximum pressure seen at each point in the domain during the detonation propagation. After a short period of time an irregular detonation wave cellular structure is obtained in the simulation. The detonation cell size measurement method is shown in Fig. 4b. The typical cell size marked by the arrows is 8 mm, matching the cell size 8.2 mm in the previous study [22]. However, the cell sizes in Fig. 4 are not uniform and they vary from 6 mm to 10 mm. Actually, if the detonation tube is long enough, detonation cells will tend to be uniform. 3.2. Numerical accuracy for rotating detonation A mesh-convergence study shown in Table 1 and Fig. 5 is conducted for the adiabatic case and the case 1200 K to compare the computational results caused by three sets of quadrilateral mesh cells. The RDW velocity and temperature in Table 1 are on the outer wall at 24 μs for the adiabatic case. The maximum value of y∗ along the walls is situated closely behind the RDW front and on the outer wall, and it is shown in Table 1 for the adiabatic case at 24 μs. The maximum value 290 of y∗ proves the cell size 0.05 mm is feasible for numerical simulations since the logarithmic law is employed when 11.225 < y∗ < 300 and laminar stressstrain relationship is employed when y∗ < 11.225 in ANSYS FLUENT. The heat flux peaks of walls situated closely behind the RDW front and on the outer wall for the case 1200 K at 9.2 μs are also presented in Table 1. Heat flux peaks over 100 MW/m2 were obtained in the experiments [23], which may be referred to. The little difference of heat flux peaks between the cell sizes 0.02 mm and 0.05 mm also proves the cell size 0.05 mm is roughly feasible for the computations. By considering the variation for cell sizes 0.02 mm and 0.05 mm, the theoretical velocity variation is estimated as Dv = Dx/t = 0.03/0.4∗ 10 0 0 m/s = 75 m/s, where Dx is the difference 0.03 mm of mesh cell sizes and t is the time interval 0.4 μs for calculating the RDW velocity. If the numerical velocity variation is smaller than the theoretical velocity variation, it is considered that the computed results are roughly convergent. The RDW velocity difference between the cell sizes 0.05 and 0.02 mm is 28 m/s, which is lower than the theoretical velocity variation 75 m/s. Thus, the numerical results by using the cell size 0.05 mm can be used to analyze the flow field qualitatively as the results are approximate. The numerical rotating detonation velocity 2378 m/s is higher than the velocity 1240 m/s–1630 m/s [6] of the experi-
62
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Pressure, Pa
Temperature, K Fig. 3. Pressure and temperature contours at 0.1422 ms for the 2D detonation tube simulation.
(a) Global view
(b) Close-up view Fig. 4. Maximum pressure history for the 2D detonation tube simulation.
mental rotating detonation for the stoichiometric hydrogen and air mixture as the numerical rotating detonation is two-dimensional while the experimental rotating detonation is three-dimensional. The axial flow perpendicular to the rotating detonation moving direction in the experiment consumes part of the combustion heat from the rotating detonation, causing a lower rotating detonation velocity. Insufficient mixing of fuel and oxidant in the experiment reduces the combustion heat, also causing a lower rotating detonation velocity. The rotating detonation velocity 2378 m/s is also higher than the velocity 2166 m/s for the three-dimensional rotating detonation [17]. The combustion heat loss used to accelerate the axial flow for the three-dimensional RDE causes that. In ad-
dition, the velocity of Chapman–Jouguet detonation for the stoichiometric hydrogen-air mixture in the conditions of 0.1 MPa and 300 K is 1979 m/s, which is lower than the rotating detonation velocity 2378 m/s. The rotating detonation part near the outer wall is compressed and intensified by the outer wall. Thus, it usually has a higher velocity than the rotating detonation part near the inner wall and the Chapman–Jouguet detonation velocity. Pressure and temperature contours at 24 μs for the cell sizes 0.05 mm and 0.02 mm are shown in Fig. 5. The flow field behind the rotating detonation (pointed out by the arrow) is finer for the cell size 0.02 mm and more transverse waves (pointed out by arrows) are obtained. However, the rotating detonation fronts at
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
63
Table 1 Mesh convergence study. Cell size, mm
Cell amount
RDW velocity, m/s
RDW temperature, K
y∗
Heat flux peak, MW/m2
0.02 0.05 0.1
785,500 125,800 31,400
2350 2378 2475
4336 4351 4283
190 290 410
80 74 35
Fig. 5. Pressure (Pa, left) and temperature (K, right) contours at 24 μs, cell size 0.05 mm above, cell size 0.02 mm below, the adiabatic case.
almost the same position are similar to each other. The main results in the study are obtained by using the cell size 0.05 mm, as is limited by the computational conditions. 3.3. Y-shaped RDW A Y-shaped RDW is mainly comprised of an RDW front (RDWF) and a transverse wave (TW). The RDWF usually includes a normal detonation wave (ND) or an oblique detonation wave (OD) and a curved detonation wave (CD), as are divided by the TW. Actually, the Y-shaped RDW is formed from the interaction between the RDWF and TW, and the TW is caused by unsteady or uneven heat release, or pressure gradients in the flow field. Collisions between the RDW and walls result in uneven heat release and large pressure gradients, which are the main causes of the TW in this study. The development of Y-shaped RDW for the wall temperature 1200 K is shown in Fig. 6. At the beginning when the explosion wave develops, it collides with the outer wall, leading to a
reflected shock wave (i.e., TW) at 2.8 μs. While the TW is moving to the inner wall, an obvious Y-shaped RDW shown at 6 μs is formed gradually. Because the outer wall has effects of compression and collision on the RDW, the detonation part near the outer wall is usually intense enough to be normal detonation. The detonation part near the inner wall is weak curved detonation, caused by the expansion effects near the inner wall. A reflected shock wave (RS) occurs at 10 μs after a part of the TW collides with the inner wall. Because the TW moves inward and chases after the RDWF, the region between the RDWF and TW becomes narrower. That region disappears when the TW catches up with the RDWF at 15.6 μs, which causes the Y-shaped RDW to become an RDWF without the TW and intensifies the detonation part near the inner wall. Meanwhile an RS from the inner wall collides with the outer wall to produce a TW and a new Y-shaped RDW will be formed. Thus, the Y-shaped RDW appears after a shock wave collides with the outer wall and disappears when the TW catches up with the RDWF, which occurs periodically. Boundary oblique deto-
64
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Fig. 6. Development of Y-shaped RDW for the case 1200 K, pressure contours.
nation (BOD) caused by high wall temperature has been formed at 18.8 μs and the ND becomes oblique detonation at 24 μs. Close-up views of pressure contours (Pa) and contours of heat of reaction (W) in Fig. 7 show the Y-shaped RDW clearly, which was also discovered in the three-dimensional study of RDEs [24]. The Y-shaped RDW is mainly comprised of an RDWF and a TW at 6 μs. The RDWF is formed from normal detonation and curved detonation. Affected by the walls of high temperature, the RDWF becomes more complicated at 18.8 μs and is formed from two boundary oblique detonations, a normal detonation and a curved detonation. The normal detonation is intensified by the compression effects from the TW and the outer wall and it is more intense than the curved detonation and BODs, as is shown in Fig. 7. Therefore, one conclusion is obtained that a normal detonation is more intense than a curved detonation or BOD in the same reactant conditions. The BODs caused by high wall temperature were also discovered in three-dimensional RDE simulations of isothermal wall conditions [15], as shown in Fig. 8b. It demonstrates that the twodimensional simulation in this study is a valuable reference for the RDE. Contours of heat of reaction with two legends of different maximum levels at 18.8 μs are shown to figure out the transverse wave property and verify BOD. Dashed and dashdot lines respectively for the transverse wave and curved detonation or BOD are drawn where heat of reaction is so little that it cannot be displayed. It is seen that the transverse wave part and the BOD near the normal detonation have large heat of reactions and the BOD near the inner wall has a small heat of reaction, as indicates the TW part near the normal detonation and the BODs are detonation.
The BOD near the inner wall is weaker than the one near the outer wall because of the expansion and compression effects respectively near the inner and outer walls. However, the transverse wave part away from the normal detonation has no heat of reaction, indicating it is a shock wave. It is concluded that the detonation part of the transverse wave leads to the shock part, as is similar to the relation between the rotating detonation and oblique shock in an RDE. It is known in the same way that the transverse wave at 6 μs is also comprised of a detonation wave and a shock wave. The conclusion also applies to the cases of 40 0, 60 0, 80 0, 10 0 0 and 1500 K. However, the BODs for the cases below 800 K are so short that it is difficult to distinguish them from the RDWFs because the low temperatures of the walls have small effects on the kinetic rate of reaction. It is noted that with the BOD development the normal detonation near the outer wall turns gradually to be an oblique detonation wave (OD) shown at 24 μs in Fig. 6. To figure out that what has happened before 24 μs should be known. The inner peak in Fig. 9b is caused by the transverse wave and the other peak is caused by the interaction between the normal detonation and outer BOD. The transverse wave moving to the inner wall intensifies the normal detonation and makes the pressure variation along the normal detonation front except the two sides not large shown in Fig. 9b, causing the detonation velocity variation along the normal detonation to be small. Thus the angular velocity of the normal detonation part near the inner wall is higher and the ND becomes an OD at 24 μs with the inner side stretching forward into the reactants.
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Fig. 7. Close-up view of Y-shaped RDW for the case 1200 K and three-dimensional RDE [24], pressure contours (Pa), contours of heat of reaction (W).
65
66
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Fig. 8. Temperature contours with the wall temperature 800 K, (a) two-dimensional simulation in this study, (b) axial slice through the detonation zone for the threedimensional simulation [15].
5E+06 4.5E+06
Inner side
Outer side
4E+06
Pressure, Pa
3.5E+06 3E+06
2.5E+06 2E+06
1.5E+06 1E+06 500000 0
(a) Pressure contours.
0.005
x,m
0.0055
0.006
(b) Pressure distributions along the x component of ND. Fig. 9. Pressure distributions at 23.2 μs for the case 1200 K.
3.4. Effects of thermal wall conditions on rotating detonation Hot walls heating the reactants near the walls increase the reaction rates and make the detonation near the walls propagate more quickly than other parts, inducing the BOD mentioned above. As is seen in Fig. 10, obvious BODs pointed out by the arrows are formed for cases above 800 K, and the BOD length near the outer wall in Fig. 11a is increased with the increased wall temperature and time. Here the BOD length is a tangential length of the BOD front. The radial thickness h of the non-detonated boundary region heated by hot walls shown in Figs. 11b and 12 is increased with the time, promoting the BOD development. The h is not increased
in some time intervals but the average temperature in the heated region in Fig. 12 is being increased all the time. In fact chemical reactions in the boundary region heated by walls are in progress from 0 ms, which raises the average temperature in the boundary region. However, the chemical reactions are isobaric deflagration and they are much slower than that of detonation. For example, the RDW at 24 μs has rotates around the outer wall for nearly one cycle 56.8 mm but deflagration caused by the outer wall spreads for only 0.25 mm radially for the case 1200 K. Thus it is roughly estimated that detonation is 227 times faster than deflagration for the case 1200 K, which is the reason why an RDW rotating continuously in an RDE without a cooling system is not interrupted
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Fig. 10. Temperature (K) contours at 24 μs for different wall temperatures.
Fig. 11. BOD length near the outer wall and radial thickness h of the non-detonated region heated by the outer wall for the cases 10 0 0 K, 120 0 K and 1500 K.
67
68
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Fig. 12. Development of radial thickness h of the heated region by the outer wall for the case 1200 K, temperature contours (K).
Table 2 The average RDW velocities on the outer wall in different thermal wall conditions. Temperature, K Velocity, m/s
400 2022
600 2043
800 2034
by the deflagration induced by high wall temperature. Higher wall temperature causes a higher RDW velocity and a thicker deflagration region near the walls, which counteract with each other as for detonation proportion in combustion. That makes fuel consumption by detonation not decreased too much, which is helpful for the propulsion performance. Because the walls do not heat the reactants for the adiabatic case shown in Fig. 5, there are no BODs occurring. Pressure distributions along the outer wall at 24 μs are shown in Fig. 13. The angle with the vertex (0, 0) starts from the positive x axis and the positive rotation is counter-clockwise. For example, the point (0, 10) mm is corresponding to the angle π /2. The angle at which the pressure begins to rise in the close-up view (such as that pointed out by the arrow for the case 1500 K) is smaller when the wall temperature is higher, indicating that the RDW is closer to the point (10, 0) mm and the average velocity (from 0 to 24 μs) of the RDW is higher. That is also verified by Table 2. The angles of the pressure rising for the adiabatic case and the case 800 K are both nearly 0.704 rad, as is in accord with the nearly equal velocities in Table 2. The average velocity 2030 m/s of the RDW is lower than the velocity 2378 m/s at 24 μs for the adia-
10 0 0 2059
1200 2143
1500 2258
Adiabatic 2030
batic case because the developing RDW at the beginning is slower than the roughly steady one. There are main peaks and sub-peaks for the RDW of the cases 1500 K, 10 0 0 K and adiabatic shown in Fig. 13, which are also discovered in many experiments such as that in the references [25–27]. Obviously the main peak is caused by the boundary oblique detonation (a part of RDW). A transverse wave TW2 in Fig. 14 is running after the RDW along the outer wall and it causes the sub-peak. The TW2 is caused by the reflected shock wave from the inner wall, which stems from the unsteady heat release due to the collision between the RDW and the outer wall. If one transverse wave dies out due to the viscosity or pressure equilibrium, other TWs will still occur randomly due to the collision between the RDW and the outer wall. Thus, there are always TWs keeping moving back and forth between the inner and outer walls. When the region near the outer wall and behind the boundary oblique detonation is affected by a transverse wave, the sub-peak occurs. That also means if the transverse wave is near the inner wall, there will be no sub-peaks recorded by the sensor in the outer wall, as is seen in Fig. 13 for the cases 400 K, 600 K, 800 K and Fig. 14a. In this two-dimensional study the TWs move radially and tangentially while in an RDE combustor the transverse
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
69
Fig. 13. Pressure distributions along the outer wall at 24 μs for all the cases.
waves propagating along the axial flow will also result in the subpeaks, as is concluded. In the experiments the TWs moving axially have no direct relation to the RDW and therefore the main peak and sub-peak have no time sequence. Either of the two peaks may occur earlier, as shown in Fig. 14a. Because the RDW is more intense than the transverse waves, the main peak is always higher than the sub-peak. It is noted that there may be a few sub-peaks such as those pointed out by arrows in Fig. 14a for one RDW cycle in an RDE because of the multiple transverse waves resulted from the complicated flow field in the experiments. 4. Conclusions A rotating detonation wave in a two-dimensional annular channel with different thermal wall conditions is numerically studied by the standard k-epsilon model. It is found that boundary oblique detonation near the inner and outer walls occurs after some time from ignition for isothermal boundary conditions of high wall temperature above 800 K. The length of the boundary oblique detonation increases with the increasing wall temperature and time because higher wall temperature and more time cause wider heated boundary regions with higher temperature, which increase the kinetic rate of reaction when detonation occurs. Higher wall temperature results in higher average RDW velocity, which is
related to the boundary oblique detonation. As the RDW is much faster (227 times for the case 1200 K) than deflagration due to the high wall temperature, the RDW in an RDE can rotate continuously and most of the reactants will not be combusted by boundary deflagration. The Y-shaped RDW comprised of an RDW front and a transverse wave including detonation and shock is also discovered. The RDW front includes normal detonation and curved detonation. After the RDW develops for a period of time, boundary oblique detonation occurs and is included into the RDW front. The Y-shaped structure appears after the shock wave collides with the outer wall to produce a transverse wave and disappears when the transverse wave catches up with the RDW front. Transverse waves behind or before the RDW produce sub-peaks of pressure traces in the RDE experiments.
Acknowledgments This work was supported by Doctor Research Foundation of Southwest University of Science and Technology [14zx7141], National Natural Science Foundation of China [11502247] and [11602207].
70
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71
Fig. 14. Generation mechanism of Sub-peak of the RDW pressure traces, (b-c) for the case 1500 K.
Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.compfluid.2016.09. 008. References [1] Falempin F, Naour BL. In: R&T effort on pulsed and continuous detonation wave engines; 2009. p. 7284. [2] Bykovskii FA, Vedernikov EF. Continuous detonation of a subsonic flow of a propellant. Combust. Explo. Shock 2003;39:323–34. [3] Bykovskii FA, Zhdan SA, Vedernikov EF. Continuous spin detonation in annular combustors. Combust. Explo. Shock 2005;41:449–59. [4] Wolanski P. Detonative propulsion. Proc. Combust. Inst. 2013;34:125–58. [5] Lin W, Zhou J, Liu S, Lin Z. An experimental study on CH4/O2 continuously rotating detonation wave in a hollow combustion chamber. Exp. Therm. Fluid Sci. 2015;62:122–30. [6] Rankin BA, Richardson DR, Caswell AW, Naples A, Hoke JL, Schauer FR. Imaging of OH∗ chemiluminescence in an optically accessible nonpremixed rotating detonation engine. In: AIAA; 2015. p. 1604. [7] Nakayama H, Moriya T, Kasahara J, Matsuo A, Sasamoto Y, Funaki I. Stable detonation wave propagation in rectangular-cross-section curved channels. Combust. Flame 2012;159(2):859–69. [8] Fotia M.L., Schauery F., Kaemming T., Hoke J. Study of the experimental performance of a rotating detonation engine with nozzled exhaust flow. AIAA 2015– 0631.
[9] Bykovskii FA, Vedernikov EF. Heat fluxes to combustor walls during continuous spin detonation of fuel-air mixtures. Combust. Explo. Shock 2009;45:70–7. [10] Russo RM, King PI, Schauer FR, Thomas LM. Characterization of pressure rise across a continuous detonation engine. In: AIAA; 2011. p. 6046. [11] Wang YH, Wang JP, Li YS, Li Y. Induction for multiple rotating detonation waves in the hydrogen-oxygen mixture with tangential flow. Int. J. Hydrogen Energ. 2014;39:11792–7. [12] Theuerkauf S, King P, Schauer F, et al. Thermal management for a modular rotating detonation engine. In: AIAA; 2013. p. 1176. [13] Theuerkauf SW, Schauer FR, Anthony R, Hoke JL. Average and instantaneous heat release to the walls of an RDE. In: AIAA; 2014. p. 1503. [14] Theuerkauf SW, Schauer FR, Anthony R, Hoke JL. Experimental characterization of high-frequency heat flux in a rotating detonation engine. In: AIAA; 2015. p. 1603. [15] Cocks PAT, Holley AT, Greene CB, Haas M. Development of a high fidelity RDE simulation capability. In: AIAA; 2015. p. 1823. [16] User’s Manual, ANSYS Fluent 15.0. [17] Tang XM, Wang JP, Shao YT. Three-dimensional numerical investigations of the rotating detonation engine with a hollow combustor. Combust. Flame 2015;162:997–1008. [18] Wang YH. Rotating detonation in a combustor of trapezoidal cross section for the hydrogen-air mixture. Int. J. Hydrogen Energy 2016;41:5605–16. [19] Viegas JR, Rubesin MW, Horstman CC. On the use of wall functions as boundary conditions for two-dimensional separated compressible flows. Technical Report AIAA-85-0180, AIAA 23rd Aerospace Sciences Meeting Reno, Nevada; 1985. [20] Liou MS, Steffen CJ. A new flux splitting scheme. J. Comput. Phys. 1993;107(1):23–39.
Y. Wang et al. / Computers and Fluids 140 (2016) 59–71 [21] Liou MS. A sequel to AUSM: AUSM+. J. Comput. Phys. 1996;129:364–82. [22] Ciccarelli G, Ginsberg T, Boccio J, Finfrock C, Gerlach L, Tagawa H, Malliakos A. Detonation cell size measurements in high-temperature hydrogen-air-stream mixtures at the BNL high-temperature combustion facility. Brookhaven National Laboratory; 1997. Technical Report NUREG/CR-6391, BNL-NUREG-52482. [23] Theuerkauf SW. Heat exchanger design and testing for a 6-inch rotating detonation engine. Graduate School of Engineering and Management, Air Force Institute of Technology (AU). Wright-Patterson AFB OH; 2013. MS Thesis, AFIT-ENY-13-M-33.
71
[24] Zhou R, Wang JP. Numerical investigation of shock wave reflections near the head ends of rotating detonation engines. Shock Waves 2013;23(5):461–72. [25] Wang YH, Wang JP. Effect of equivalence ratio on the velocity of rotating detonation. Int. J. Hydrogen Energy 2015;40:7949–55. [26] George AS, Randall S, Anand V, Driscoll R, Gutmark E. Characterization of initiator dynamics in a rotating detonation combustor. Exp. Therm. Fluid. Sci. 2016;72:171–81. [27] Peng L, Wang D, Wu XS, Ma H, Yang CL. Ignition experiment with automotive spark on rotating detonation engine. Int. J. Hydrogen Energy 2015;40(26):8465–74.