1. Introduction
An oblique detonation engine (ODE) is an air-breathing hypersonic engine using oblique detonation waves (ODWs) in a combustor (Wolański Reference Wolański2013; Jiang et al. Reference Jiang, Zhang, Liu, Wang and Luo2021; Rosato et al. Reference Rosato, Thornton, Sosa, Bachman, Goodwin and Ahmed2021). When deflagration is replaced by detonation, this propulsion system can achieve a wider flight Mach number and high thermal-cycle efficiency. The initiation structures and mechanisms of ODWs have been widely investigated with gaseous fuels (Teng, Ng & Jiang Reference Teng, Ng and Jiang2017; Martínez-Ruiz et al. Reference Martínez-Ruiz, Huete, Sánchez and Williams2020; Bachman & Goodwin Reference Bachman and Goodwin2021; Teng et al. Reference Teng, Tian, Zhang, Zhou and Ng2021b; Domínguez-González et al. Reference Domínguez-González, Martínez-Ruiz, Scotzniovsky, Sanchez and Williams2022). However, liquid fuel has the advantages of higher energy density and easier storage, and adapting liquid fuels is a critical step towards practical ODE applications.
Many studies have been carried out on the morphology and mechanisms of gaseous ODW initiation. A parametric study is performed to analyse the effect of inflow pressure and Mach number on initiation structure and length (Teng et al. Reference Teng, Ng and Jiang2017). A geometric analysis of two characteristic heights was performed to clarify the mechanism of wave system variation (Teng et al. Reference Teng, Tian, Zhang, Zhou and Ng2021b) and the structure of wedge-induced oblique detonations with small heat release is also investigated (Domínguez-González et al. Reference Domínguez-González, Martínez-Ruiz, Scotzniovsky, Sanchez and Williams2022). The effects of boundary layers are investigated and an ignition criterion to predict the formation of an ODW for a given inflow temperature are established (Bachman & Goodwin Reference Bachman and Goodwin2021). However, these studies used a global or detailed reaction model of gaseous fuels without focusing on ODWs in liquid fuels.
Few studies have been performed on the propagation stability of two-phase ODWs. The initiation and stabilization of two-phase ODWs in kerosene–air mixtures over a wedge are numerically studied (Ren et al. Reference Ren, Wang, Xiang and Zheng2018). The effects of the spray equivalence ratio have also been investigated (Ren et al. Reference Ren, Wang, Xiang and Zheng2018, Reference Ren, Wang, Xiang and Zheng2019). These studies used a two-step reaction mechanism without considering the complex chemical reaction kinetics. A detailed mechanism has been used to find a new ODW initiation structure in acetylene fuel (Zhang et al. Reference Zhang, Fang, Ng and Teng2019). Therefore, it is essential to investigate the propagation stabilities of ODWs in liquid fuels detailed mechanisms. Furthermore, in previous studies (Ren et al. Reference Ren, Wang, Xiang and Zheng2018, Reference Ren, Wang, Xiang and Zheng2019) ODWs have focused on droplets with small diameters (less than 
$10\ \mathrm {\mu }{\rm m}$), and investigations with large droplet diameters are lacking.
In this work, we numerically study the influences of liquid fuel droplet properties on the stability of oblique detonation propagation using the hybrid Eulerian–Lagrangian method. The liquid considered is n-heptane, which is injected into a model ODE through a lean prevaporized n-heptane–air mixture. A parametric analysis is performed to reveal the effects of the droplet diameter of inflow liquid n-heptane on ODW initiation, and unsteady and steady ODWs are observed in the numerical results. This work aims to contribute by: (i) determining how droplet diameter affects the initiation structure and the unsteady mode; (ii) explaining why initiation length and steadiness change with droplet diameter; and (iii) proposing a criterion to identify unsteady and steady ODWs and the validations of this criterion. This paper is organized as follows: the computational method and physical models are introduced in § 2; the results are presented and discussed in § 3; and conclusions are given in § 4.
2. Mathematical and physical models
2.1. Numerical methods
2.1.1. Governing equations
The hybrid Eulerian–Lagrangian method is adopted to study two-phase oblique detonation combustion. For the gas phase, the assumptions are adopted and listed as follows:
(i) the Soret effects and Dufour effects are neglected;
(ii) the gravity force is neglected because the gravity force is smaller than the gas and particle momentum forces,
(iii) thermal radiation is not considered.
Therefore, the Navier–Stokes equations are solved together with the species mass fraction equations and ideal gas equation of state (Huang et al. Reference Huang, Zhao, Xu, Li and Zhang2021; Zhao, Cleary & Zhang Reference Zhao, Cleary and Zhang2021b),
$$\begin{gather} \frac{\partial \rho}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho \boldsymbol{u}) = S_{m}, \end{gather}$$
$$\begin{gather}\frac{\partial(\rho \boldsymbol{u})}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}[\boldsymbol{u}(\rho\boldsymbol{u})] + \boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{T} = \boldsymbol{S}_{\boldsymbol{F}}, \end{gather}$$
$$\begin{gather}\frac{\partial(\rho E)}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}[\boldsymbol{u}(\rho E)] + \boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{u}p) + \boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{T} \boldsymbol{\cdot}\boldsymbol{u}) + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q} = \dot{\omega}_{T}+S_{e}, \end{gather}$$
$$\begin{gather}\frac{\partial (\rho{Y}_{m})}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}[\boldsymbol{u}(\rho Y_m)]+\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{s}_{m} =\dot{\omega}_{m}+S_{Y_m},\quad (m=1,\ldots, M-1), \end{gather}$$
$$\begin{gather}p=\rho R T, \end{gather}$$
where t is time, 
$\rho$ and 
$\boldsymbol {u}$ are the density and velocity vector of gas phase, respectively. Here T is the temperature of gas phase, and p is the pressure of gas phase which can be updated from the equation of state, (2.5), where R is the specific gas constant. Here 
$E$ in (2.3) is the total energy of gas phase. Here 
$\boldsymbol {T}$ in (2.2) and (2.3) is the viscous stress tensor and 
$\boldsymbol {q}$ in (2.3) is the diffusive heat flux. Here 
$\boldsymbol {s}_m$ in (2.4) is the species mass flux. Here 
$Y_{m}$ in (2.4) is the mass fraction of the mth species, where 
$M$ is the total species number. Here 
$\dot {\omega }_T$ in (2.3) is the combustion heat release. Here 
$\dot {\omega }_m$ in (2.4) is the net production rate of the mth species by all N chemical reactions, where 
$N$ is the total reaction number. The definitions of the above parameters are elucidated in the subsequent sections.
The source terms 
$S_m$, 
$\boldsymbol {S}_{\boldsymbol {F}}$, 
$S_e$ and 
$S_{Y_m}$ in (2.1)–(2.4) correspond to interphase exchanges of mass, momentum, energy and species, respectively, thus taking into consideration the full coupling between the continuous gas phase and dispersed liquid phase. Their expressions are given in (2.10)–(2.13). If purely gaseous flows are studied, then these source terms are zero.
2.1.2. Thermodynamics and chemical reaction
The gas total energy 
$E$ is defined as 
$E \equiv e + 0.5 |\boldsymbol {u}|^2$. Here 
$e$ is the specific internal energy, and is defined as 
$e \equiv h_s -p/ \rho$. Here, 
$h_s \equiv \int _{T_0}^{T} C_p \,{\rm d}T+h_{0}$ is the sensible enthalpy, and 
$T_{0}$ and 
$h_{0}$ are reference temperature and sensible enthalpy. Here, 
${C}_p= \sum _{m=1}^{M} Y_{m} {C}_{p,m}$ is the heat capacity at constant pressure, and 
${C}_{p,m}$ is the heat capacity at constant pressure of mth species, estimated from JANAF (joint Army–Navy–Air Force) polynomials (McBride Reference McBride1993). The temperature 
$T$ is calculated from the solved 
$h_s$ by the Newton–Raphson method. Here R in (2.5) is the specific gas constant, calculated from 
${R} = R_u \sum _{m=1}^{M} Y_{m} W_m^{-1}$. Here 
$W_m$ is the molar weight of mth species and 
$R_u = 8.314\ \mathrm {J}\ (\mathrm {mol}\ \mathrm {K})^{-1}$ is the universal gas constant.
In (2.4), 
$\dot {\omega }_m$ is the net production rate of the mth species by all N chemical reactions, and can be calculated from the reaction rate of each elementary reaction, i.e. 
$\dot {\omega }_m = W_m \sum _{j=1}^{N} {\omega }_{m,j}$. Here, 
$N$ is the total reaction number and 
${\omega }_{m,j}$ is calculated from 
${\omega }_{m,j} = (v_{m,j}^{\prime \prime }-v_{m,j}^{\prime })[K_{f,j}\prod _{m=1}^{M} (X_m)^{v_{m,j}^{\prime }}-K_{r,j}\prod _{m=1}^{M} (X_m)^{v_{m,j}^{\prime \prime }}]$. Here 
$v_{m,j}^{\prime \prime }$ and 
$v_{m,j}^{'}$ are the molar stoichiometric coefficients of 
$m$th species in 
$j$th reaction, respectively. Here 
$X_m = \rho Y_m / W_m$ is molar concentration. 
$K_{f,j}$ and 
$K_{r,j}$ are the forward and reverse rates of 
$j$th reaction, respectively. Here 
$K_{f,j}$ can be written in the Arrhenius Law (Law Reference Law2010) as 
$K_{f,j}= A_j T^{\beta _j} \exp [-E_{a,j}/(RT)]$. Here 
$A_j$, 
$\beta _j$ and 
$E_{a,j}$ are, respectively, the pre-exponential factor, temperature exponent and activation energy of the 
$j$th forward reaction. The backwards reaction rate 
$K_{r,j}$ can be evaluated based on 
$K_{f,j}$ and the equilibrium constant of the 
$j$th elementary reaction (Law Reference Law2010). Here 
$\dot {\omega }_T$ in (2.3) is the combustion heat release, and is estimated as 
$\dot {\omega }_T=-\sum _{m=1}^M {\omega }_{m,j}\Delta h_{f,m}$, and 
$\Delta h_{f,m}$ is the formation enthalpy of 
$m$th species. Only (
${M} - 1$) equations are solved in (2.4), and the mass fraction of inert species can be calculated from 
$\sum _{m=1}^{M} Y_{m} =1$.
A skeletal mechanism for 44 species and 112 reactions (Liu et al. Reference Liu, Hewson, Chen and Pitsch2004) was used for n-heptane oxidization. This mechanism has been well validated and agrees well with experimental data such as the ignition delay time over a pressure range of 2.8 to 44 bars (Liu et al. Reference Liu, Hewson, Chen and Pitsch2004) and Chapman–Jouguet velocity (Qi et al. Reference Qi, Dai, Yu and Chen2017). This mechanism has also been applied to study gas detonation propagation and two-phase spray detonation (Qi et al. Reference Qi, Dai, Yu and Chen2017; Zhao, Ren & Zhang Reference Zhao, Ren and Zhang2021c; Zhao & Zhang Reference Zhao and Zhang2021).
2.1.3. Transport coefficients
In (2.3), 
$\boldsymbol {T}$ is the viscous stress tensor, and is modelled by 
$\boldsymbol {T}=-2 \mu \,\mathrm {dev}(\boldsymbol {D})$. Here 
$\mu$ is dynamic viscosity, and is predicted using Sutherland's law, 
$\mu =A_s \sqrt {{T}}/(1+{T}_S/{T})$, which is a temperature-dependent parameter. Following the previous studies (Maragkos, Beji & Merci Reference Maragkos, Beji and Merci2017; Avdonin, Meindl & Polifke Reference Avdonin, Meindl and Polifke2019; Zhao et al. Reference Zhao, Ren and Zhang2021c), independent of species composition, 
$A_s = 1.672 \times 10^{-6}\ {\rm kg}\ ({\rm m} \ {\rm s} \ \sqrt {{\rm K}})^{-1}$ is the Sutherland coefficient, while 
$T_S = 170.672\ {\rm K}$ is the Sutherland temperature. Moreover, 
$\boldsymbol {D}\equiv [\boldsymbol {\nabla } \boldsymbol {u} +(\boldsymbol {\nabla } \boldsymbol {u})^{\mathrm {T}}]/2$ is the deformation gradient tensor and its deviatoric component, i.e. 
${\rm dev}(\boldsymbol {D})$, is defined as 
$\mathrm {dev}(\boldsymbol {D})\equiv \boldsymbol {D}-\mathrm {tr}(\boldsymbol {D})\boldsymbol {I}/3$ with 
$\boldsymbol {I}$ being the unit tensor.
The diffusive heat flux 
$\boldsymbol {q}$ is modelled with Fourier's law, i.e. 
$\boldsymbol {q} = -{k} \boldsymbol {\nabla } T$. Thermal conductivity k is calculated using the Eucken approximation (Poling, Prausnitz & O'Connell Reference Poling, Prausnitz and O'Connell2001), i.e. 
${k} = \mu {C}_v(1.32+1.37{R}/ {C}_v)$. Here 
${C}_v$ is the heat capacity constant volume and derived from 
${C}_v = {C}_p-{R}.$ In (2.4), 
$\boldsymbol {s}_m = -{D}_m\boldsymbol {\nabla }{(\rho Y_{m})}$ is the species mass flux. Here 
${D}_m$ is calculated through 
${D}_m = {k} / (\rho {C}_p$) with a unity Lewis number (
$Le = 1$) assumption.
The validation of above models can be found in a previous study (Huang et al. Reference Huang, Zhao, Xu, Li and Zhang2021). Besides, the Eucken model has also been used in previous detonation studies (Hayashi, Tsuboi & Dzieminska Reference Hayashi, Tsuboi and Dzieminska2020; Ren & Zheng Reference Ren and Zheng2021).
2.1.4. Assumptions of the liquid phase
The basic assumptions of the liquid phase are listed as below.
(i) The liquid phase is modelled as a spray of spherical droplets using the Lagrangian method, with two-way coupling between the gas and liquid phases implemented through the particle-source-in-cell approach (Crowe, Sharma & Stock Reference Crowe, Sharma and Stock1977).
(ii) Because dilute sprays (volume
${\rm fraction} < 0.001$) are considered, point droplet assumption is used. The flows around the droplets are not resolved, and the heat, mass and momentum exchanges between two phases are modelled. The interdroplet interactions are neglected because the momentum response time is much shorter than the particle collision time for dilute sprays (Crowe, Sommerfeld & Tsuji Reference Crowe, Sommerfeld and Tsuji1998).(iii) The temperature inside the droplets is assumed to be uniform because the Biot number is small.
(iv) Since the ratio of gas density to droplet material density is well below one, Basset force, history force and gravity force are neglected (Crowe et al. Reference Crowe, Sharma and Stock1977).
(v) It is assumed that chemical reaction occurs only in the gas phase.
2.1.5. Governing equations for the individual droplets
Therefore, the equations of mass, momentum and energy for the liquid phase are, respectively (Crowe et al. Reference Crowe, Sommerfeld and Tsuji1998),
$$\begin{gather} \frac{{\rm d}m_{d}}{{\rm d}t}=-\dot{m}_{d}, \end{gather}$$
$$\begin{gather}\frac{{\rm d} {\boldsymbol{u}}_{d}}{{\rm d} t}=\frac{{\boldsymbol{F}}_{d} + {\boldsymbol{F}}_p} {m_d}, \end{gather}$$
$$\begin{gather}c_{p, d} \frac{{\rm d} T_d}{{\rm d} t}=\frac{\dot{Q}_c+\dot{Q}_{lat}}{m_d}, \end{gather}$$
where t is time, and 
$m_{d}={\rm \pi} \rho _{d}d^{3}/6$ is the mass of a single droplet, where 
$\rho _{d}$ and d are the droplet density and diameter, respectively. Here 
${\boldsymbol {u}}_{d}$ is the droplet velocity, 
$c_{p, d}$ is the droplet heat capacity and 
$T_{d}$ is the droplet temperature. Infinite thermal conductivity is assumed for each droplet because small droplets are investigated. Here 
$\dot {m}_d$ is the evaporation rate; 
${\boldsymbol {F}}_{d}$ is the Stokes drag force and 
${\boldsymbol {F}}_{p}$ is the pressure gradient force; 
$\dot {Q}_{c}$ is the convective heat transfer between the gas and liquid phases and 
$\dot {Q}_{lat}$ accounts for the heat transfer caused by the latent heat of evaporation. The definitions of the above parameters are elucidated in the subsequent sections.
Moreover, density 
$\rho _{d}$ and heat capacity 
$c_{p, d}$ of the droplet are functions of droplet temperature 
$T_{d}$ (Daubert & Danner Reference Daubert and Danner1985; Perry & Green Reference Perry and Green2007). Specifically, 
$\rho _{d}(T_d)=a_1/[a_2^{1+(1-T_d/a_3)^{a_4}}]$, and 
$c_{p, d}(T_d)=(b_1)^2/\tau +b_2-\tau [2.0 b_1 b_3+\tau (b_1 b_4 +\tau (1/3 \times (b_3)^2 +\tau (0.5 b_3 b_4 +0.2\tau (b_4)^2)))]$, where 
$\tau$ is defined as 
$\tau = 1.0 - \min (T_d,T_c)/T_c$. Here, 
$a_1$, 
$a_2$, 
$a_3$ and 
$a_4$ or 
$b_1$, 
$b_2$, 
$b_3$ and 
$b_4$ are the species-specific constants, and 
$T_c$ is the critical temperature (Daubert & Danner Reference Daubert and Danner1985; Perry & Green Reference Perry and Green2007).
2.1.6. Evaporation model
The modelling of droplet evaporation in this study incorporates several assumptions. The droplet is considered to exhibit spherical symmetry throughout the evaporation process. The evaporation is assumed to occur quasisteadily. The thermal equilibrium at the droplet surface is presumed, with the Lewis number set to unity to signify this equilibrium condition.
The evaporation rate in (2.6) is calculated with the Abramzon and Sirignano model (Abramzon & Sirignano Reference Abramzon and Sirignano1989),
\begin{equation} \dot{m}_d={\rm \pi} d\rho_f D_f \widetilde{S h} \ln \left(1+B_M\right)\!.\end{equation}
Here, 
$\rho _{f} = p_S W_m/(RT_{d,S})$ and 
$D_{f} = 3.6059 \times 10^{-3}(1.8T_{d,S})^{1.75}[\alpha _d/( p_S\beta _d)]$ are the density and mass diffusivity at the film, and are estimated using the one-third rule between the respective gas and liquid quantities (Abramzon & Sirignano Reference Abramzon and Sirignano1989). Here 
$\alpha _d$ and 
$\beta _d$ are the constants for specific species (Fuller, Schettler & Giddings Reference Fuller, Schettler and Giddings1966). The surface vapour pressure 
$p_S$ is estimated from 
$p_S={\rm exp}[c_1+c_2/T_{d,S}+c_3\ln T_{d,S}+c_4{(T_{d,S})}^{c_5}]$. For 
$n$-heptane, the constants, 
$c_1$, 
$c_2$, 
$c_3$, 
$c_4$ and 
$c_5$, can be obtained from a previous study (Perry & Green Reference Perry and Green2007). Moreover, the droplet surface temperature 
$T_{d,S}$ is estimated from 
$T_{d,S} = (T+2T_d)/3$ (Abramzon & Sirignano Reference Abramzon and Sirignano1989). Here 
$T$ is the gas temperature. The modified Sherwood number 
$\widetilde {S h}$ is calculated as 
$\widetilde {Sh}=2+[(1+R e_d S c)^{1 / 3} \max (1, R e_d)^{0.077}-1] / F(B_M)$, with the Schmidt number 
$Sc$ = 1.0. Here 
$Re_{d} \equiv {\rho }d | \boldsymbol {u}_{d}-\boldsymbol {u} |/{\mu }$ is the droplet Reynolds number. Here, 
$\boldsymbol {u}$, 
$\rho$ and 
$\mu$ refer to the velocity vector, density and dynamic viscosity of the gas phase. Here 
$F(B_M)=(1+B_M)^{0.7} \ln (1+B_M) / B_M$ is introduced to consider the variation in film thickness due to Stefan flow (Abramzon & Sirignano Reference Abramzon and Sirignano1989). Here 
$B_{M}\equiv (Y_{Fs}-Y_{F\infty })/(1-Y_{Fs})$ is the Spalding mass transfer number, where 
$Y_{Fs}$ and 
$Y_{F\infty }$ are the fuel vapour mass fractions of the droplet surface and gas phase, respectively. Here 
$Y_{Fs}$ is calculated from 
$Y_{Fs}=W_F X_{Fs}/[W_F X_{Fs}+W_{eF}(1-X_{Fs})]$, where 
$W_F$ is the molecular weight of the vapour, 
$W_{eF}$ is the averaged molecular weight of the mixture excluding the fuel vapour and 
$X_{Fs} = X_F p_{sat} / p$ is the mole fraction of the vapour at the droplet surface. Here, 
$X_F$ is the molar fraction of the condensed species in the gas phase. Here 
$p_{sat}$ is the saturation pressure and calculated based on Raoult's Law (Daubert & Danner Reference Daubert and Danner1985), i.e. 
$p_{sat}=\exp [c_1+c_2/T_{d}+c_3\ln T_{d}+c_4{(T_{d})}^{c_5}]$. Here, the constants, 
$c_1$, 
$c_2$, 
$c_3$, 
$c_4$ and 
$c_5$, are the same value for the surface vapour pressure 
$p_S$.
The Abramzon and Sirignano model (Abramzon & Sirignano Reference Abramzon and Sirignano1989) has also been used in other two-phase detonation simulations (Watanabe, Matsuo & Matsuoka Reference Watanabe, Matsuo and Matsuoka2019).
2.1.7. Drag force and pressure gradient force
Following a previous study (Meng et al. Reference Meng, Zhao, Xu, Zhang and Zhang2023), the effects of compressibility and droplet volume fraction are not considered in the drag force due to dilute small droplets. The Stokes drag is used in (2.7) and modelled as 
$\boldsymbol {F}_{d}=({18 \mu }/{\rho _d d^2})({C_d R e_d}/{24}) m_d{(\boldsymbol {u}-\boldsymbol {u}_{d})}$ (Liu, Mather & Reitz Reference Liu, Mather and Reitz1993), where 
$C_{d}$ is the drag coefficient and estimated using the Schiller and Naumann model (Naumann & Schiller Reference Naumann and Schiller1935). The validation can be found in a previous study (Huang et al. Reference Huang, Zhao, Xu, Li and Zhang2021).
A large pressure gradient can be observed when droplets pass through the shocks. Therefore, the pressure gradient force is calculated as 
$\boldsymbol {F}_{p}=-V_d \boldsymbol {\nabla } p$, where 
$V_d$ is the volume of a single fuel droplet.
2.1.8. Heat transfer between gas and liquid phases
In (2.8), 
$\dot {Q}_{c}=h_{c}A_{d}(T-T_{d})$ denotes the convective heat transfer between the gas and liquid phases, where 
$A_{d}$ is the surface area of a single droplet; 
$h_{c}$ is the convective heat transfer coefficient (Marshall & Ranz Reference Marshall and Ranz1952), and is estimated using the correlation of Ranz and Marshall through the modified Nusselt number, 
$\widetilde {N u}=2+[(1+R e_d P r)^{1 / 3} \max (1, R e_d)^{0.077}-1] / F(B_T)$, where 
$P_{r}$ is the gas Prandtl number (assumed to be unity here) and 
$B_{T} = {(1+B_M)}^{(C_{p,v}/C_{p,d})/Le}$ is the Spalding heat transfer number; 
$Le$ refers to the Lewis number of the gas mixture and 
$C_{p,v}$ is the constant pressure specific heat of vapour. The validation of this convective heat transfer model can be found in a previous study (Huang et al. Reference Huang, Zhao, Xu, Li and Zhang2021).
Furthermore, 
$\dot {Q}_{lat}=-\dot {m}_dh(T_d)$ in (2.8) accounts for the heat transfer caused by the latent heat of evaporation, and 
$h(T_d)$ is the fuel vapour enthalpy at the droplet temperature 
$T_d$ (Perry & Green Reference Perry and Green2007). 
$h(T_d)$ can be estimated from 
$h(T_d) = d_1(1-T_r)^{[(d_2 T_r+d_3)T_r+d_4]T_r+d_5}$, where 
$d_1$, 
$d_2$, 
$d_3$, 
$d_4$ and 
$d_5$ are the constants and can be obtained from a previous study (Perry & Green Reference Perry and Green2007), and 
$T_r$ is defined as 
$T_r=T_d/T_c$.
2.1.9. Breakup model
The breakup condition in the simulation is established based on several assumptions. It is assumed that the mass of the droplet remains constant both before and after breaking up. A uniform temperature is considered for the droplet just before and just after breaking up. The assumption is made that the droplet maintains a spherical shape in both stages of the breakup process. These assumptions help define the conditions for the breakup phenomenon in the simulation.
Large droplets breakup into finer droplets primarily near the shock front (Watanabe et al. Reference Watanabe, Matsuo, Chinnayya, Matsuoka, Kawasaki and Kasahara2021). Finer droplets are observed to evaporate faster and absorb more heat from the gas phase. As in previous two-phase detonation studies (Shi et al. Reference Shi, Xu, Ren and Zhang2022; Xu & Zhang Reference Xu and Zhang2022), we model droplet breakup in terms of the aerodynamic forces resulting from shock impact, following the method proposed by Pilch & Erdman (Reference Pilch and Erdman1987). Typically, droplets breakup when the Weber number (
$We$) exceeds 12 (Craig et al. Reference Craig, Moharreri, Schanot, Rogers, Anderson and Dhaniyala2013), although the Ohnesorge number can also play a role in this process. Various breakup modes can be modelled on the basis of the total breakup time (Pilch & Erdman Reference Pilch and Erdman1987). Besides, this model also shows a good agreement with experiment data when modelling the interaction between a cloud of water droplets and a planar shock wave (Chauvin et al. Reference Chauvin, Daniel, Chinnayya, Massoni and Jourdan2016). This model has been used in successive detonation researches (Watanabe et al. Reference Watanabe, Matsuo, Chinnayya, Matsuoka, Kawasaki and Kasahara2020; Jourdaine, Tsuboi & Hayashi Reference Jourdaine, Tsuboi and Hayashi2022; Shi et al. Reference Shi, Xu, Ren and Zhang2022; Xu & Zhang Reference Xu and Zhang2022), representing reasonable results.
2.1.10. Source terms of the liquid phase
Two-way coupling between the gas and liquid phases is considered in terms of exchanges of mass 
$S_m$, momentum 
${\boldsymbol {S}}_{\boldsymbol {F}}$, energy 
$S_e$ and species 
$S_{Y_m}$. Therefore, the source terms for the foregoing equation in gas phase read
$$\begin{gather} S_m=\frac{1}{V_c} \sum_1^{N_d} \dot{m}_d, \end{gather}$$
$$\begin{gather}{\boldsymbol{S}}_{F}=-\frac{1}{V_c} \sum_1^{N_d} (\boldsymbol{F}_{p} + \boldsymbol{F}_{d}), \end{gather}$$
$$\begin{gather}S_e=-\frac{1}{V_c} \sum_1^{N_d}\left(\dot{Q}_c+\dot{Q}_{l a t}\right)\!, \end{gather}$$
$$\begin{gather}S_{Y_m}= \begin{cases} S_m & \text{for the liquid fuel species, } \\ 0 & \text{for other species,} \end{cases} \end{gather}$$
where 
$V_{c}$ is the computational fluid dynamics cell volume defined as the volume of a discretized computational cell in a numerical simulation, and 
$N_{d}$ is the droplet number in a cell.
2.2. Liquid-fuelled ODE model
Figure 1 shows the schematic of an ODE and a two-dimensional computational domain, which is used here to mimic an ODE combustor. Following previous ODE work (Teng et al. Reference Teng, Tian, Zhang, Zhou and Ng2021b), the engine inlet wave configuration is proposed by Dudebout, Sislian & Oppitz (Reference Dudebout, Sislian and Oppitz1998) and adopted in later research (Sislian et al. Reference Sislian, Schirmer, Dudebout and Schumacher2001). A difference from Sislian et al. (Reference Sislian, Schirmer, Dudebout and Schumacher2001) is that the air inflow at an altitude of 30 km is assumed to be compressed by two equal-strength oblique shock waves (OSWs) with a total deflection angle of 
$20^\circ$ to minimize the entropy increase, which is essential to improve the propulsion performance. Following previous studies (Bian, Zhou & Teng Reference Bian, Zhou and Teng2021; Teng et al. Reference Teng, Tian, Zhang, Zhou and Ng2021b), the injection process is not modelled, and the fuel–air system is assumed to be well premixed. This assumption is made because these complexities of physical phenomenon are influenced by various gas-dynamic and geometric parameters, which are challenging to determine for lack of referential engines so far. An ODE typically operates within a flight Mach number range of 7 to 12. This study assumes that the liquid-fuelled ODE operates at a flight Mach number 
$M_0$ of 9, owing to the challenges of initiating liquid-fuelled ODEs at low flight Mach numbers and the potential mixing issues caused by high-speed air flow in the inlet duct, which may limit their application at high Mach numbers. Under these assumptions, the inflow parameters of the ODE combustor can be deduced from the Rankine–Hugoniot relation, and the inflow temperature 
$T_3$, pressure 
$P_3$, velocity 
$V_3$ and Mach number 
$M_3$ are 697 K, 28 554 Pa, 
$2535\ {\rm m}\ {\rm s}^{-1}$ and 4.8, respectively. The supersonic homogeneous inflow reflects on the two-dimensional wedge to generate an OSW. To address the thermal requirements for droplet evaporation, a relatively large wedge angle 
$\theta$ of 
$27^\circ$ is chosen to generate an OSW with high postshock temperature. The computational domain is displayed and enclosed by a dashed line in figure 1, and the grid used in this study is uniform and structured.
Figure 1. Schematic of oblique detonation engine and computational domain.
For the boundary conditions, inflow conditions are used for the left-hand boundary of the computational zone. The right-hand and upper boundaries are interpolated under the assumption of zero gradient for all flow parameters. The slip reflecting boundary condition is used on the wedge surface. It should be noted that some recent studies (Fang, Zhang & Hu Reference Fang, Zhang and Hu2019; Bachman & Goodwin Reference Bachman and Goodwin2021) have demonstrated that the boundary layer may change the ODW structures, but the effects are weak and limited near the wedge in some cases with thick boundary layers. For the cases in this study, the thickness of the boundary layer near the ignition position is below 5.0 % of the induction zone height, and the Reynolds number (
$Re$) is of the order of 
$10^6$. Early studies (Li, Kailasanath & Oran Reference Li, Kailasanath and Oran1993) pointed out that boundary layer effects are negligible under very high Reynolds number, and most theoretical studies of the ODW (Teng et al. Reference Teng, Ng and Jiang2017, Reference Teng, Tian, Zhang, Zhou and Ng2021b; Zhang et al. Reference Zhang, Fang, Ng and Teng2019) are based on the slip reflecting boundary. Therefore, to inherit previous ODW theories, the wall boundary on the wedge is modelled with the geometrical constraints of the grid system using the slip reflecting boundary condition. Moreover, we have considered the effects of different boundary conditions, and the effect of the boundary layer is presented in Appendix B.
For all cases in this study, mixtures of liquid or gaseous n-heptane and air were used as inflow mixtures. Three case groups parameterized by the liquid phase equivalence ratio 
$\phi _{l}$ were considered, as listed in table 1. The total equivalence ratio was assumed as 
$\phi _{t}=\phi _{g}+\phi _{l}=1$ for all the two-phase cases, whilst the gas equivalence ratio 
$\phi _{g}$ varied from 0.3 to 0.7, corresponding to 
$\phi _{l} = 0.7\unicode{x2013}0.3$ in cases B, C and D. Here, 
$\phi _{l}$ is defined as the mass ratio of the droplets to the oxidizer normalized by the mass ratio of n-heptane vapour to air under stoichiometric conditions. This study assumed that the droplet velocity is equal to the gas velocity, and the droplet temperature was set to 300 K. The droplets were injected uniformly, and complex distribution models of droplets were not considered. Following the previous studies of liquid ODW (Ren et al. Reference Ren, Wang, Xiang and Zheng2018; Guo et al. Reference Guo, Xu, Li and Zhang2022), a rebound condition is used to model the behaviour of droplets upon collision with the wall. This condition entails reversing the normal component of the droplet velocities after the collision. Both the computational domain and mesh scale were adjusted according to the multiscale nature of the phenomenon. The numerical grid resolution used for different cases varied from the coarsest mesh of 
$400\ \mathrm {\mu } {\rm m}$ to the finest mesh of 
$50\ \mathrm {\mu } {\rm m}$. Droplet diameters of 
$2\ \mathrm {\mu } {\rm m}$, 
$10\ \mathrm {\mu } {\rm m}$, 
$20\ \mathrm {\mu } {\rm m}$ and 
$30 \ \mathrm {\mu } {\rm m}$ were simulated with mesh sizes of 
$400\ \mathrm {\mu } {\rm m}$, 
$400 \ \mathrm {\mu } {\rm m}$, 
$400\ \mathrm {\mu } {\rm m}$ and 
$100\ \mathrm {\mu } {\rm m}$, respectively. These give ratios of droplet diameter to grid resolution of 0.005, 0.025, 0.05 and 0.3, respectively. We have conducted several tests with different ratios, and these tests have demonstrated that the simulation results remain insensitive to this ratio when it is smaller than unity, which aligns well with previous study (Sontheimer, Kronenburg & Stein Reference Sontheimer, Kronenburg and Stein2021). The initiation structure was the main concern in this study. There were at least 500 grids in the induction zone, and as we examined, initiation lengths slightly changed with finer mesh resolution. In each group, various values of the initial diameter of a monodispersed droplet, 
$d_0$, were considered, ranging from 2 to 
$40\ {\mathrm {\mu }}{\rm m}$. This study used various droplet diameters to investigate the effect of droplet diameter, and droplet volume changes in a large range from 
$8\ \mathrm {\mu } {\rm m}^3$ to 
$64\,000\ \mathrm {\mu } {\rm m}^{3}$. Following a previous study (Crowe et al. Reference Crowe, Sommerfeld and Tsuji1998), the mesh grid volumes were adjusted to ensure that the adequate number of Lagrangian particles were used in simulation and confirm the dilute characteristics of the spray ODE. The volume fraction 
$\alpha$ of the injected spray was below 0.001 in all cases, as shown in table 1. The number of tracked Lagrangian particles ranged from 55 000 to 550 000 because of the complex droplet evaporation with ODWs and the scales of liquid ODWs. The study of the particle number dependence is presented in Appendix A to ensure the results of this study are not influenced by the particle number used. Using the variables declared earlier, the number of droplets per unit volume can be calculated in a non-specified manner. For example, the number of droplets is 
$5.26 \times 10^{14}$ per cubic metre when 
$\phi _l = 0.5$, 
$\phi _g = 0.5$, and 
$d_0 = 2\ \mathrm {\mu } {\rm m}$ and 
$1.56 \times 10^{11}$ per cubic metre when 
$\phi _l = 0.5$, 
$\phi _g = 0.5$, and 
$d_0 = 30\ \mathrm {\mu } {\rm m}$. Purely gaseous cases with 
$\phi _{g}= 0.5$ and 1.0 were also considered for comparison.
Table 1. Liquid fuel spray information.
2.3. The solver used in this study
In the gas phase, the equations (2.1)–(2.4) are discretized using a cell-centred finite volume method. Temporal discretization is achieved through a second-order implicit backward method and the time step is adjusted to ensure that the maximum Courant number is less than 0.2. The convective terms in the momentum equations utilize the second-order MUSCL-type (‘monotonic upstream-centred scheme for conservation laws’-type) Riemann-solver-free scheme by a previous study (Kurganov, Noelle & Petrova Reference Kurganov, Noelle and Petrova2001), while the total variation diminishing scheme is applied to the convective terms in the energy and species equations. Additionally, the diffusion terms in (2.2)–(2.4) are discretized using a second-order central differencing scheme. The chemistry integration employs an Euler implicit method, and its accuracy has been validated against other ordinary differential equation solvers in previous work (Huang et al. Reference Huang, Zhao, Xu, Li and Zhang2021). For the liquid phase, water droplets are tracked based on their barycentric coordinates, and (2.6)–(2.8) are solved using a first-order Euler method. The right-hand terms, such as 
$\dot {m}_d$ in (2.6), are integrated in a semi-implicit approach. Gas properties at the droplet location are interpolated from gas phase simulation results.
Both the gas and liquid equations were solved using a multicomponent, two-phase and reactive solver, RYrhoCentralFoam, with two-way interphasic coupling in terms of mass, momentum, energy and species through (2.10)–(2.13). The RYrhoCentralFoam solver has been validated for gaseous supersonic flows and detonative combustion problems (Jasak Reference Jasak1996; Zhao et al. Reference Zhao, Li, Teo, Khoo and Zhang2020, Reference Zhao, Ren and Zhang2021c,Reference Zhao, Chen, Zhang and Swaminathana). Satisfactory accuracies are achieved in terms of capturing the supersonic flow discontinuity, detonation propagation speed and detonation cell size in hydrogen–air mixtures (Jasak Reference Jasak1996; Zhao et al. Reference Zhao, Li, Teo, Khoo and Zhang2020). It has also been validated and used to study two-phase detonation problems (Zhao et al. Reference Zhao, Cleary and Zhang2021b, Reference Zhao, Li, Teo, Khoo and Zhang2020; Zhao & Zhang Reference Zhao and Zhang2021; Guo et al. Reference Guo, Xu, Li and Zhang2022). More information about the numerical schemes and solution strategies can be found in the previous studies (Zhao et al. Reference Zhao, Li, Teo, Khoo and Zhang2020; Huang et al. Reference Huang, Zhao, Xu, Li and Zhang2021).
3. Results and discussion
3.1. Oblique detonation waves in a pure gas mixture
Figure 2 shows the distributions of temperature, pressure, density gradient, heat release rate (
$HRR$) and 
${\rm CH}_{2}{\rm O}$ and 
${\rm OH}$ mass fractions for case A with 
$\phi _{g} = 0.5$ and 
$\phi _{l} = 0.0$. It is a typically smooth OSW–ODW transition similar to the ‘type II’ wave system in a previous study (Teng et al. Reference Teng, Tian, Zhang, Zhou and Ng2021b) for the purely gaseous case. The basic ODE structures are captured in figures 2(a)–2(c), including the main ODW, OSW with a deflection angle of 
$10.2^\circ$, smooth OSW–ODW transition, compression waves, reflection waves and slip line. The temperature behind the OSW is approximately 1600 K (figure 2a), where some heat (figure 2e) is subsequently released owing to the induction reaction and deflagration. The distributions of some key species such as 
${\rm CH}_{2}{\rm O}$ and 
${\rm OH}$ are also shown in figures 2(e) and 2(f). 
${\rm CH}_{2}{\rm O}$ is mainly distributed behind the OSW, while the small-molecule intermediate products such as 
${\rm OH}$ radicals (figure 2f) are mainly distributed behind the main ODW.
Figure 2. Distributions of (a) temperature, (b) pressure, (c) density gradient, (d) heat release rate, (e) OH mass fraction and (f) 
${\rm CH}_{2}{\rm O}$ mass fraction for case A with 
$\phi _{g} = 0.5$ and 
$\phi _{l} = 0.0$. The dashed line in (d–f) represents the OSW.
Figure 3 shows the distributions of temperature, pressure, density gradient, heat release rate and 
${\rm CH}_{2}{\rm O}$ and OH mass fractions for case A with 
$\phi _{g} = 1.0$ and 
$\phi _{l} = 0.0$. The basic ODE structures are also captured in figures 3(a)–3(c), including the main ODW, secondary ODW, OSW with a deflection angle of 
$9.1^\circ$, compression wave, reflection waves, triple point and slip line. These structures show that the ODW morphology in this work is similar to the ‘type III’ wave systems in a previous study (Teng et al. Reference Teng, Tian, Zhang, Zhou and Ng2021b), for which the structures are stable. For 
$\phi _{g} = 0.5$ and 
$\phi _{l} = 0.0$, the morphology is more like a type II wave system. It should be noted that a high equivalence ratio is expected to increase the heat release and then generate an ODW surface with a higher degree so that the OSW–ODW transition in case A with 
$\phi _{g} = 1.0$ and 
$\phi _{l} = 0.0$ becomes abrupt with a triple point. The temperature behind the OSW is approximately 1400 K (figure 3a), which is lower than that in the case A with 
$\phi _{g} = 0.5$ and 
$\phi _{l} = 0.0$. Therefore, the location of the OSW–ODW transition in the case A with 
$\phi _{g} = 1.0$ and 
$\phi _{l} = 0.0$, which is approximately around 0.06 m as shown in figure 3(a), moves farther than that in case A with 
$\phi _{g} = 1.0$ and 
$\phi _{l} = 0.0$, which is approximately 0.03 m in figure 3(a).
Figure 3. Distributions of (a) temperature (K), (b) pressure (Pa), (c) density gradient (
${\rm kg}\ {\rm m}^{-4}$), (d) heat release rate (
${\rm J}\ {\rm m}^{-3}\ {\rm s}^{-1}$), (e) OH mass fraction and (f) 
${\rm CH}_{2}{\rm O}$ mass fraction for case A with 
$\phi _{g} = 1.0$ and 
$\phi _{l} = 0.0$. The dashed line in (d–f) represents the OSW.
In addition to the detonative combustion behind the main ODW and secondary ODW, an intense deflagration with a high heat release rate is also observed in figure 3(d) owing to the compression and reflection waves. However, no obvious release of chemical reaction heat occurs behind the OSW. This absence of heat release can be attributed to the fact that the chemical reaction time scale of n-heptane is larger than the flow time scale under the given conditions. Although there is no obvious heat release behind the OSW in figure 3(d), the chemical reaction does occur in this region, which is clearly seen with the distributions of key species concentrations such as the 
${\rm CH}_{2}{\rm O}$ mass fraction in figure 3(e). We also find that the 
${\rm CH}_{2}{\rm O}$ is mainly distributed behind the OSW, while the small-molecule intermediate products such as 
${\rm OH}$ radicals are mainly distributed behind the main and secondary ODWs in figure 3(f).
The goal of this work is to study the influences of n-heptane droplets on oblique detonation with a skeletal chemical mechanism. Therefore, the droplet evaporation, heat transfer and two-phase momentum exchange all affect the detonation structures and oblique detonation propagation when fuel droplets are present. These will be further analysed in the following sections. Moreover, the mesh resolution sensitivity for purely gaseous and partially prevaporized liquid ODWs is analysed in Appendix A.
3.2. Oblique detonation wave structures in n-heptane sprays
As shown in table 1, initial droplet diameters of 
$2\unicode{x2013}40 \ \mathrm {\mu } {\rm m}$ were considered. This section discusses the effects of n-heptane spray on ODW structures for the cases with 
$\phi _{l} = 0.5$ and 
$d_{0} = 2\ \mathrm {\mu } {\rm m}$, 
$20\ \mathrm {\mu } {\rm m}$ and 
$30\ \mathrm {\mu } {\rm m}$, which represent typical working conditions.
Figures 4 and 5 show the gas temperature and pressure distributions, respectively, for the case with 
$\phi _{l} = 0.5$, 
$\phi _{g} = 0.5$ and 
$d_0= 2\ \mathrm {\mu } {\rm m}$. Six different times, 
$t = 1.02$, 1.03, 1.04, 1.05, 1.07 and 1.10 ms, are presented to study the unsteady characteristics of an ODE fuelled with n-heptane sprays. The main ODW structures are similar to ‘type IV’ in a previous study (Teng et al. Reference Teng, Tian, Zhang, Zhou and Ng2021b), for which there is a normal detonation wave (NDW) under the triple point. However, new features arise owing to the presence of dispersed droplets. In figures 4(a) and 5(a), an explosion point appears in front of the NDW in the initiation zone, which evolves into a new NDW. Figures 4(b)–4(d) and 5(b)–5(d) show the formation of a new NDW from the explosion point. The explosion point initiates a new combustion wave, which is initiated from the wedge and moves both upwards to the OSW and backwards to downstream. Finally, the combustion wave evolves into a new NDW followed with the quenching of the previous NDW. Figures 4(e) and 5(e) show that the new NDW moves backwards. When the NDW moves farther backwards in figures 4(f) and 5(f), a new explosion point appears again in front of the NDW in the initiation zone. This new explosion point subsequently generates a new NDW at the succeeding times. Generally, this is an unsteady process in which the NDW always appears within a certain distance, moves backwards, disappears and re-forms in the case with 
$\phi _{l} = 0.5$ and 
$d_0 = 2\ \mathrm {\mu } {\rm m}$.
Figure 4. Temperature distributions at different times for 
$\phi _{l} = 0.5$, 
$\phi _{g} = 0.5$ and 
$d_{0} = 2\ \mathrm {\mu } {\rm m}$.
Figure 5. Pressure distributions at different times for 
$\phi _{l} = 0.5$, 
$\phi _{g} = 0.5$ and 
$d_{0} = 2\ \mathrm {\mu } {\rm m}$.
To further explain the unsteady behaviours in figures 4 and 5, figure 6 shows the distributions of the droplet volume fraction 
$\alpha$, droplet evaporation rate, droplet heat transfer rate, temperature gradient, mole fraction of 
${\rm C}_{7}{\rm H}_{16}$ and heat release rate for the case with 
$\phi _{l} = 0.5$, 
$\phi _{g} = 0.5$ and 
$d_0 = 2\ \mathrm {\mu } {\rm m}$. It is seen from figure 6(a) that the volume fraction of liquid droplets is much less than 0.001, which satisfies the sparse Lagrangian hypothesis (Crowe et al. Reference Crowe, Sommerfeld and Tsuji1998). In figures 6(a)–6(c), the droplets evaporate near the entrance of the combustor, and the droplets undergo both mass and heat transfer before the OSW. Because of the droplet evaporation, the mean temperature in the induction zone decreases to approximately 1320 K, which is lower than that of the purely gaseous case with 
$\phi _{g} = 1.0$, and therefore, the location of the OSW–ODW transition (approximately 0.27 m in figure 4) moves farther backwards than that in case A with 
$\phi _{g} = 1.0$ and 
$\phi _{l} = 0.0$. There may be two mechanisms underlying this unsteady behaviour. First, as shown in figure 6(d), the non-uniform distribution of temperature gradients suggests that there are some cool spots with positive temperature gradients, which would cause a detonation wave according to previous study (Dai et al. Reference Dai, Chen, Chen and Ju2015). Second, the non-uniform radical distribution in the induction zone causes a spatially inhomogeneous heat release reaction, as shown in figure 6(e). Both mechanisms would lead to the explosion point with intense deflagration appearing in the induction zone (figure 6e) in front of the NDW. Notably, the reason for the unsteady structures caused by the spray droplets is different from that in previous work (Yang, Ng & Teng Reference Yang, Ng and Teng2019), in which the periodic inflow condition was used.
Figure 6. Distributions of the (a) droplet volume fraction 
$\alpha$, (b) droplet evaporation rate, (c) droplet heat transfer rate, (d) temperature gradient, (e) the mole fraction of 
${\rm C}_{7}{\rm H}_{16}$, (f) heat release rate, (g) droplet velocity, (h) droplet temperature and (i) droplet diameter at time 
${t} = 1.02\ {\rm ms}$ for 
$\phi _{l} = 0.5$, 
$\phi _{g} = 0.5$ and 
$d_{0} = 2\ \mathrm {\mu } {\rm m}$.
The droplet velocity, temperature and diameter distributions are presented in figures 7(g)–7(i) to show the behaviour of the droplets with small diameters. Before the OSW, most of their velocities remain relatively constant. The droplet temperature increases by approximately 20 K over a short distance. After this quick heating, the droplets evaporate within a certain distance. A few droplets pass through the OSW and experience a significant velocity decrease, and significant temperature increase.
Figure 7. Distributions of (a) temperature at different times, and (b) droplet volume fraction 
$\alpha$, (c) droplet heat transfer rate, (d) heat release rate, (e) droplet velocity, (f) droplet temperature and (g) droplet diameter at time 
${t} = 0.80$ ms for 
$\phi _{l} = 0.5$, 
$\phi _{g} = 0.5$ and 
$d_{0} = 20\ \mathrm {\mu } {\rm m}$.
Figure 7 shows the distributions of the gas temperature, droplet volume fraction 
$\alpha$, droplet heat transfer rate, and heat release rate for a larger droplet diameter, where 
$\phi _{l} = 0.5$, 
$\phi _{g} = 0.5$ and 
$d_{0} = 20\ \mathrm {\mu } {\rm m}$. Generally, the basic structure and unsteady behaviour of the ODW with 
$d_{0} = 20\ \mathrm {\mu } {\rm m}$ is similar to that with 
$d_{0} = 2\ \mathrm {\mu } {\rm m}$. Additionally, a deflagration is observed behind the ODW in figures 7(a) and 7(d). This is due to the distinct location of droplet evaporation. In figure 7(b), the droplets do not evaporate significantly before the OSW or ODW, but evaporate completely within a certain distance behind the OSW and ODW. The colour bar of figure 7(b) shows that the heat transfer rate before the OSW and ODW ranges between 
$1\times 10^6$ and 
$1\times 10^7$, which is lower than that for 
$d_{0} = 2\ \mathrm {\mu } {\rm m}$ in figure 6(c), indicating only a small amount of evaporation before the OSW and ODW. It is noteworthy that the heat transfer rates behind the ODW are greater than those behind the OSW, mainly owing to the higher temperature generated by the combustion. We define the region from the start of the OSW and ODW to the end of droplet evaporation as the ‘droplet evaporation layer’ in figure 7(b). Therefore, the oblique detonation wavefront mainly comes from the reaction of initial pre-evaporated n-heptane, while the deflagration comes from the reaction of evaporated n-heptane droplets in the droplet evaporation layer behind the ODW. Furthermore, the location of the OSW–ODW transition (approximately 0.02 m) is shifted forwards compared with that for 
$d_{0} = 2\ \mathrm {\mu } {\rm m}$. This shift is attributed to the fact that the OSW for 
$d_{0} = 20\ \mathrm {\mu } {\rm m}$ is mainly induced by the initial pre-evaporated n-heptane. As a result, the temperature in the evaporation layer (approximately 1500 K) becomes higher than in the rest of the induction zone, where it is approximately 1300 K (figure 7a). The high temperature leads to higher chemical reaction rates in the evaporation layer in figure 7(d), which makes the OSW–ODW transition move forwards (Teng et al. Reference Teng, Ng and Jiang2017, Reference Teng, Tian, Zhang, Zhou and Ng2021b).
The droplet velocity, temperature and diameter distributions are presented in figures 7(e), 7(f) and 7(g) to show droplet behaviour. The droplet velocity decreases significantly after the droplets pass through the OSW and ODW, especially the ODW, which has a larger pressure gradient force on the droplets. The droplet heating and breakup is presented in figures 7(f) and 7(g). The temperature of droplets with 
$d_{0} = 20\ \mathrm {\mu } {\rm m}$ increases slowly before the OSW and ODW because individual droplets have large mass. Behind the OSW and ODW, the high temperature of the postshock mixtures raises the droplet temperature immediately, and the big droplets break up into small droplets. After the heating and breakup, the droplets evaporate in a certain distance. Compared with the distances of velocity loss, heating and breakup, the distance over which small droplets evaporate following the breakup dominates the scale of the evaporation layer. Generally, droplets with an initial diameter of 
$20\ \mathrm {\mu } {\rm m}$ experience a rapid decrease in velocity, rapid heating and breakup behind the OSW and ODW. After this process, the droplets break up into smaller droplets and then evaporate over a certain distance.
Figure 8 shows the distributions of the temperature, pressure, evaporation rate, droplet heat transfer rate, 
${\rm C}_{7}{\rm H}_{16}$ mole fraction and heat release rate for 
$\phi _{l} = 0.5$ and 
$\phi _{g} = 0.5$ when the initial droplet diameter further increases to 
$d_0 = 30\ \mathrm {\mu } {\rm m}$. The basic ODE structures are captured in figures 8(a)–8(f), including the main ODW with a cellular-like structure on its wavefront, OSW, OSW–ODW transition, slip line, reflection shock wave and deflagration behind the ODW. The OSW and ODW are connected by a curved shock wave, therefore leading to a smooth OSW–ODW transition (Teng et al. Reference Teng, Ng and Jiang2017). This smooth transition is caused by a high temperature of approximately 1584 K (figure 8a) in the droplet evaporation layer in the induction zone, which will also be discussed in § 3.3. Compared with 
$d_0 = 20 \ \mathrm {\mu } {\rm m}$ in figure 7, an initial diameter of 
$30\ \mathrm {\mu } {\rm m}$ shows a relatively steady state. The formation and mechanism of a steady ODW with a large droplet diameter will be further discussed in §§ 3.3 and 3.4. Compared with the pure-gas case with the same gas equivalence ratio in figure 2(b), both cases exhibit a similarly smooth wave structure. However, a notable difference between wave structures is observed for 
$d_0 = 30\ \mathrm {\mu } {\rm m}$, where a significant increase in pressure and temperature is observed behind the transition from the OSW to the main ODW in the region enclosed by the dashed oval in figures 8(a) and 8(b). The pressure and temperature increases in the study are mainly caused by local heat release from the combustion of unburned gas mixture, which forms due to the evaporation of liquid fuels. This heat release leads to a thermal expansion effect that interacts with the supersonic inflow, forming compression waves, resulting in the observed pressure and temperature rise.
Figure 8. Distributions of (a) temperature, (b) pressure, (c) evaporation rate, (d) droplet heat transfer rate, (e) 
${\rm C}_{7}{\rm H}_{16}$ mole fraction, (f) heat release rate, (g) droplet velocity, (h) droplet temperature, (i) droplet diameter for 
$\phi _{l} =0.5$, 
$\phi _{g} = 0.5$ and 
$d_{0} = 30 \ \mathrm {\mu } {\rm m}$.
In figures 8(c) and 8(d), the occurrence of droplet evaporation before the OSW and ODW is minimal, while a larger number of droplets are observed behind the OSW and ODW. This can be attributed to the fact that larger liquid droplets require a longer time to undergo complete evaporation. This leads to a longer finite thickness of the evaporation layer, and a longer distance is observed between the oblique detonation wavefront and the deflagration behind the ODW. In figures 8(d) and 8(f), there is a high heat transfer behind the ODW, and this deflagration maintains the temperature behind the ODW. It can be observed in figure 8(c) that some droplets move ahead of the OSW near the leading edge owing to the upwards bounce of droplets on the wedge. This phenomenon has also been reported and explained by Guo et al. (Reference Guo, Xu, Li and Zhang2022), who suggest that it occurs because the momentum exchange between the droplets and the gas is not fully achieved in the narrow region behind the OSW. Furthermore, a sharp triangular region with lower temperature forms in the initiation zone in figure 8(a). This phenomenon is also reported by this previous study (Guo et al. Reference Guo, Xu, Li and Zhang2022), who suggest that the heating and evaporation of high-concentration n-heptane droplets change the local thermo-chemical state in zone B (figure 8e), thereby bifurcating the reaction. Intense deflagration occurs at the edge of the sharp triangular region, which is clearly seen in figure 8(f), whereas there is almost no heat release from chemical reactions inside this region.
The droplet velocity, temperature and diameter distributions are presented in figures 8(g), 8(h) and 8(i) to show droplet behaviour. The velocity of droplets with 
$d_0 = 30\ \mathrm {\mu } {\rm m}$ decreases more slowly than that of droplets with 
$d_0 = 20\ \mathrm {\mu } {\rm m}$ owing to the large Stokes number behind the OSW and ODW. The droplet temperature increases even more slowly before the OSW and ODW, and the droplets do not significantly break up. Behind the OSW and ODW, droplets with a diameter of 
$30\ \mathrm {\mu } {\rm m}$ require a longer distance to increase their temperature and breakup than those with a diameter of 
$20\ \mathrm {\mu }{\rm m}$. The break up of droplets behind the OSW and ODW significantly reduces the thickness of the evaporation layer. Apparently, the droplets after breakup also require a longer distance to fully evaporate. For 
$d_0 = 30\ \mathrm {\mu } {\rm m}$, the thickness of the evaporation layer generally depends on the distance over which the droplets heat, break up and evaporate. The break up of droplets behind the OSW and ODW significantly reduces the thickness of the evaporation layer.
To further clearly present the evolution of two-phase ODW, the supplementary movies available at https://doi.org/10.1017/jfm.2024.194 present an animation of the cases with 
$\phi _{g} = 0.5$, 
$\phi _{l} = 0.5$, 
$d_0 = 20\ \mathrm {\mu } {\rm m}$ and 
$30\ \mathrm {\mu } {\rm m}$.
3.3. Effect of initial droplet diameters on the initiation length and steadiness of liquid ODWs
The initiation length is one of the most important parameters in an ODW and is mainly affected by the inflow conditions and wedge angle in pure-gas cases (Teng et al. Reference Teng, Tian, Zhang, Zhou and Ng2021b). Following a previous study (Teng et al. Reference Teng, Ng and Jiang2017), the initiation length here is defined along the streamline close to the wedge, from the front tip to the end of the induction zone, where the temperature increases to 110 % of the postshock temperature. Figure 9 shows the numerical initiation lengths for the two-phase ODW for 
$\phi _{l} = 0.3$, 0.5 and 0.7 with different initial droplet diameters. The total equivalence ratio 
$\phi _{l}+\phi _{g}$ is equal to 1. The initiation length is found to increase first (
$d_{0} <10\ \mathrm {\mu } {\rm m}$) and then decrease (
$d_{0} >10\ \mathrm {\mu } {\rm m}$) with the initial droplet diameter for each pre-evaporation equivalence ratio 
$\phi _{g}$. The variation in the fluctuation range of initiation lengths with different droplet diameters follows a trend similar to that of the initiation length, increasing and then decreasing to a constant value. The fluctuation ranges of initiation length decreases to nearly zero when the ODW is steady. When the initial droplet diameter is approximately 
$10\ \mathrm {\mu } {\rm m}$, the ODW has the longest initiation length. However, as the initial droplet diameter further increases from 30 to 
$40\ \mathrm {\mu } {\rm m}$, the initiation length does not change with the initial droplet diameter. Similar unsteady detonation behaviours are also observed in the cases with fluctuating initiation lengths in figure 9, including the explosion point, combustion wave propagation and formation of a new NDW. This leads to a cyclical pattern in the initiation length, where it first slowly increases when the NDW is moving downward, then suddenly decreases when the new NDW forms, then slowly increases again and finally decreases suddenly again. However, for the cases with triangular markers in figure 9, a relatively steady detonation propagation is observed with little change in the initiation length. It is worth noting that, in cases with small droplets in figure 9, the fluctuation in initiation length decreases as the initial droplet diameter decreases. If the initial droplet diameter further decreases to an extremely low level (such as 
$0.1\ \mathrm {\mu } {\rm m}$ or 
$0.01 \ \mathrm {\mu } {\rm m}$), the ODW should be steady. However, this study is mainly concerned about liquid ODWs with the common droplet diameter range from 
$1\ \mathrm {\mu } {\rm m}$ to 
$40 \ \mathrm {\mu } {\rm m}$, and the steadiness of ODWs with extremely small droplets will be investigated in the future.
Figure 9. Numerical initiation lengths for 
$\phi _{l} = 0.3$, 0.5 and 0.7 with different initial droplet diameters.
To investigate how initial droplet diameter affects initiation length, it is necessary to first understand how droplet evaporation changes gas mixture parameters. This is demonstrated in figure 10 through the temperature, pressure and velocity profiles along various streamlines for droplet diameters of 
$2\ \mathrm {\mu } {\rm m}$, 
$10 \ \mathrm {\mu } {\rm m}$, 
$20\ \mathrm {\mu } {\rm m}$ and 
$30\ \mathrm {\mu } {\rm m}$. The different streamlines in figure 10 correspond to different flow regions of the ODW field. The red lines travel through the OSW–ODW transition, the black ones travel along the main ODW surface, and lines with other colours travel through the induction OSW. The location of the evaporation layer in each case is marked in figures 10(ai), 10(bi) and 10(ci).
Figure 10. Temperature, pressure and velocity profiles of the gas phase along the different streamlines for 
$\phi _{l} = 0.5$ and (a) 
$d_{0} = 2\ \mathrm {\mu } {\rm m}$, (b) 
$d_{0} = 10 \ \mathrm {\mu } {\rm m}$, (c) 
$d_{0} = 20\ \mathrm {\mu } {\rm m}$ and (d) 
$d_{0} = 30\ \mathrm {\mu } {\rm m}$.
We first discuss the effect of droplet evaporation on gas temperature for different droplet diameters. For a droplet diameter of 
$2\ \mathrm {\mu } {\rm m}$ in figure 10(ai), preshock temperatures are approximately 500 K for all streamlines because the droplets evaporate quickly before the OSW and ODW. As the droplets completely evaporate before the OSW and ODW, the temperature in the induction zone remains nearly constant before the heat release and the postshock temperature decreases from approximately 1400 K for the pure-gas inflow case in figure 3(a) to approximately 1320 K. As the droplet diameter increases to 
$10\ \mathrm {\mu } {\rm m}$ in figure 10(bi), the preshock temperature of the bottom streamline is higher than for a droplet diameter of 
$2 \ \mathrm {\mu } {\rm m}$ owing to the more challenging evaporation process for larger droplets. As the streamlines progress from the bottom to the top of the computational region, the preshock temperature decreases and then becomes relatively constant. The droplets with a diameter of 
$10\ \mathrm {\mu } {\rm m}$ injected from the upper region travel a longer distance, reaching the OSW or ODW from the upstream to the downstream, leading to more droplets evaporating before the wavefront and a decrease in preshock temperature. The postshock temperature of the streamline along the wedge first jumps to a high peak at approximately 1500 K owing to the high preshock temperature. The temperature in induction zone then decreases to a low level for droplet evaporation in the evaporation layer behind the OSW. For a streamline (
${y} = 0.1$ m and 0.2 m in figure 10bi) travelling through the upper region of the induction zone, the temperature peak and evaporation layer disappear because the droplets completely evaporate before the wavefront, and the temperature in the induction zone remains constant at a low level.
As the droplet diameter increases to 20 and 
$30\ \mathrm {\mu } {\rm m}$, there is no cooling of the preshock temperature because the droplets do not evaporate significantly before the wavefront. As streamlines traverse the induction zone for the case of 
$20\ \mathrm {\mu } {\rm m}$ in figure 10(ci), the temperatures in the induction zone exhibit a similar trend, jumping to a high level at approximately 1500 K, remaining constant over a short distance, and decreasing more slowly to a lower level in the evaporation layer than that in the case of 
$10\ \mathrm {\mu } {\rm m}$. It is worth noting that for 
$20\ \mathrm {\mu } {\rm m}$, the high temperature in the evaporation layer behind the OSW is sustained for a longer distance than that for 
$10\ \mathrm {\mu } {\rm m}$. For 
$d_0 = 30\ \mathrm {\mu } {\rm m}$ in figure 10(di), the temperature behind the OSW jumps to 1500 K and the temperature in the evaporation layer increases after remaining constant over a short distance. The reason for this temperature change for 20 and 
$30\ \mathrm {\mu } {\rm m}$ is discussed in § 3.4, and involves the competition between the droplet heat transfer and chemical reaction heat release. In general, the preshock temperature decreases for small droplet diameters, and the postshock temperature jumps to a high peak if droplets incompletely evaporate before encountering the induction OSW.
The effect of droplet evaporation on gas pressure and velocity for different droplet diameters will now be discussed. Generally, it is observed in figures 10(b) and 10(c) that the preshock velocity and pressure in different streamlines change slightly in all cases, leading to similar postshock pressure and velocity. The pressure and velocity in the induction zone remain constant before the heat release. This demonstrates that droplet evaporation has little effect on gas pressure and velocity.
As discussed above, the preshock and postshock parameters change with the various droplet diameters. Figure 11 shows the deflection angle of the induction OSW for different initial droplet diameters with 
$\phi _{l} = 0.5$. The deflection angle has a large decrease when the inflow fuels mixture changes from gas to liquid owing to the decrease in preshock temperature. The deflection angle of the induction OSW slightly increases as the droplet diameter increases to 
$15\ \mathrm {\mu } {\rm m}$ as a result of an overall increase in preshock temperature. As the droplet diameter further increases from 15 to 
$25\ \mathrm {\mu } {\rm m}$, the deflection angle significantly increases. The high temperature peak in the evaporation layer can sustain a longer distance behind the OSW than that in the case of 
$10\ \mathrm {\mu } {\rm m}$, as shown in figure 10(c), leading to a decrease in density and increase in deflection angle. For droplet diameters from 
$15\ \mathrm {\mu } {\rm m}$ to 
$25 \ \mathrm {\mu } {\rm m}$, the preshock parameters in figure 10(c) are nearly the same as those in case A with 
$\phi _{g} = 0.5$ and 
$\phi _{l} = 0.0$ in figure 2, but the deflection shock is smaller than that in case A. This is due to the increase in density resulting from the temperature drop in the induction zone after the temperature peak caused by droplet evaporation, which leads to a smaller deflection angle. As the droplet diameter is further increased to 
$30\ \mathrm {\mu } {\rm m}$, there is no temperature drop in the induction zone, and the deflection angle remains constant at approximately 
$10.2^\circ$, which is similar to the deflection angle in case A. As discussed above, droplet evaporation has little effect on gas pressure and velocity, and it has been suggested that the initiation temperature has an important effect on the initiation length (Teng et al. Reference Teng, Tian, Zhang, Zhou and Ng2021b). Figure 12(a) shows the average temperature in the initiation zone for 
$\phi _{l}= 0.3$, 0.5 and 0.7 with different initial droplet diameters to explain the variation in the initiation lengths. For initial droplet diameters less than 
$10 \ \mathrm {\mu } {\rm m}$, the amount of droplet evaporation causes the temperature in the initiation zone to decrease with the initial droplet diameter. When the initial droplet diameter is greater than 
$10\ \mathrm {\mu } {\rm m}$, the temperature in the initiation zone increases with the initial droplet diameter. Figure 12(b) shows the temperature profiles along the streamline near the wall for initial droplet diameters of 
$d_{0} = 10$, 20 and 
$30 \ \mathrm {\mu } {\rm m}$ for comparison. For 
$d_{0} = 20\ \mathrm {\mu } {\rm m}$, fewer droplets evaporate before reaching the OSW and ODW than for 
$d_{0}= 10\ \mathrm {\mu } {\rm m}$. This results in a thicker droplet evaporation layer, which corresponds to a region of higher temperature. As a consequence, there is a greater heat release due to the higher heat release rate of the induction reaction of the initially pre-evaporated n-heptane. This leads to higher average temperatures in the initiation zone, resulting in shorter initiation lengths and fluctuations in their ranges. Therefore, the droplet evaporation in the initiation zone has important effects on the oblique detonation propagation in an ODE.
Figure 11. Deflection angle of induction OSW for 
$\phi _{l} = 0.5$ with different initial droplet diameters.
Figure 12. (a) Average temperatures in the initiation zone for 
$\phi _{l} = 0.3$, 0.5 and 0.7 with different initial droplet diameters and (b) temperature profiles along the streamline near the wall for 
$\phi _{l} = 0.5$, 
$d_{0} = 10$, 20 and 
$30\ \mathrm {\mu } {\rm m}$.
3.4. Formation of the steady ODW with large droplet diameters
Figure 13 shows the profiles of the droplet heat transfer rate and heat release rate along the streamline near the wall for 
$\phi _{l} = 0.5$ and 
$d_{0} = 20$ and 
$30 \ \mathrm {\mu } {\rm m}$ to clarify the formation of steady ODWs with large droplet diameters. For 
$d_{0} = 20\ \mathrm {\mu } {\rm m}$, most droplet evaporation occurs in front of the location of maximum reaction heat release, as shown in figure 13(a). However, this is contrary to the case with 
$d_{0} = 30\ \mathrm {\mu } {\rm m}$. The droplet evaporation has little effect on the initiation zone because most droplet evaporation occurs downstream. This is why the temperature does not significantly decrease in the initiation zone for 
$d_{0} = 30\ \mathrm {\mu } {\rm m}$ in figure 13(b). Generally, unsteady ODW behaviour occurs when dispersed droplets evaporate before the induction reaction finishes. When droplets evaporate after the induction reaction finishes, the ODW is steady and the deflagration behind the OSW can be observed (figure 8f).
Figure 13. Profiles of droplet heat transfer rate (
${\rm J}\ {\rm m}^{-3}\ {\rm s}^{-1}$) and heat release rate (
${\rm J}\ {\rm m}^{-3}\ {\rm s}^{-1}$) along the streamline near the wall for 
$\phi _{l} = 0.5$ and (a) 
$d_{0} = 20\ \mathrm {\mu } {\rm m}$ and (b) 
$d_{0} = 30 \ \mathrm {\mu } {\rm m}$.
The observed changes in temperatures in figure 10(ci) and 10(di) can be attributed to the interplay between droplet heat transfer and heat release rate. In figure 13(a), significant heat transfer occurs at some distance behind the OSW, sustaining the high temperature behind the OSW for a longer distance than that for 
$10\ \mathrm {\mu } {\rm m}$. Furthermore, it is evident that the heat transfer rate peaks at a longer distance for 
$d_{0} = 30\ \mathrm {\mu } {\rm m}$ than for 
$d_{0} = 20\ \mathrm {\mu } {\rm m}$. This observation suggests that larger droplets require more time to evaporate and break up. This can be attributed to the fact that droplets with larger initial diameters have a smaller total surface area, resulting in reduced interaction between the phases. As the droplet heat transfer for 
$d_{0} = 20\ \mathrm {\mu } {\rm m}$ occurs before the induction reaction finishes, significant heat transfer can cause a decrease in the induction temperature, resulting in a delay in induction reaction time. For 
$d_{0} = 30\ \mathrm {\mu } {\rm m}$, after the induction reaction, the liquid fuel droplets that have not yet evaporated continue to evaporate and react with postshock gaseous mixtures, leading to the second heat release peak. This maintains the increase in the postshock temperature, as shown in figure 10(di). Additionally, this increase in temperature causes an increase in droplet heat transfer rate, as shown in figure 13. The first heat transfer rate peak for 
$d_{0} = 20 \ \mathrm {\mu } {\rm m}$ is approximately 
$1.5\times 10^{10}\ {\rm J}\ {\rm m}^{-3}\ {\rm s}^{-1}$, which is smaller than that for 
$d_{0} = 30\ \mathrm {\mu } {\rm m}$.
As discussed above, the propagation instabilities of the ODW with partially prevaporized sprays is related to the parameters in the induction zone. The variation in non-dimensional numbers with droplet diameter in the induction zone are discussed. Notably, we utilized data spanning over 1000 instants to calculate the average parameters relevant for the non-dimensional numbers in this study. The behaviour of droplets suspended in the gas flow in the induction zone behind the OSW is assessed through the Stokes number. This Stokes number 
$St$ is calculated as 
$St={\bar {t}_{a,i}}/{\bar {t}_i},$ where 
$\bar {t}_{i}$ is the characteristic main flow time of the induction zone, which is the time-averaged ODW induction time of the streamline near wedge, and 
$\bar {t}_{a,i}$ is the average droplet acceleration time, which is defined as 
${{\bar {t}_{a,i}} = { \bar {\rho }_{d,i} ({\bar {d}_{i}})^2} / ({18 {\bar {\mu }_{i}}}),}$ where 
$\bar {\rho }_{d,i}$ and 
${\bar {d}_{i}}$ are the average droplet density and average droplet diameter in the induction zone, and 
${\bar {\mu }_{i}}$ are the average viscosity of gas mixture calculated through the average temperature in the induction zone. As shown in figure 14(a), this Stokes number increases with increasing droplet diameter, and the large Stokes numbers observed for large droplets are responsible for the slow rate of decrease in droplet velocity behind the OSW in those cases, as shown in figure 8(g). Therefore, this Stokes number can provide a qualitative reference for the behaviour of droplets suspended in the gas flow behind the OSW.
Figure 14. (a) Stokes numbers and (b) Damköhler numbers along the streamline near the wedge for 
$\phi _{l} = 0.5$ with different droplet diameters.
The Damköhler number for a particular diameter is defined as 
$Da = {\bar {t}_i} / {\bar {t}_e},$ where 
$\bar {t}_e$ is the time-averaged main flow time of droplet evaporation, which is the time from when the droplets first pass through the OSW to when they are completely evaporated. As shown in figure 14(b), the Damköhler number is observed to decrease as the droplet diameter increases. The cases with unsteady ODWs have large Damköhler numbers, indicating that the unsteady behaviour occurs when the chemical reaction time is greater than the droplet evaporation time, and vice versa. The Damköhler numbers can provide a qualitative reference for the steadiness of ODWs.
3.5. Criterion for identifying unsteady and steady ODWs
Steady combustion is essential to ODE application, and knowing the mechanisms of steady and unsteady ODWs is important. For two-phase oblique detonation with large droplets, we find that the occurrence of unsteady and steady detonation propagation is mainly related to the location of droplet evaporation. Figure 15 shows the schematic of the evaporation layer and initiation length in a two-phase ODW to explain quantitatively why the initiation reaction finishes before droplet evaporation when the droplet diameter is large. The mechanisms of unsteady and steady two-phase ODWs with large droplets that evaporate behind the OSW and ODW can be clarified with the following criterion:
\begin{equation} \left.\begin{aligned}\delta_{E} < L_{E}\ & \text{unsteady ODW, } \\ \delta_{E} > L_{E}\ & \text{steady ODW. }\end{aligned}\right\} \end{equation}
Here, 
$\delta _{E}$ is the thickness of the droplet evaporation layer, which is defined along the streamline close to the wedge, from the start of the induction OSW to the end of droplet evaporation, where the evaporation rate is zero.
Figure 15. Schematics of evaporation layer and initiation length in the two-phase ODW.
Here 
$L_{E}$ is the corresponding theoretical initiation length, which is calculated using the constant-volume combustion (CVC) calculation (Teng et al. Reference Teng, Ng and Jiang2017). Unlike a previous study (Teng et al. Reference Teng, Ng and Jiang2017), this study used the average temperature, pressure and equivalence ratio in the evaporation layers along the streamline close to the wedge as inputs to calculate the induction times according to CVC. For large droplet diameters, the parameters in the induction zone was used to avoid the effect of reaction heat release, as shown in table 2. The induction time is defined as the reaction time required to attain a mixture temperature that is 10 % higher than the average temperature of the evaporation layer. The corresponding theoretical initiation length 
$L_{E}$ can be deduced by multiplying induction time with velocity of the evaporation layer. It is worth noting that to evaluate whether a postshock cooled gas can be ignited within the evaporation layer, the input equivalence ratio is considered to be the equivalence ratio of an inflow gas mixture due to unfinished droplet evaporation.
Table 2. Average parameters of the evaporation layer thickness 
$\delta _{E}$, and the corresponding theoretical initiation length 
$L_{E}$ in cases with 
$\phi _{l} = 0.5$ and 
$d_{0} = 10$, 20 and 
$30\ \mathrm {\mu } {\rm m}$.
The numerical and theoretical initiation lengths are obtained from the parameters of the streamline near the wedge according to previous studies (Teng et al. Reference Teng, Ng and Jiang2017, Reference Teng, Bian, Zhou and Zhang2021a). They agree well, indicating that ODW initiation is directly related to the induction reaction at the bottom, and ODW induction near the wedge is similar to that of CVC. Considering that the steadiness of liquid ODWs is a result of both induction reaction and droplet evaporation, it is reasonable to use the streamline parameters for this criterion. Furthermore, this unsteady behaviour is mostly caused by the hot spots near the wedge, as shown in figures 6 and 7. The streamline selection is based on this physical phenomenon.
The postshock pressure and velocity are slightly affected by droplet evaporation. Therefore, 
$L_E$ is mainly related to the temperature of the evaporation layer. Figure 16 shows the average temperature 
$T_E$ of the evaporation layers along the streamlines near the wedge for 
$\phi _{l} = 0.5$ with different initial droplet diameters. The temperature of the evaporation layer increases when the droplet diameter increases from 5 to 
$25\ \mathrm {\mu } {\rm m}$ because of the increase of preshock temperature (discussed in § 3.3 and shown in figure 10). Notably, there is a sudden increase in 
$T_E$ when the droplet diameter changes from 25 to 
$30\ \mathrm {\mu } {\rm m}$. This can be attributed to the local heat release rate of the reaction, which surpasses the droplet heat transfer rate (as discussed in § 3.4 and illustrated in figure 10). To better assess the ignition ability of gases in the evaporation layer and avoid the effects of heat release, the induction temperature, which is the temperature behind the induction OSW in figure 10(di), is chosen as 
$T_E$ in cases where the heat release reaction occurs before the evaporation behind the OSW. As shown in table 2, 
$L_E$ decreases with increasing 
$T_E$ as droplet diameter increases.
Figure 16. Average temperature 
$T_E$ of the evaporation layers along the streamlines near the wedge for 
$\phi _l = 0.5$ with different initial droplet diameters.
In table 2, 
$\delta _{E}$ increases first and then decreases. Here 
$\delta _{E}$ increases when the droplet diameter changes from 
$10\ \mathrm {\mu } {\rm m}$ to 
$20\ \mathrm {\mu } {\rm m}$ owing to the increased droplets mass transfer behind the OSW. Because heat is released before significant heat is transferred by the big droplets (see figure 10) and the increase of in temperature accelerates droplet evaporation, 
$\delta _{E}$ decreases when the droplet diameter changes from 
$20\ \mathrm {\mu } {\rm m}$ to 
$30\ \mathrm {\mu } {\rm m}$.
The variation range of 
$\delta _{E}$ with diameter is smaller than that of 
$L_E$ in table 2. When the droplet diameter is small, the characteristic droplet evaporation thickness 
$\delta _{E}$ is shorter than the characteristic initiation reaction length 
$L_E$. The droplets completely evaporate before the ODW is initiated, and the ODW is unsteady. Similarly, when the droplet diameter is large, 
$\delta _{E}$ is larger than 
$L_E$, and droplet evaporation occurs after ODW initiation and has little effect on the induction reaction, and the ODW initiation structure is steady.
Figure 17 shows the comparison of the evaporation layer thickness 
$\delta _{E}$ and the corresponding theoretical initiation length 
$L_E$ for 
$\phi _{l} = 0.3$, 0.5 and 0.7 with different initial droplet diameters to validate the criterion for distinguishing steady and unsteady ODWs. Despite the droplet percentages, 
$L_E$ decreases and 
$\delta _{E}$ increases with increasing droplet diameter in most cases. An unsteady ODW occurs with small droplet diameter when 
$L_E$ is larger than 
$\delta _{E}$. When 
$L_E$ is smaller, a steady ODW occurs. Therefore, unsteady and steady detonation propagation in an ODW can be identified by comparing the droplet evaporation layer thickness 
$\delta _{E}$ and the corresponding initiation length 
$L_E$. Moreover, the fluctuating range of 
$\delta _{E}$ is smaller than that of 
$L_E$. Therefore, the steadiness of the two-phase ODW is mainly contribute to 
$L_E$.
Figure 17. Comparison of the thickness of the droplet evaporation layer 
$\delta _{E}$ and corresponding theoretical initiation length 
$L_{E}$ for 
$\phi _{l} = 0.3$, 0.5 and 0.7 with different initial droplet diameters (
$\phi _{l}+ \phi _{g}=1.0$).
Notably, the diameter boundary of a steady ODW decreases with increasing droplet percentage 
$\phi _{l}$, as shown in figure 17, where the sum of the droplet and gas percentages is equal to 1.0. For example, the ODW with 
$d_{0} = 25\ \mathrm {\mu } {\rm m}$ is unsteady when 
$\phi _{l} = 0.3$ and steady when 
$\phi _{l} = 0.7$. When the gas equivalence ratio is decreased, the ODW can be initiated in a shorter distance according to previous study (Teng et al. Reference Teng, Bian, Zhou and Zhang2021a). This issue is also addressed by the pure-gas cases in this study. In figure 2, the initiation length for 
$\phi _{g} = 0.5$ and 
$\phi _{l} = 0.0$ is approximately 0.023 m, which is shorter than the initiation length for 
$\phi _{g} = 1.0$ and 
$\phi _{l} = 0.0$, which is approximately 0.044 m in figure 3. Here 
$L_E$ values are compared for 
$d_{0} = 25\ \mathrm {\mu } {\rm m}$ with 
$\phi _{l} = 0.3$ and 0.7 in figures 17(a) and 17(c). For 
$\phi _{l} = 0.7$, 
$L_E$ is shorter because of the lower equivalence ratio of the gas mixture deduced from the gas percentage (
$\phi _{g} = 0.3$ when 
$\phi _{l} = 0.7$). Moreover, 
$\delta _{E}$ changes little with different droplet percentages compared with 
$L_E$, and the diameter boundary of the steady ODW decreases with increasing droplet percentage 
$\phi _{l}$ mainly owing to the reduction in 
$L_E$.
Because of the difficulties in conducting ODW experiments, no corresponding experimental data are available to validate this criterion. Cases with droplet percentages from 20 % to 80 % and droplet diameters from 5 to 
$40\ \mathrm {\mu } {\rm m}$ were tested to validate this criterion over a wider range. The black and red squares in figure 18 represent the unsteady and steady ODW, respectively. The dashed line is the boundary between steady and unsteady ODWs predicted by this criterion, where 
$L_E$ is equal to 
$\delta _{E}$. The results agree well, demonstrating that this criterion is also valid. Moreover, as discussed above, the diameter boundary of a steady ODW decreases with increasing droplet percentage 
$\phi _{l}$ mainly owing to the reduction in 
$L_E$.
Figure 18. Validation of this method for different droplet percentages (
$\phi _{l}$) and different droplet diameters (
$d_{0}$), where (
$\phi _{l}+ \phi _{g}=1.0$).
4. Conclusions
Two-dimensional ODWs in partially prevaporized n-heptane sprays were simulated with the Eulerian–Lagrangian method using a skeletal chemical mechanism. The emphasis was laid on the influences of initial droplet diameter and liquid equivalence on oblique detonation. The results show that the initial droplet diameter significantly affects ODW propagation in a partially prevaporized n-heptane spray. When the initial droplet diameter is less than a certain value, the droplet evaporation in the initiation zone causes the NDW to appear within a certain distance, move backwards, disappear and re-form, therefore leading to unsteady NDW propagation. However, when the initial droplet diameter is further increased beyond the critical value, a steady ODW is produced, and the second reaction front behind the ODW becomes an intense deflagration front.
Further analysis suggests that the ODW is unsteady when the evaporation layer thickness along the streamline is larger than the corresponding theoretical initiation length. A steady ODW is obtained when the evaporation layer thickness along the streamline is less than the corresponding theoretical initiation length. Moreover, the initiation length increases first and then decreases with increasing initial droplet diameter.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.194.
Funding
This research was supported by the National Natural Science Foundation of China (nos. 12202014 and 12002041) and Advanced Jet Propulsion Innovation Center, Aero Engine Academy of China (project ID. HKCX2022-01-018).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Study on mesh and particle number dependences
Figures 19(a) and 20(a) show the qualitative comparison of temperature fields with the present and fine meshes for the pure-gas case with 
$\phi _{g} = 1.0$ and 
$\phi _{l} = 0.0$ and the gas–liquid two-phase case with 
$\phi _{l} = 0.5$, 
$\phi _{g} = 0.5$, and 
$d_{0} = 30\ \mathrm {\mu } {\rm m}$, respectively. The grid sizes of the fine and coarse meshes in figures 19 and 20 are 50 and 
$100\ \mathrm {\mu } {\rm m}$, respectively. In addition, figures 19(b)–19(c) and 20(b)–20(c) show the quantitative comparisons by plotting the temperature and pressure along three typical streamlines, 
${y} = 0.02$ m, 0.04 m and 0.06 m in figures 19(a) and 20(a), with different mesh resolutions. These lines correspond to different flow regions of the ODW field, including the wedge surface, secondary ODW structures and steady ODW surface. Insignificant differences are observed between the results of the present and fine meshes. Therefore, the chosen mesh resolution of 
$100\ \mathrm {\mu } {\rm m}$ provides converged global structures that are sufficient to guarantee the reliability of the conclusions in this study. Moreover, figure 21 presents a similar qualitative and quantitative comparison for different particle numbers. The negligible differences observed indicate that the results of this study are not influenced by particle number.
Figure 19. Distributions of (a) temperature with different mesh resolutions for pure gas case with 
$\phi _{g} = 1.0$ and 
$\phi _{l} = 0.0$. White lines are streamlines at 
${y} = 0.02$ m, 0.04 m and 0.06 m. Profiles of (b) temperature and (c) pressure with different mesh resolutions.
Figure 20. Distributions of (a) temperature with different mesh resolutions for 
$\phi _{l} = 0.5$, 
$\phi _{g} = 0.5$ and 
$d_{0} = 30\ \mathrm {\mu } {\rm m}$. White lines are streamlines at 
${y} = 0.02$ m, 0.04 m and 0.06 m. Profiles of (b) temperature and (c) pressure with different mesh resolutions.
Figure 21. Distributions of (a) temperature with different mesh resolutions for 
$\phi _{l} = 0.5$ and 
$\phi _{g} = 0.5$ and 
$d_{0} = 30\ \mathrm {\mu } {\rm m}$. White lines are streamlines at 
${y} = 0.02$ m, 0.04 m and 0.06 m. Profiles of (b) temperature and (c) pressure with different particle numbers.
Appendix B. Effect of boundary layers
Figure 22 shows the ODWs simulated with 
$\phi _{l} = 0.5$, 
$\phi _{g} = 0.5$, and 
$d_{0} = 2$, 10 and 
$30\ \mathrm {\mu } {\rm m}$ to clarify the effect of the boundary layers. The no-slip reflecting boundary condition is used on the wedge surface, and the mesh near the wall is densified to capture the boundary layers. The wall mesh is densified. According to a previous study (Sontheimer et al. Reference Sontheimer, Kronenburg and Stein2021), the mesh size should be larger than local droplet size. Therefore, the sizes of wall mesh in the y-direction for 
$d_{0} = 2$, 10 and 
$30\ \mathrm {\mu } {\rm m}$ are chosen to be 10, 20 and 
$50\ \mathrm {\mu } {\rm m}$, respectively. The size is increased outward with an equal ratio of 1.05 until it reaches the same scale as the outer mesh, which has sizes of 200, 400 and 
$100\ \mathrm {\mu } {\rm m}$, respectively. Previous study (Fang et al. Reference Fang, Zhang and Hu2019) reports that despite the recirculation zone in the vicinity of the wall, the type of global ODW configuration with 
$d_{0} = 2$, 10 and 
$30\ \mathrm {\mu } {\rm m}$ remains unaffected. Moreover, the ODWs remain unsteady for 
$d_{0} = 2$ and 
$10\ \mathrm {\mu } {\rm m}$ and steady for 
$d_{0} = 30\ \mathrm {\mu } {\rm m}$. The unsteady behaviours of ODWs remain the same, including the explosion points being generated and NDWs re-forming and moving downstream, as shown in figures 22(a) and 22(b). The initiation lengths in most cases considering boundary layers decrease owing to the high temperature caused by the recirculation zone. Furthermore, the initiation lengths increase for small droplet diameters in figures 22(a) and 22(b), where the initiation lengths for 
$d_{0} = 2$ and 
$10\ \mathrm {\mu } {\rm m}$ are 0.12 and 0.20 m, and decrease for large droplet diameters in figure 22(c), where the initiation length for 
$d_{0} = 30\ \mathrm {\mu } {\rm m}$ is approximately 0.06 m. Generally, most conclusions in this study are still valid.
Figure 22. Temperature distributions considering boundary effects with 
$\phi _{l} = 0.5$, 
$\phi _{g} = 0.5$ and (a) 
$d_{0} = 2\ \mathrm {\mu } {\rm m}$, (b) 
$d_{0}= 10\ \mathrm {\mu } {\rm m}$ and (c) 
$d_{0} = 30\ \mathrm {\mu } {\rm m}$.
