2. Statistical model of transition

Consider a steady continuous compressible medium in which solid spherical particles of radius $r_{p}^{*}$ are suspended with an average volumetric concentration $c_{{\infty }}^{\mathrm{*}}$ (an extension to the case of particles distributed in radius is given below). Hereafter, notations are supposed dimensionless unless denoted with a superscript asterisk. The subscripts ‘w’, ‘ $\infty$ ’, ‘e’, ‘p’ and ‘с’ indicate quantities on the wall, in the free stream, on the upper edge of the boundary layer, related to particulates and to the point of particle–surface collision, respectively. To make quantities dimensionless, the following scales are used: length $L^{\mathrm{*}}=1$ m, time $\tau ^{\mathrm{*}}=L^{\mathrm{*}}/U_{e}^{\mathrm{*}}$ , speed $U_{e}^{\mathrm{*}}$ , temperature $T_{e}^{\mathrm{*}}$ and pressure $\rho _{e}^{\mathrm{*}}U_{e}^{\mathrm{*}2}$ . The particles are still relative to the medium and randomly distributed in the volume with a constant distribution function. A sharp wedge of infinite span and a half-angle $\theta$ moves through the medium at zero angle of attack at a supersonic speed. The Cartesian coordinate systems $Or_{1}r_{2}r_{3}$ and $Oxyz$ are shown in figure 1. In the system $Or_{1}r_{2}r_{3}$ , the free-stream velocity is $\mathbf{U}_{{\infty }}^{\mathrm{*}}=(U_{{\infty }}^{\mathrm{*}},0,0)$ , and the Mach number is $M_{{\infty }}=U_{{\infty }}^{\mathrm{*}}/a_{{\infty }}^{\mathrm{*}}\gt 1$ .

Particulates penetrate from the free stream into the shock layer and collide with the wedge surface. Crossing the laminar boundary layer, they excite wave packets associated with the first or second mode according to the terminology of Mack (Reference Mack1969). The particles disappear after collision, i.e. the collisions are perfectly inelastic and their traces on the wall are negligible. The wave packets grow downstream and form a random disturbance field. According to the amplitude method of Mack (Reference Mack1977) (see also Fedorov & Tumin Reference Fedorov and Tumin2021), the transition onset occurs at $x=x_{t}$ where the root-mean-square (RMS) value of fluctuations (hereafter we deal with the wall pressure fluctuations, $p_{w}^{\prime}$ ) reaches a critical level. The problem is to find $x_{t}$ , which is further referred to as a transition onset station.

It is assumed that:

1. (i) Each particle excites a wave packet independently from other particles. The packets have a low initial intensity and propagate downstream without interacting with each other (the linear theory is used). In this case, the average frequency of particle collisions per unit area, $\nu ^{*}$ of dimension s⁻¹ m⁻², equals the average frequency of the wave-packet birth in the boundary layer. Since the undisturbed flow does not depend on the spanwise coordinate z, we get

   (2.1) \begin{equation}\nu ^{*}\left(x\right)=c_{\infty }^{*}U_{\infty }^{*}\sin \theta _{c}\left(x\right),\end{equation}
   where $\theta _{c}$ is the angle between the $Ox$ direction (figure 1) and the particle velocity vector at the collision point. This assumption is justified by the fact that, as a rule, the concentration of particles of required radius is low.

2. (ii) Collisions of particles with the wedge surface are random; i.e. the excited wave packets are not coherent. In this case, the average squared amplitude of disturbances in the far field ( $x\gg x_{c}$ ) equals the sum of squared amplitudes of individual wave packets.

3. (iii) Initially, wave packets are localized in $t$ and $z$ near the collision points $(t_{c},z_{c},x_{c})$ .

To solve the formulated problem, it is necessary to:

1. (i) Solve the receptivity and linear stability problems to obtain the dependence of wave-packet amplitude on $x$ . This was done theoretically by Fedorov (Reference Fedorov2013) and numerically by Chuvakhov et al. (Reference Chuvakhov, Fedorov and Obraz2019), Al Hasnine et al. (Reference Al Hasnine2020), Browne at al. (Reference Browne, Al Hasnine and Brehm2021) and Habeck et al. (Reference Habeck2024).

2. (ii) Calculate an average statistical contribution of all wave packets to the disturbance field at the observation stations $x\gt x_{c}$ :

   (2.2) \begin{equation}\left\langle p_{w}^{\prime 2}\right\rangle \left(x\right)\equiv \lim _{\begin{array}{l} T\rightarrow {\infty }\\ Z\rightarrow {\infty } \end{array}}\left[\frac{1}{T\cdot Z}\int _{-T/2}^{T/2}{\rm d}\tilde{t}\int _{-Z/2}^{Z/2}{\rm d}\tilde{z}p_{w}^{\prime 2}(\tilde{t},\tilde{z},x)\right].\end{equation}

3. (iii) Formulate the amplitude criterion $p_{w,cr}^{\prime}$ and find the transition onset station using the relationship $\sqrt{\langle p_{w}^{\prime2}\rangle (x_{t})}=p_{w,cr}^{\prime}$ .

A disturbance propagating in a weakly non-parallel boundary layer can be expanded in the biorthogonal eigenfunction system as a sum of normal modes (Fedorov & Tumin Reference Fedorov and Tumin2003; Tumin Reference Tumin2007). Following Fedorov (Reference Fedorov2013) and Chuvakhov et al. (Reference Chuvakhov, Fedorov and Obraz2019), we consider the case of $M_{\infty }=4$ , $\theta =14^{\circ}$ . Since the local Mach number at the upper boundary-layer edge is $M_{e}=3$ , only the first mode is unstable. Its wave packet propagates from the collision point $(x_{c},z_{c},t_{c})$ and is expressed as

(2.3) \begin{equation}p_{w}^{\prime}\left(t,z,x,t_{c},z_{c},x_{c}\right)=\frac{1}{\left(2\pi \right)^{2}}\int _{-{\infty }}^{{\infty }}{\rm d} \omega \int _{-{\infty }}^{{\infty }}{\rm d}\beta A\left(\omega ,\beta ;x_{c},x\right){\rm e}^{i\beta (z{-z_{c}})-i\omega (t{-t_{c}})},\end{equation}

(2.4) \begin{equation}A(\omega ,\beta ,x_{c},x)=\varepsilon (\omega ,\beta ,x_{c}){\rm e}^{iS({x_{c}},x,\omega ,\beta )},\end{equation}

(2.5) \begin{equation}S\left(\omega ,\beta ,x_{c},x\right)=\int _{x_{c}}^{x}\alpha \left(x^{\prime},\omega ,\beta \right){\rm d}x^{\prime},\end{equation}

where $\omega$ and $\beta$ are real frequency and spanwise wavenumber, $\alpha =\alpha _{r}+{\rm i}\alpha _{i}$ is complex eigenvalue, $S$ is eikonal and $A$ is spectral density of the wave packet observed at station $x$ . The function $\varepsilon (\omega ,\beta ,x_{c})$ results from a theoretical or numerical solution of the receptivity problem and determines the wave-packet spectral density at the collision point:

(2.6) \begin{equation}\varepsilon (\omega ,\beta ,x_{c})=\int _{-{\infty }}^{{\infty }}\int _{-{\infty }}^{{\infty }}p^{\prime}_{w}(\tilde{t},\tilde{z},x_{c})e^{i\omega \tilde{t}-i\beta \tilde{z}}{\rm d}\tilde{t}{\rm d}\tilde{z}, \tilde{t}=t-t_{c}, \tilde{z}=z-z_{c}.\end{equation}

Hereafter we use the theoretical solution of Fedorov (Reference Fedorov2013). If the eigenfunction $\hat{p}_{m}(\omega ,\beta ,x,y)$ of the unstable mode $m$ is normalized to the static pressure $P_{e}$ in the undisturbed boundary layer $(\hat{p}_{m}(\omega ,\beta ,x,y=0)/P_{e}=1)$ , then, in notations of Fedorov (Reference Fedorov2013), the spectral density is $\varepsilon =r_{p}^{2}D_{m}(\omega ,\beta ,x_{c})$ , where $D_{m}$ is calculated using equations (3.35) and (3.36) of Fedorov (Reference Fedorov2013) for the first mode ( $m=1$ ).

By assumption 2, the mean-square value of pressure fluctuations at observation station $x$ is obtained by averaging over $t$ and $z$ the corresponding squared amplitudes of all wave packets excited by particles on the surface strip $0\lt x_{c}\lt x$ , $-Z/2\lt z_{c}\lt Z/2$ over the time interval $-T/2\lt t_{c}\lt T/2$ , where $Z\rightarrow {\infty }, T\rightarrow {\infty }$ . Thus,

(2.7) \begin{align}\left\langle p_{w}^{\prime 2}\right\rangle &\equiv \lim _{\begin{array}{l} T\rightarrow {\infty }\\ Z\rightarrow {\infty } \end{array}}\left[\frac{1}{T\cdot Z}\int _{0}^{x}\left[\nu (x_{c})\cdot T\cdot Z\int _{-T/2}^{T/2}{\rm d}\tilde{t}\int _{-Z/2}^{Z/2}{\rm d}\tilde{z}\left(p_{w}^{\prime 2}\right)\right]{\rm d}x_{c}\right]\\ &=\int _{0}^{x}{\rm d}x_{c}\left[\nu (x_{c})\int _{-{\infty }}^{{\infty }}\int _{-{\infty }}^{{\infty }}p_{w}^{\prime 2}(\tilde{z},\tilde{t})\,{\rm d}\tilde{t}\,{\rm d}\tilde{z}\right].\nonumber\end{align}

The same result can be obtained using an extension of Campbell’s theorem to random space–time processes (Cohen Reference Cohen2014), and considering a process that consists of a sequence of pulses $p(z,t)=\sum _{k=1}^{{\infty }}v(z-z_{c},t-t_{c})$ , where random times $t_{c}$ and coordinates $z_{c}$ are Poisson distributed with an average sample rate ${\rm d}\gamma =\nu (x_{c}){\rm d}x_{c}.$

Using Parseval’s theorem

(2.8) \begin{equation}\int _{-{\infty }}^{{\infty }}\int _{-{\infty }}^{{\infty }}p_{w}^{\prime 2}(\tilde{z},\tilde{t})\,{\rm d}\tilde{t}\,{\rm d}\tilde{z}=\frac{1}{(2\pi )^{2}}\int _{-{\infty }}^{+{\infty }}\int _{-{\infty }}^{+{\infty }}\left| A(\omega ,\beta ,x_{c})\right| ^{2}\,{\rm d}\omega\, {\rm d}\beta,\end{equation}

and substituting the amplitude (2.4) into (2.8), we reduce (2.7) to the form

(2.9) \begin{equation}\left\langle p_{w}^{\prime 2}\right\rangle (x)=\frac{1}{(2\pi )^{2}}\int _{0}^{x}{\rm d}x_{c}\left[\nu (x_{c})\int _{-{\infty }}^{+{\infty }}\int _{-{\infty }}^{+{\infty }}\left| \varepsilon (\omega ,\beta ,x_{c})\right| ^{2}{\rm e}^{2N({x_{c}},x,\omega ,\beta )}\,{\rm d}\omega \, {\rm d}\beta \right],\end{equation}

(2.10) \begin{equation}N(x_{c},x,\omega ,\beta )=-S_{i}(x_{c},x,\omega ,\beta )=-\int _{x_{c}}^{x}\alpha _{i}\left(x^{\prime},\omega ,\beta \right){\rm d}x^{\prime},\end{equation}

where $N$ is an integral amplification factor of elementary wave $(\omega ,\beta )$ . The functions $N$ and $\alpha$ are calculated using the linear stability theory. Such calculations can be carried out in the local-parallel approximation using a standard stability code. The non-parallel mean flow effects can be taken into account using a code for solving the parabolic stability equations (see e.g. Chang Reference Chang2003) or the method of multiple scales (Gaponov Reference Gaponov1980; Nayfeh Reference Nayfeh1980; Tumin & Fedorov Reference Tumin and Fedorov1982).

The integral (2.9) can be computed numerically. However, for typical flight conditions the values of $N$ are large at the transition onset: $N\gt 10$ according to data reported by Malik (Reference Malik, Bushnell and Hefner1990). Therefore, the integrals in (2.9) can be estimated asymptotically. Assuming that $N\gg 1$ at observation station $x$ , we apply the saddle point method. Since the integrand is real-valued, (2.9) is cast in the form

(2.11) \begin{equation}\left\langle p_{w}^{\prime 2}\right\rangle (x)=\frac{1}{2\pi ^{2}}\int _{0}^{x}{\rm d}x_{c}\left[\nu (x_{c})\int _{0}^{+{\infty }}{\rm d}\omega \int _{-{\infty }}^{+{\infty }}\left| \varepsilon (\omega ,\beta ,x_{c})\right| ^{2}{\rm e}^{2N({x_{c}},x,\omega ,\beta )}\,{\rm d}\beta \right].\end{equation}

The saddle points are determined from the conditions $N_{\beta }=0, N_{\omega }=0$ . For the first mode, there are two saddle points $(\omega _{s},\pm \beta _{s})$ , $\omega _{s}\gt 0, \beta _{s}\gt 0$ , which give equal contributions to the integral. Then

(2.12) \begin{equation}\left\langle p_{w}^{\prime 2}\right\rangle (x)=\frac{1}{\pi }\int _{0}^{x}{\rm d}x_{c}\left[\nu (x_{c})\frac{\left| \varepsilon (\omega _{s},\beta _{s},x_{c})\right| ^{2}}{\sqrt{N_{\omega \omega }N_{\beta \beta }-N_{\omega \beta }^{2}}}{\rm e}^{2{N_{s}}({x_{c}},x)}\right],\end{equation}

where $N_{s}=N(\omega _{s},\beta _{s})$ , and the second derivatives $N_{\omega \omega }, N_{\beta \beta }\, \unicode{x0438}\, N_{\omega \beta }$ are calculated at $(\omega ,\beta )=(\omega _{s},\beta _{s})$ . The integral (2.12) is also estimated using the saddle point method. The saddle point, $x_{cs}$ , is obtained from the equation

(2.13) \begin{equation}\frac{\partial N_{s}\left(x_{cs}\right)}{\partial x_{c}}=\alpha _{i}\left(x_{cs},\omega _{s},\beta _{s}\right)=0.\end{equation}

Hence, $x_{cs}=x_{0}(\omega _{s}, \beta _{s})$ is a neutral point of the dominant wave $(\omega _{s}, \beta _{s})$ , and

(2.14) \begin{equation}\left\langle p_{w}^{\prime2}\right\rangle (x)=\left(\nu \left(x_{c}\right)\frac{\left| \varepsilon (\omega _{s},\beta _{s},x_{c})\right| ^{2}}{\pi \sqrt{N_{\omega \omega }N_{\beta \beta }-N_{\omega \beta }^{2}}}\sqrt{\frac{\pi }{N_{{x_{c}}{x_{c}}}}}{\rm e}^{2{N_{s}}\left(x_{c},x\right)}\right)_{{x_{c}}={x_{0}}\left(\omega _{s},\beta _{s}\right)}.\end{equation}

The physical meaning of (2.14) is that the major contribution to the disturbance intensity at observation station $x$ is given by the wave packets excited in the narrow strip $(x_{0}-\Delta x_{0}/2, x_{0}+\Delta x_{0}/2)$ , where $\Delta x_{0}=\sqrt{\pi /N_{xx}(\omega _{s},\beta _{s},x_{0})}$ . The spectral bandwidth of these packets is $\Delta \omega \Delta \beta =\pi /\sqrt{N_{\omega \omega }N_{\beta \beta }-N_{\omega \beta }^{2}}$ .

In the second-mode case, the most unstable are plane waves of $\beta _{s}=0$ . There is only one saddle point, $\omega _{s}$ , and $N_{\omega \beta }\equiv 0$ , and the wave-packet intensity is estimated as

(2.15) \begin{equation}\left\langle p_{w}^{\prime2}\right\rangle (x)=\left(\nu \left(x_{c}\right)\frac{\left| \varepsilon (\omega _{s},0,x_{c})\right| ^{2}}{2\pi \sqrt{N_{\omega \omega }N_{\beta \beta }}}\sqrt{\frac{\pi }{N_{{x_{c}}{x_{c}}}}}{\rm e}^{2{N_{s}}\left(x_{\mathrm{c}},x\right)}\right)_{{x_{c}}={x_{0}}\left(\omega _{s},0\right)}.\end{equation}

The corresponding amplitude criterion is presented by Marineau (Reference Marineau2017). Note that this solution has not yet been tested using direct numerical simulation. As shown by Al Hasnine et al. (Reference Al Hasnine, Russo, Tumin and Brehm2023b), in high-speed boundary layers, where the second mode is most unstable, the interaction between first and second modes due to non-parallel effects may affect the initial spectrum of the second-mode instability. This case requires further theoretical analysis.

Thus, the foregoing statistical model leads to the compact analytical solution (2.14) for the first mode or (2.15) for the second mode. These solutions allow for quick calculations of the RMS fluctuations, $p^{\prime}_{w,\mathrm{RMS}}(x)=\sqrt{\langle p_{w}^{\prime2}\rangle }$ , in the far field where $N$ -factors are sufficiently large.

Since the dependence on the particle radius is contained only in the collision frequency $\nu$ and the spectral density $\varepsilon$ , the result is easily extended to a more realistic case where the particle concentration $c_{\infty }$ is distributed in radius $r_{p}$ . Introducing $g(r_{p})={\rm d}c_{\infty}(r_p)/{\rm d}r_p$ , the collision frequency for particles having the radius from the interval $(r_{p},r_{p}+ {\rm d}r_{p})$ is

(2.16) \begin{equation}{\rm d}\nu (x,r_{p})=U_{\infty }\sin \theta_{c}\left(x,r_{p}\right)g(r_{p}){\rm d}r_{p}\equiv F(x,r_{p}){\rm d}r_{p}.\end{equation}

Replacing $\nu (x_{0})$ in (2.14) by ${\rm d}\nu (x_{0},r_{p})$ and integrating over $r_{p}$ , we obtain

(2.17) \begin{equation}\left\langle p_{w}^{\prime 2}\right\rangle (x)=\frac{\varPhi (\omega _{s},\beta _{s},x_{0})}{\pi \sqrt{({N_{\omega \omega }}{N_{\beta \beta }}-N_{\omega \beta }^{2})_{{x_{c}}={x_{0}}}}}\sqrt{\frac{\pi }{N_{{x_{c}}{x_{c}}}(x_{0},x)}} {\rm e}^{2{N_{s}}({x_{0}},x)},\end{equation}

(2.18) \begin{equation}\Phi \left(\omega _{s},\beta _{s},x_{0}\right)=\int _{r_{p,min}}^{r_{p,max}}F(x_{0},r_{p})\left| \varepsilon (\omega _{s},\beta _{s},x_{0},r_{p})\right| ^{2}{\rm d}r_{p},\end{equation}

where $r_{p,min }$ is the minimum radius of particles penetrating into the boundary layer at station $x_{0} $ ; $r_{p,max }$ is the maximum radius, starting from which the function $g(r_{p})$ drops sharply (Turco Reference Turco1992) that allows us to neglect the contribution from particles of ${r_{p}}\gt r_{p,max }$ .

Although small particles of $r_{p}\lt r_{p,min }$ do not reach the surface, they induce acoustic waves, which radiate to the boundary layer and can excite the wave packets of unstable modes. Numerical modelling of Habeck et al. (Reference Habeck2024) shows that such a secondary mechanism is two orders of magnitude weaker than the primary one. In addition, numerical simulations by Chuvakhov et al. (Reference Chuvakhov, Fedorov and Obraz2019) indicate that the amplitude of boundary-layer disturbances depends weakly (within 5 % of its value) on the part of a particle trajectory being outside of the boundary layer. Therefore, the assumption of negligible contribution from particulates of $r_{p}\lt r_{p,min }$ seems reasonable.

We now discuss the amplitude criterion for transition onset. Figure 2 shows the critical RMS value of wall pressure fluctuations, $p^{\prime}_{w,cr}/P_{e}$ , as a function of the local Mach number $M_{e}$ . These data were obtained for the transition onsets measured in experiments or determined from direct numerical simulations. Note that in the notations of Marineau (Reference Marineau2017), $p^{\prime}_{w,cr}/P_{e}=A_{T}=0.6A_{max }$ . The red points were determined from the numerical solutions for flat plates and a parabolic airfoil of 10 % thickness ratio (Chuvakhov & Egorov Reference Chuvakhov and Egorov2021, Reference Chuvakhov and Egorov2023; Chuvakhov & Pogorelov Reference Chuvakhov and Pogorelov2023; Chuvakhov & Fedorov Reference Chuvakhov and Fedorov2024). In the range $1\lt M_{e}\lt 6$ , the data are well approximated by an exponential function (dashed line). In the case of $M_{e}=3$ to be considered in § 3, the transition onset corresponds to $p_{w,cr}^{^{\prime}}/P_{e}=0.6\mathrm{\, \% }$ .

Figure 2.

The critical RMS value of wall pressure fluctuations, $p^{\prime}_{w,cr}/P_{e}$ , as a function of the local Mach number, $M_{e}$ . The red points were determined using the numerical solutions for flat plates and a parabolic airfoil of 10 % thickness ratio (Chuvakhov & Egorov Reference Chuvakhov and Egorov2021, Reference Chuvakhov and Egorov2023; Chuvakhov & Pogorelov Reference Chuvakhov and Pogorelov2023; Chuvakhov & Fedorov Reference Chuvakhov and Fedorov2024).

It is convenient to perform computations in two steps. First, $N$ -factors of the dominant waves, $N_{s}(x)=N(\omega _{s},\beta _{s};x)$ , and their derivatives, $N_{\omega \omega }, N_{\beta \beta }, N_{\omega \beta }$ , are calculated and stored for a sufficiently wide range of observation stations $x$ . This is the most time-consuming part of the computations. Second, the initial spectral density $\varepsilon (\omega _{s},\beta _{s},x_{0},r_{p})$ , the integral (2.18), the distribution (2.17) and the transition onset $x_{t}$ are computed for a specific distribution $g(r_{p})$ . Since this step takes an order of magnitude less time than the first step, it is feasible to perform parametric calculations for particulates of different densities $\rho _{p}$ and distributions $g(r_{p})$ .

3. Application of the statistical model

Using (2.17)–(2.18), calculations are performed for the free stream corresponding to the standard atmosphere at an altitude of 20 km: $M_{{\infty }}=4$ , $T_{{\infty }}^{\mathrm{*}}=216.65$ K, $Re_{1,{\infty }}= 7.381\times 10^{6}$ m⁻¹. The wedge temperature is adiabatic, while the wall-temperature disturbance is zero. The gas is perfect with specific heat ratio $\gamma =1.4$ and Prandtl number $Pr =0.72$ . Calculations of inviscid flow over a wedge of half-angle $\theta =14^{\circ}$ provide the flow parameters at the upper boundary-layer edge: $M_{e}=3$ , $T_{e}^{\mathrm{*}}=324.86$ K, $U_{e}^{*}=1084$ m s⁻¹, $Re_{1,e}=11.255\times 10^{6}$ m⁻¹.

To calculate the undisturbed laminar boundary layer a self-similar solution of zero pressure gradient is used; i.e. the viscous–inviscid interaction is neglected. The dynamic viscosity is calculated using the Sutherland law with a constant of 110.4 K; the bulk viscosity is assumed to be zero according to the Stokes hypothesis.

Particle trajectories in the shock layer are calculated using the method reported by Fedorov (Reference Fedorov2013). The particle drag coefficient is calculated using the empirical correlation of Crowe (Reference Crowe1967). The distributions $g^{*}(r_{p}^{*})$ , the chemical compositions of particles and their densities, $\rho _{p}^{\mathrm{*}}$ , are given in table 1 and illustrated in figures 3 and 4. To determine the minimum radius $r_{p,min }$ in (2.18), the following criterion is used: $r_{p}\lt r_{p,min }$ , if the angle of the particle entry into the boundary layer satisfies the condition $\sin \theta _{c}\lt 10^{-3}$ .

Figure 3.

Radial distribution of particulates $g^{\mathrm{*}}(r_{p}^{\mathrm{*}})$ at an altitude of 16.75 km according to Habeck et al. (Reference Habeck, Flaten and Candler2020). Data are taken from figure 1a of Al Hasnine et al. (Reference Al Hasnine, Russo and Brehm2023a); $r_{p0}^{\mathrm{*}}=1$ µm.

Figure 4.

Radial distribution of particulates in a volcanic cloud after an eruption. Data are taken from figure 3 of Turco (Reference Turco1992).

Table 1.

Parameters of atmospheric particulates.

Data of Turco (Reference Turco1992) indicate that the distribution of $g^{*}(r_{p}^{*})$ drops sharply for $r_{p}^{\mathrm{*}}\gt 10$ µm. Moreover, large particles with $r_{p}^{\mathrm{*}}\gt 10$ µm were recorded very rarely in experiments of Habeck et al. (Reference Habeck, Flaten and Candler2020) and Lawrence et al. (Reference Lawrence2023). Therefore, the upper limit of integration in (2.18) is assumed to be equal to $r_{p,max }^{\mathrm{*}}=10$ μm ( $r_{p,max }=10^{-5}$ ).

The distributions of the $N$ -factor and its derivatives do not depend on the particle radius and other parameters of the receptivity problem. The calculations of these quantities are performed using the linear stability theory in the local-parallel approximation. Figure 5 shows the increments $\sigma (x)=-\alpha _{i}(x)$ , $N$ -factors $N(x_{0},x)$ and frequencies $f^{*}$ of the waves that are dominant at the observation stations marked by circles in figure 5(b). The envelope of the $N(x_{0},x)$ curves is shown by the thick line and corresponds to the function $N_{s}(x_{0},x)$ in (2.17). Data are shown only for the first nine observation stations relevant to the range $2\lt N_{s}\lt 14$ . Computations were performed for 15 stations relevant to $2\lt N_{s}\lt 21$ . The wave packets are dominated by oblique waves having wavefront angles $\psi =\arctan (\beta _{r}/\alpha _{r})\simeq 70^{\circ}$ and frequencies $f^{\mathrm{*}}\sim 10\text{ kHz}$ , which are typical for the first-mode instability. The corresponding neutral points, $x_{0}(\omega _{s},\beta _{s})$ , are located close to the wedge leading edge, with most of them being at $x_{0}\lt 0.1$ , where submicrometre-sized particles can penetrate into the boundary layer.

Figure 5.

(a) Linear stability characteristics of the dominant waves of the first mode for various observation stations marked by circles in (b).

We consider the functions included in the integral (2.18) and characterizing the process of wave-packet excitation near the neutral points in the ‘Minnesota’ case of normal atmospheric conditions (see table 1). In the other cases, the functions behave similarly. Figure 6 shows the functions $F(r_{p}^{\mathrm{*}})$ for the nine neutral points $x_{0}$ where $\sigma =0$ (see figure 5 a). These functions tend to zero as $r_{p}^{\mathrm{*}}\rightarrow 0$ because smaller particles decelerate more rapidly. The curves end at the stations where $r_{p}^{*}=r_{p,min }^{*}$ . For large particles $F(r_{p}^{\mathrm{*}})\rightarrow 0$ because their concentration decreases with growing $r_{p}^{\mathrm{*}}$ . As a result, maximums of the order of $F\sim 10^{7}$ are observed in the vicinity of $r_{p}^{\mathrm{*}}\simeq 1$ μm.

Figure 6.

Functions $F(r_{p}^{\mathrm{*}})$ for various observation stations marked in figure 5(b).

Figure 7 shows the functions $| \varepsilon (\omega _{s},\beta _{s},x_{0},r_{p})|$ characterizing receptivity. Large particles weakly decelerate in the shock layer and enter the boundary layer at speeds close to the free-stream speed. The drag force of such particles increases roughly as $r_{p}^{2}$ that leads to the increase of $\varepsilon (r_{p})$ . Since this trend is not completely compensated by the decrease of $F(x_{0},r_{p})$ (see figure 6), a slow increase in the integrand function

(3.1) \begin{equation}\phi \left(x_{0},r_{p}\right)=F(x_{0},r_{p})\left| \varepsilon (\omega _{s},\beta _{s},x_{0},r_{p})\right| ^{2}\end{equation}

is observed in figure 8. This behaviour leads to the increase in the integral $\Phi (\omega _{s},\beta _{s},x_{0})$ and, consequently, to the increase in $\langle p_{w}^{\prime 2}\rangle$ as $r_{p,max }$ grows. However, this trend is relatively weak as shown in figure 9 where the distributions $p^{\prime}_{w,RMS}(x)$ are plotted for $r_{p,max }^{*}=5$ and $10$ μm. Exploiting this fact and taking into account that the particle concentration drops sharply for $r_{p}^{\mathrm{*}}\gt 10$ μm, we set $r_{p,max }^{\mathrm{*}}=10$ μm in computations of the distributions $p^{\prime}_{w,RMS}(x)$ for all the cases considered.

Figure 7.

Spectral density modulus $| \varepsilon (\omega _{s},\beta _{s},x_{0},r_{p})|$ for the wave packets associated with different observation stations marked in figure 5(b).

Figure 8.

Integrands $\phi (x_{0},r_{p}^{\mathrm{*}})$ of (3.1) for various observation stations marked in figure 5(b).

Figure 9.

The distributions $p^{\prime}_{w,RMS}(x)/P_{e}$ calculated for two values of $r_{p,\max }^{*}$ .

These distributions are shown in figure 10. For the case of a volcanic cloud, two subcases are illustrated: flight in the cloud one week and one month past the eruption. To evaluate the role of particles in the transition process, the distribution of $p^{\prime}_{w,RMS}(x)/P_{e}$ for the perfectly ‘clean’ case, where transition is caused by thermal fluctuations, is shown by the green line. This distribution was calculated using the theoretical model developed by Fedorov & Averkin (Reference Fedorov, Averkin, Schlatter and Henningson2010) and Fedorov & Tumin (Reference Fedorov and Tumin2017). The black line shows the envelope $N_{s}(x_{0},x)$ . Arrows and numbers indicate the values of $N_{s}(x_{0},x_{t})$ corresponding to the critical level, $p^{\prime}_{w,cr}/P_{e}=0.6\mathrm{\, \% }$ , that was determined using the data in figure 2. As expected, the transition onset due to thermal fluctuations ( $x_{t}\approx 4.5$ ) corresponds to the largest value of $N_{s}\approx 16.8$ . In the ‘Minnesota’ case of normal atmospheric conditions (red line), the transition onset is predicted at $x_{t}\approx 3.8$ and $N_{s}\approx 15.3$ . Particles from the rocket exhaust (pink curve) shift transition to station $x_{t}\approx 3.0$ where $N_{s}\approx 13.5$ . Note that the data of Turco (Reference Turco1992) were obtained more than 23 years ago. Since that time, the frequency of space flights has increased significantly. Therefore, the present-day concentration of such particles may be many times higher, and it will continue to grow.

Figure 10.

Distributions of $p^{\prime}_{w,RMS}(x)/P_{e}$ for the cases of table 1. Green line, perfectly ‘clean’ case where transition is due to thermal fluctuations; thick black line, the $N$ -factor envelope $N_{s}(x_{0},x)$ ; thin black line, the critical level $p^{\prime}_{w,cr}/P_{e}=0.6\,\mathrm{\% }$ ; arrows and numbers indicate the transition onset $x_{t}$ and the values of $N_{s}(x_{0},x_{t})$ .

As expected, volcanic clouds cause earlier transition. The transition onset is predicted at $x_{t}\approx 2.8$ and $N_{s}\approx 12.9$ (cyan line) for flight through the cloud one month past the eruption, and $x_{t}\approx 2.5$ and $N_{s}\approx 12.2$ (blue line) one week past the eruption. The one-week cloud contains many large particles of radius $3\lt r_{p}^{\mathrm{*}}\lt 10$ μm (figure 4). They make a major contribution to the initial disturbance, as demonstrated by the function $\phi (x_{0},r_{p})$ in figure 11. In the one-month cloud, the concentration is relatively high for small particles ( $0.05\lt r_{p}^{\mathrm{*}}\lt 3$ μm) and low for large particles ( $r_{p}^{\mathrm{*}}\gt 3$ μm) (see figure 4). Nevertheless, this difference weakly affects the transition onset: $| \Delta x_{t}| /x_{t}\approx 10\, \%$ .

Figure 11.

Integrand functions $\phi (x_{0},r_{p}^{\mathrm{*}})$ for various observation stations marked in figure 5(b) (station numbers increase along the arrows). Flight though the cloud one week (blue lines) and one month (orange lines) past the eruption.