## 2 Model Description  
We first present the constraints on the combustion models. Then, we present models that incorporate PDFs that depend only on the maximal radius. Several standard assumptions are commonly used for models that describe the thermal explosion of vaporised fuel droplets. Our model is considered adiabatic, because ignition processes occur faster than heat losses by diffusion processes. The pressure variations were assumed to be negligible [21]. In addition, we assume that the diffusion of the gas phase is insignificant compared to that of the liquid phase, and the heat transfer coefficient is determined by the gas features [25]. We further assume that the liquid temperature is identical to the temperature of the saturated liquid because a quasi-steady-state approximation is reasonable for vaporising droplets [26]. Finally, the reaction is modelled in the simplest possible way as a first-order, highly exothermic chemical reaction. This means that the model is represented as a system of nonlinear ordinary differential equations. Here, we present the mathematical model of polydisperse fuel spray, and the discrete model (more details about the assumption of the model can be found in [23]):  
\[\alpha_{g}\rho_{g}c_{pg}\frac{dT_{g}}{dt}=\alpha_{g}\mu_{f}Q_{f}AC_{f}e^{(- \frac{E}{Bt_{g}})}-4\pi\lambda_{g}(T_{g}-T_{d})\sum_{i=1}^{m}R_{d_{i}}n_{d_{i}}, \tag{1}\]  
\[\frac{d(R_{di}^{2})}{dt}=\frac{2\lambda_{g}}{\rho_{L}L}(T_{d}-T_{g}),\;\;i=1,\dots,m,\;\;(m\;\;equations) \tag{2}\]  
\[\frac{dC_{f}}{dt}=-AC_{f}e^{(-\frac{E}{Bt_{g}})}+\frac{4\pi\lambda_{g}}{L \alpha_{g}\mu_{f}}(T_{g}-T_{d})\sum_{i=1}^{m}R_{d_{i}}n_{d_{i}}. \tag{3}\]  
The initial conditions are  
\[at\;\;t=0:\;\;T_{g}(t=0)=T_{g0},\;\;R_{d_{i}}=R_{d_{i0}},\;\;\forall i:T_{d_{i }}=T_{g0},\;\;C_{f}=C_{f0}. \tag{4}\]  
For the model above, the droplet radius distribution can be approximated by a smooth PDF. The droplet radii sum is replaced by the integral of the corresponding PDF [17; 24]. The motivation for this approximation is based on a standard representation of the Riemann integral as a limit of the Riemann sum. This means that this approximation works better for a high number of parcels and can be less accurate for small numbers of parcels. Formally, we can use any smooth function that satisfies the equality. \[\int_{0}^{\infty}R_{0}\hat{P}_{0}(R_{0})dR_{0}=\sum_{i=1}^{m}n_{d_{i}}R_{d_{i 0}}. \tag{5}\]  
Here, we use the sub-index \(0\) for \(t=0\). We use commonly known distributions such as the Rosin-Rammler, Nukiyama-Tanasawa, and gamma distributions as examples. The corresponding normalised PDF is given by the following equations:  
\[P_{0}(R_{0})=\frac{\hat{P}_{0}(R_{0})}{\langle\hat{P}_{0}\rangle},\;\;\text{ {\em were}};\;\;\langle\hat{P}_{0}\rangle\equiv\int_{0}^{\infty}\hat{P}_{0}(R_{0 })dR_{0}. \tag{6}\]  
Note that the integration by variable \(R_{0}\) is an approximation that corresponds to the sum of all the radii of the droplets at the initial time. The droplet estimation by the PDF is superior to the 'parcel' approach, because the smooth PDF allows us to extrapolate experimental data for radii that do not appear in the standard measurement. In addition, the PDF is not affected by the limitations of the parcel method. The parcel method describesthe droplet distribution poorly, because it usually uses a comparatively small number of parcels due to calculation complexity. Using (2) conditions, the radius \(R\) can be expressed as follows:  
\[R=\left(R_{m}^{2}-R_{m,0}^{2}+R_{0}^{2}\right)^{(1/2)}. \tag{7}\]  
Thus, \(R_{0}\) can be viewed as a function of \(R\). \[R_{0}=G(R)\equiv\left(R^{2}-(R_{m}^{2}-R_{m,0}^{2})\right)^{(1/2)}. \tag{8}\]  
We define the rectangle function:  
\[\Pi_{a}(x)\equiv U_{0}(x)-U_{a}(x). \tag{9}\]  
Here, \(U_{0}(x),U_{a}(x)\) are unit step functions.