2.1 DNS

The DNS of the spatially evolving flat-plate boundary layer is performed by solving the three-dimensional compressible Navier–Stokes (N–S) equations

(2.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{U}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\boldsymbol{F}_{j}}{\unicode[STIX]{x2202}x_{j}}-\frac{\unicode[STIX]{x2202}\boldsymbol{V}_{j}}{\unicode[STIX]{x2202}x_{j}}=0,\end{eqnarray}$$

with

(2.2) $$\begin{eqnarray}\boldsymbol{U}=\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D70C}\\ \unicode[STIX]{x1D70C}u_{1}\\ \unicode[STIX]{x1D70C}u_{2}\\ \unicode[STIX]{x1D70C}u_{3}\\ E\end{array}\right),\quad \boldsymbol{F}_{j}=\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D70C}u_{j}\\ \unicode[STIX]{x1D70C}u_{1}u_{j}+p\unicode[STIX]{x1D6FF}_{1j}\\ \unicode[STIX]{x1D70C}u_{2}u_{j}+p\unicode[STIX]{x1D6FF}_{2j}\\ \unicode[STIX]{x1D70C}u_{3}u_{j}+p\unicode[STIX]{x1D6FF}_{3j}\\ (E+p)u_{j}\end{array}\right),\quad \boldsymbol{V}_{j}=\left(\begin{array}{@{}c@{}}0\\ \unicode[STIX]{x1D70E}_{1j}\\ \unicode[STIX]{x1D70E}_{2j}\\ \unicode[STIX]{x1D70E}_{3j}\\ \unicode[STIX]{x1D70E}_{jk}u_{k}+\unicode[STIX]{x1D705}\displaystyle \frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}x_{j}}\end{array}\right).\end{eqnarray}$$

Here, the Einstein summation is used, and the subscript $j=1,2,3$ denotes the index in three-dimensional Cartesian coordinates. The coordinates $x_{1}$ , $x_{2}$ and $x_{3}$ are equivalent to $x$ , $y$ and $z$ , respectively. The velocity components are denoted by $u_{j}$ in the $j$ th coordinate direction, or $u$ , $v$ and $w$ in the $x$ , $y$ and $z$ directions, respectively; $\unicode[STIX]{x1D70C}$ is the density, $p$ the static pressure, $\unicode[STIX]{x1D6FF}_{ij}$ the Kronecker delta function, $\unicode[STIX]{x1D705}$ the thermal conductivity and $T$ the temperature. The N–S equations (2.1) are non-dimensionalized with the free-stream quantities $\unicode[STIX]{x1D70C}_{\infty }$ , $T_{\infty }$ , $U_{\infty }$ and $\unicode[STIX]{x1D707}_{\infty }$ . The length and time scales are non-dimensionalized by the reference length $L=1$ inch and time $L/U_{\infty }$ , respectively.

It is noted that the subscript ‘ $\infty$ ’ denotes the free-stream properties, and ‘ $w$ ’ denotes the quantities on the wall. The free-stream Mach number is $Ma_{\infty }\equiv U_{\infty }/c_{\infty }$ where $c_{\infty }$ is the free-stream sound speed. The free-stream Reynolds number is $Re_{\infty }\equiv \unicode[STIX]{x1D70C}_{\infty }U_{\infty }L/\unicode[STIX]{x1D707}_{\infty }$ , and the wall temperature $T_{w}$ is normalized by the free-stream temperature $T_{\infty }=169.4$  K.

In the perfect gas assumption, the total energy is

(2.3) $$\begin{eqnarray}E=\unicode[STIX]{x1D70C}c_{v}T+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}u_{i}u_{i},\end{eqnarray}$$

where $c_{v}$ is the specific heat at constant volume. For a compressible Newtonian flow, the viscous stress is

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{ij}=\unicode[STIX]{x1D707}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}\right)-\frac{2}{3}\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}u_{k}}{\unicode[STIX]{x2202}x_{k}}\unicode[STIX]{x1D6FF}_{ij},\end{eqnarray}$$

where the dynamic viscosity $\unicode[STIX]{x1D707}$ is assumed to obey Sutherland’s law.

As shown in figure 2, the computational domain is bounded by an inlet boundary and a non-reflecting outlet boundary in the streamwise $x$ -direction, a wall boundary and a non-reflecting upper boundary in the wall-normal $y$ -direction and two periodic boundaries in the spanwise $z$ -direction. The inlet boundary is given by a laminar compressible boundary-layer similarity solution. The wall temperature is prescribed by an approximate adiabatic (or recovery) temperature

(2.5) $$\begin{eqnarray}T_{w}=1+r\frac{\unicode[STIX]{x1D6FE}-1}{2}Ma_{\infty }^{2},\end{eqnarray}$$

so that the wall can be considered as quasi-isothermal (see Duan, Beekman & Martin Reference Duan, Beekman and Martin2010). Here, $\unicode[STIX]{x1D6FE}=1.4$ denotes the ratio of specific heats at constant pressure and volume, and $r=Pr^{1/3}$ denotes the flat-plate turbulent recovery factor (see White Reference White2006) with the Prandtl number $Pr=0.72$ .

Figure 2.

A schematic diagram of the computational domain with boundary conditions. The inclined structure is sketched by the dashed line with the inclination angle $\unicode[STIX]{x1D6FC}$ .

In order to trigger the laminar–turbulent transition, a region of blowing and suction at $x_{a}\leqslant x\leqslant x_{b}$ is imposed. The wall-normal disturbing velocity is generated as (see Pirozzoli, Grasso & Gatski Reference Pirozzoli, Grasso and Gatski2004; Gao et al. Reference Gao, Fu, Ma and Li2005)

(2.6) $$\begin{eqnarray}v_{bs}(x,z,t)=Af(x)g(z)h(t),\quad x_{a}\leqslant x\leqslant x_{b}.\end{eqnarray}$$

Here, $A$ is the amplitude of disturbance, and the functions are

(2.7) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}f(x)=4\sin \unicode[STIX]{x1D6E9}(1-\cos \unicode[STIX]{x1D6E9})/\sqrt{27},\\ \unicode[STIX]{x1D6E9}=2\unicode[STIX]{x03C0}(x-x_{a})/(x_{a}-x_{b}),\end{array}\right\} & & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle g(z)=\mathop{\sum }_{l=1}^{l_{max}}Z_{l}\sin [2\unicode[STIX]{x03C0}l(z/L_{z}+\unicode[STIX]{x1D719}_{l})],\\ \displaystyle \mathop{\sum }_{l=1}^{l_{max}}Z_{l}=1,\quad Z_{l}=1.25Z_{l+1},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle h(t)=\mathop{\sum }_{m=1}^{m_{max}}T_{m}\sin [2\unicode[STIX]{x03C0}(m\unicode[STIX]{x1D6FD}t+\unicode[STIX]{x1D719}_{m})],\\ \displaystyle \mathop{\sum }_{m=1}^{m_{max}}T_{m}=1,\quad T_{m}=1.25T_{m+1},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $L_{z}$ is the size of the computational domain in the spanwise direction, $\unicode[STIX]{x1D6FD}$ is the base frequency of disturbance, $\unicode[STIX]{x1D719}_{l}$ and $\unicode[STIX]{x1D719}_{m}$ are random numbers ranging from 0 to 1 and $l_{max}=10$ and $m_{max}=5$ are maximum numbers of disturbing modes. In general, the transition is delayed by decreasing $A$ , and the transitional region is extended and the overshoot of $c_{f}$ at the end of transition is diminished by decreasing $\unicode[STIX]{x1D6FD}$ .

The DNS is implemented using the OpenCFD code (Li, Fu & Ma Reference Li, Fu and Ma2010), which has been widely used and validated in compressible transitional and turbulent flows (see Li et al. Reference Li, Fu and Ma2010; Zhang et al. Reference Zhang, Bi, Hussain and She2014; Zheng et al. Reference Zheng, Yang and Chen2016). The N–S equations (2.1) are integrated in time by using the third-order total-variation-diminishing-type Runge–Kutta method. The convection terms $\unicode[STIX]{x2202}\boldsymbol{F}_{j}/\unicode[STIX]{x2202}x_{j}$ are approximated by a seventh-order weighted essentially non-oscillatory scheme (Jiang & Shu Reference Jiang and Shu1996), and the viscous terms $\unicode[STIX]{x2202}\boldsymbol{V}_{j}/\unicode[STIX]{x2202}x_{j}$ are approximated by an eighth-order central finite-difference scheme.

In the present DNS, two Mach numbers are considered, and the parameters are listed in table 1. The computational domain with size $L_{x}\times L_{y}\times L_{z}$ is discretized using grid points $N_{x}\times N_{y}\times N_{z}$ with grid spacing $\unicode[STIX]{x0394}x\times \unicode[STIX]{x0394}y\times \unicode[STIX]{x0394}z$ . The superscript ‘ $+$ ’ denotes a quantity normalized by a wall unit (Robinson Reference Robinson1991), i.e. the viscous length scale $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}\equiv \unicode[STIX]{x1D707}_{w}/(\unicode[STIX]{x1D70C}_{w}u_{\unicode[STIX]{x1D70F}})$ or the friction velocity $u_{\unicode[STIX]{x1D70F}}\equiv \sqrt{\unicode[STIX]{x1D70F}_{w}/\unicode[STIX]{x1D70C}_{w}}$ with the wall shear stress $\unicode[STIX]{x1D70F}_{w}$ in the fully developed turbulent region. Moreover, the real-gas effect is not considered as in previous low-enthalpy DNS studies (e.g. Pirozzoli et al. Reference Pirozzoli, Grasso and Gatski2004; Duan, Beekman & Martin Reference Duan, Beekman and Martin2011; Zhang et al. Reference Zhang, Bi, Hussain and She2014).

Table 1.

DNS parameters.

Table 2 summarizes the Mach number, Reynolds number $Re_{\unicode[STIX]{x1D703}}\equiv Re_{\infty }\unicode[STIX]{x1D703}$ based on the momentum thickness $\unicode[STIX]{x1D703}$ and grid resolution in wall units in recent DNS of compressible transitional and turbulent boundary layers. The grid resolution of the present study is comparable to those in recent DNS studies.

Table 2.

Parameters and grid resolutions in recent DNS of compressible boundary layers.

The computational domain is partitioned into four zones in the streamwise direction to reduce the computational cost. The range and number of grid points in the streamwise direction of each computational zone are listed in table 3. A coarse grid is used in the first laminar zone, containing the region of blowing and suction from $x_{a}=0.5$ to $x_{b}=1.0$ and a small fraction of the transition region. A fine grid is used in the third turbulent zone, where the flow has reached the fully developed turbulent state. A stretched grid is applied to the second transition zone to smoothly transform the coarse grid to the fine grid. A gradually coarser grid is used in the fourth sponge zone to minimize the numerical error generated by the outlet boundary condition. In the wall-normal direction, exponentially stretched grid points

(2.10) $$\begin{eqnarray}y_{j}=L_{y}\frac{\tanh [b((j-1)/(N_{y}-1)-1)]}{\tanh b}+L_{y},\quad j=1,2,\ldots ,N_{y},\end{eqnarray}$$

with $b=2.93$ , are applied to capture small-scale structures near the wall. A uniform grid is applied in the spanwise direction.

Table 3.

The range and number of grid points in the streamwise direction for four computational zones in DNS.

The present DNS results are validated by theoretical and existing numerical results in figure 3. The mean flow statistics in the fully developed turbulent region from DNS in figure 3(a) agree well with the DNS results in Pirozzoli et al. (Reference Pirozzoli, Grasso and Gatski2004). Figure 3(b) shows that the van Driest transforms of the normalized mean velocity profiles satisfy the theoretical fits in the viscous sublayer and the logarithm law region. Figure 4(a) shows that the skin-friction coefficient $c_{f}\equiv \unicode[STIX]{x1D70F}_{w}/(\frac{1}{2}\unicode[STIX]{x1D70C}_{\infty }U_{\infty }^{2})$ along the streamwise direction from DNS agrees well with the empirical formulae of $c_{f}$ in the laminar and fully developed turbulent regions. Here, the formula for the compressible laminar skin friction is

(2.11) $$\begin{eqnarray}c_{fL}(x)=\frac{0.664}{\sqrt{Re_{x}^{\ast }}},\end{eqnarray}$$

with $Re_{x}^{\ast }=\unicode[STIX]{x1D70C}^{\ast }U_{\infty }x/\unicode[STIX]{x1D707}^{\ast }$ (see Anderson Reference Anderson2010) at the reference temperature

(2.12) $$\begin{eqnarray}T^{\ast }=1+0.032Ma_{\infty }^{2}+0.58(T_{w}-1).\end{eqnarray}$$

The formula for the compressible turbulent skin friction is

(2.13) $$\begin{eqnarray}c_{fT}(x)=\frac{0.455}{S^{2}}\left[\ln \left(\frac{0.06}{S}Re_{x}\frac{1}{\unicode[STIX]{x1D707}_{w}}\sqrt{\frac{1}{T_{w}}}\right)\right]^{-2},\end{eqnarray}$$

with $Re_{x}\equiv Re_{\infty }x$ , $S=\sqrt{T_{w}-1}/\sin ^{-1}{\mathcal{A}}$ and ${\mathcal{A}}=(r((\unicode[STIX]{x1D6FE}-1)/2)Ma_{\infty }^{2}(1/T_{w}))^{1/2}$ , which is accurate over the practical range of $Re_{\unicode[STIX]{x1D703}}$ , $Ma_{\infty }$ and $T_{w}$ (see White Reference White2006).

Figure 3.

Validation of the present DNS result. (a) Mean profiles in the fully developed turbulent region at $Ma_{\infty }=2.25$ in the present DNS (lines) and in the DNS (symbols) by Pirozzoli et al. (Reference Pirozzoli, Grasso and Gatski2004). The overbar denotes the Reynolds average over the spanwise direction. (b) Normalized mean velocity profiles $\bar{u}^{+}$ (symbols) and their van Driest transforms $\bar{u}_{vd}=\int _{0}^{\bar{u}^{+}}(\bar{\unicode[STIX]{x1D70C}}/\bar{\unicode[STIX]{x1D70C}}_{w})^{1/2}\,\text{d}\bar{u}^{+}$ (solid and dashed lines) with theoretical fits $\bar{u}_{vd}=y^{+}$ in the viscous sublayer and $\bar{u}_{vd}=\ln y^{+}/0.41+5.5$ in the logarithm law region.

Figure 4.

The skin-friction coefficient along the streamwise direction. Note the spatial coordinate is normalized by the reference length $L=1$ inch in DNS. (a) Comparison of DNS results and empirical formulae (2.11) and (2.13) at two Mach numbers. (b) Convergence of $c_{f}$ with three grid resolutions listed in table 1 at $Ma_{\infty }=6$ .

The rapid increase of $c_{f}$ signals the transition, and this appears to coincide with the generation of three-dimensional vortices and the breakdown of large-scale coherent structures. The transitional regions are roughly at $2\leqslant x\leqslant 4$ and $3\leqslant x\leqslant 5$ for $Ma_{\infty }=2.25$ and $Ma_{\infty }=6$ , respectively.

The grid convergence of the skin friction for $Ma_{\infty }=6$ is assessed using three different grid resolutions in figure 4(b). Since $c_{f}$ is converged for the moderate and fine grids, the moderate grids are used in the following study. Furthermore, we find that the tails of the two-point correlations of the velocity components in the spanwise direction are sufficiently small at lateral boundaries (not shown), which ensures that $L_{z}$ is large enough so that the imposed periodic boundary condition does not affect flow statistics (Pirozzoli et al. Reference Pirozzoli, Grasso and Gatski2004).

2.2 Lagrangian scalar field

The three-dimensional Lagrangian scalar field $\unicode[STIX]{x1D719}(\boldsymbol{x},t)$ is an ideal passive tracer field, and isosurfaces of $\unicode[STIX]{x1D719}$ are material surfaces in the temporal evolution. The scalar field is governed by the pure convection equation

(2.14) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=0.\end{eqnarray}$$

Although $\unicode[STIX]{x1D719}$ is similar to the temperature or species concentration, it has no diffusivity so that it has no smallest scale in the evolution at long times. Thus it can display a rich geometry of flow structures at very fine scales, and excludes the subjective selection of the diffusivity value.

The backward-particle-tracking method (Yang et al. Reference Yang, Pullin and Bermejo-Moreno2010; Yang & Pullin Reference Yang and Pullin2011), which is numerically stable and has no numerical dissipation, is used to solve (2.14). In this method, trajectories of fluid particles are calculated by solving the kinematic equation

(2.15) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{X}(\boldsymbol{x}_{0},t_{0}|t)}{\unicode[STIX]{x2202}t}=\boldsymbol{V}(\boldsymbol{x}_{0},t_{0}|t)=\boldsymbol{u}(\boldsymbol{X}(\boldsymbol{x}_{0},t_{0}|t),t),\end{eqnarray}$$

where $\boldsymbol{X}(\boldsymbol{x}_{0},t_{0}|t)$ is the location at time $t$ of a fluid particle which was located at $\boldsymbol{x}_{0}$ at the initial time $t_{0}$ , $\boldsymbol{V}(\boldsymbol{x}_{0},t_{0}|t)$ is the Lagrangian velocity of the fluid particle and $\boldsymbol{u}(\boldsymbol{X}(\boldsymbol{x}_{0},t_{0}|t),t)$ is its local Eulerian velocity.

The initial Lagrangian field is uniquely determined as $\unicode[STIX]{x1D719}(\boldsymbol{x}_{0},t_{0})=y$ based on the criteria of the best approximation of vortex sheets and geometric invariance of $\unicode[STIX]{x1D719}$ in the laminar state (Zhao et al. Reference Zhao, Yang and Chen2016). The isosurfaces of $\unicode[STIX]{x1D719}(\boldsymbol{x}_{0},t_{0})$ are streamwise–spanwise planes at different distances from the wall at the initial time. Near the inlet, the amplitudes of initial disturbances are very small compared to the mean flow, so that the material sheets are almost invariant. As the disturbances develop with increasing $x$ , the material surfaces are stretched and folded in the streamwise and wall-normal plane. Subsequently, some parts of the material surfaces can be rolled up into complex shapes with three-dimensional geometric characteristics in the evolution (Zhao et al. Reference Zhao, Yang and Chen2016).

In the numerical implementation of the backward-particle-tracking method, $\unicode[STIX]{x1D719}$ at a given time $t$ is calculated as follows (Zheng et al. Reference Zheng, Yang and Chen2016).

(i) The full Eulerian velocity field on the grid $N_{x}\times N_{y}\times N_{z}$ within a time period from $t_{0}$ to $t>t_{0}$ is solved by DNS and then stored to disk.

(ii) As sketched in figure 2, at the given time $t$ , particles are placed at uniform grid points of $N_{x}^{p}\times N_{y}^{p}$ in a subdomain of interest in a number of $x$ – $y$ planes for further geometric analysis. The time period $\unicode[STIX]{x0394}T\equiv t-t_{0}$ for the backward-particle tracking in the implementation is set to be $\unicode[STIX]{x0394}T=1$ following the selection argument in Zheng et al. (Reference Zheng, Yang and Chen2016). Numerical experiments demonstrate that the statistical results for $\unicode[STIX]{x1D719}$ are not sensitive to the value of $\unicode[STIX]{x0394}T\geqslant 1$ .

(iii) The particles are released and their trajectories are calculated backward in time within $\unicode[STIX]{x0394}T$ or until they arrive at $x_{0}=0.2$ where the initial material surfaces are considered as planar sheets. A three-dimensional, fourth-order Lagrangian interpolation scheme is used to obtain the fluid velocity at the location of each particle. An explicit, second-order Adams–Bashforth scheme is applied for the time integration, and time steps $\unicode[STIX]{x0394}t=0.0025$ for $Ma_{\infty }=2.25$ and $\unicode[STIX]{x0394}t=0.005$ for $Ma_{\infty }=6$ are selected by the criterion $\unicode[STIX]{x0394}t\leqslant \unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}/u_{\unicode[STIX]{x1D70F}}$ for capturing the finest resolved scales in the velocity field.

(iv) After the backward tracking, we obtain initial locations $\boldsymbol{x}_{0}$ of the particles and the flow map

(2.16) $$\begin{eqnarray}F_{t}^{t_{0}}(\boldsymbol{X})\boldsymbol{ : }\boldsymbol{X}(\boldsymbol{x}_{0},t_{0}|t)\mapsto \boldsymbol{x}_{0},\quad t\geqslant t_{0}.\end{eqnarray}$$

Then the Lagrangian field $\unicode[STIX]{x1D719}(\boldsymbol{x},t)$ at any given time $t$ is obtained as

(2.17) $$\begin{eqnarray}\unicode[STIX]{x1D719}(\boldsymbol{x},t)=\unicode[STIX]{x1D719}(F_{t}^{t_{0}}(\boldsymbol{X}),t_{0}).\end{eqnarray}$$

2.3 Geometric diagnostic methodology

A systematic diagnostic tool has been developed to quantify the multi-scale geometry of flow structures in the transition of boundary layers (see Zheng et al. Reference Zheng, Yang and Chen2016). It is also a general signal processing method for characterizing preferential orientations in the scalar pattern of images from DNS (see § 3) and experiments (see appendix A).

First, a sliding window filter is used to extract the spatially evolving scalar structures and capture their major deformation in the compact frames at different streamwise locations. The filter is defined as

(2.18) $$\begin{eqnarray}f(x)=\exp \left[-n\left(\frac{x-x_{c}}{l_{w}}\right)^{n}\right],\end{eqnarray}$$

where $x_{c}$ is the position of the window centre, $l_{w}$ is a characteristic width of the window and $n$ is a positive even integer. Figure 5 illustrates the profile of the filter. From (2.18), the window width $l\equiv 2(x-x_{c})$ of the filter for a given $f$ is

(2.19) $$\begin{eqnarray}l=2\,l_{w}\exp \left[\frac{\ln (-\ln f/n)}{n}\right].\end{eqnarray}$$

The two-dimensional scalar field for further analysis is extracted in a compact region with the cutoff window width $l=l_{c}$ given $f=0.0001$ at the boundaries. The major features of the extracted scalar field with a large scalar gradient $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ are captured in the subdomain with the effective window width $l=l_{e}$ for $f\geqslant 0.9$ . As shown in figure 5, we choose $n=16$ to smooth the transition of $f$ from $f=0$ to $1$ with $l_{c}/l_{w}=1.93$ and $l_{e}/l_{w}=1.46$ .

Figure 5.

The profile of the sliding window function for extracting $\unicode[STIX]{x1D719}_{f}$ from the entire $\unicode[STIX]{x1D719}$ .

Setting $l_{e}=1$ and given $x_{c}$ , the extracted Lagrangian scalar field $\unicode[STIX]{x1D719}_{f}(\boldsymbol{x},t)$ is obtained as

(2.20) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{f}(\boldsymbol{x},t)=\unicode[STIX]{x1D719}(\boldsymbol{x},t)f+\unicode[STIX]{x1D719}(\boldsymbol{x}_{0},t_{0})(1-f).\end{eqnarray}$$

The extracted Lagrangian field in the sliding window has a smooth transition from an evolving scalar $\unicode[STIX]{x1D719}(\boldsymbol{x},t)$ to the initial scalar $\unicode[STIX]{x1D719}(\boldsymbol{x}_{0},t_{0})$ at lateral boundaries of the window, and the effect of this artificial transition on the further directional analysis is negligible. Thus the periodic boundary condition of $\unicode[STIX]{x1D719}_{f}$ in the $x$ -direction and the mirror extension of $\unicode[STIX]{x1D719}_{f}$ in the $y$ -direction can facilitate the fast Fourier transform of $\unicode[STIX]{x1D719}_{f}$ in further analysis without the numerical artefact due to possible discontinuities at the boundaries (Zheng et al. Reference Zheng, Yang and Chen2016).

The multi-scale and multi-directional geometric analysis (Yang & Pullin Reference Yang and Pullin2011) is applied to characterize the evolutionary geometry of $\unicode[STIX]{x1D719}_{f}$ on a series of streamwise and wall-normal ( $x$ – $y$ ) planes. This provides quantitative statistics on the inclination angle of scalar structures in the images at different scales and streamwise locations during the transition.

The mechanism of the multi-scale and multi-directional geometric analysis is briefly explained here, and more details and illustrative examples can be found in Yang & Pullin (Reference Yang and Pullin2011). As sketched in figure 6(a), if a two-dimensional scalar field or image $\unicode[STIX]{x1D711}(\boldsymbol{x})\in \mathbb{R}^{2}$ has a characteristic scale and a preferential orientation, then the geometric features are quantified from its Fourier transform $\hat{\unicode[STIX]{x1D711}}(\boldsymbol{k})={\mathcal{F}}\{\unicode[STIX]{x1D711}(\boldsymbol{x})\}$ as follows.

A filtered $\unicode[STIX]{x1D711}(\boldsymbol{x})$ at scale $j$ and along direction $l$ can be obtained by the inverse Fourier transform of $\hat{\unicode[STIX]{x1D711}}(\boldsymbol{k})$ filtered by the frequency window function $U_{j}(r,\unicode[STIX]{x1D703})=2^{-3j/4}W(2^{-j}r)V(t_{l}(\unicode[STIX]{x1D703}))$ in polar coordinates $r=\sqrt{k_{1}^{2}+k_{2}^{2}}$ and $\unicode[STIX]{x1D703}=\arctan (k_{2}/k_{1})$ with $j\in \mathbb{N}$ . Here, the radial window function $W(r)$ and the angular window function $V(t_{l})$ are defined by explicit functions (Yang & Pullin Reference Yang and Pullin2011).

As illustrated in figure 6(b), each frequency window function in Fourier space is supported on the region of a circular wedge in the range of wavenumbers $2^{j-1}\leqslant r\leqslant 2^{j+1}$ with $j=0,1,\ldots ,j_{e}$ . This implies that the corresponding spatial structure in physical space is a needle-shaped element or a curvelet (Candes et al. Reference Candes, Demanet, Donoho and Ying2006) with the characteristic length $2^{-j/2}$ and width $2^{-j}$ . Here, $j_{e}=\log _{2}[\min (k_{1,\max },k_{2,\max })/2]$ denotes the maximum number of scales, where $k_{1,\max }$ and $k_{2,\max }$ denote maximum resolved wavenumbers in the discrete Fourier transform of $\unicode[STIX]{x1D711}(\boldsymbol{x})$ .

Figure 6.

Schematic plots of (a) the inclined structure with characteristic angle and scale in physical space and (b) frequency window functions $U_{j}(r,\unicode[STIX]{x1D703})$ supported on circular wedges in Fourier space with the region highlighted by the red circle for the characteristic angle and scale.

The multi-scale decomposition of $\unicode[STIX]{x1D711}(\boldsymbol{x})$ gives the filtered component field $\unicode[STIX]{x1D711}_{j}(\boldsymbol{x})={\mathcal{F}}^{-1}\{\hat{\unicode[STIX]{x1D711}}(\boldsymbol{k})W(2^{-j}r)\}$ with a characteristic length scale ${\mathcal{L}}_{j}=2^{-j}$ from the inverse Fourier transform, where the band-pass filter $W(r)$ is compact in Fourier space. In the boundary-layer transition, the normalized characteristic length scales of structures in the multi-scale decomposition are listed in table 4. Subsequently, ${\mathcal{L}}_{j}\geqslant 0.5\unicode[STIX]{x1D6FF}$ will be referred to as ‘large scale’, ${\mathcal{L}}_{j}\leqslant 50\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ as ‘small scale’ and in between as ‘intermediate scale’.

Table 4.

Characteristic length scales ${\mathcal{L}}_{j}$ normalized by the boundary-layer thickness $\unicode[STIX]{x1D6FF}$ and viscous length scale $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ at $x=7$ .

The orientation information of each filtered component field of $\unicode[STIX]{x1D711}(\boldsymbol{x})$ is characterized by the normalized angular spectrum

(2.21) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{j}(\unicode[STIX]{x0394}\unicode[STIX]{x1D703})=\frac{\displaystyle \int \hat{\unicode[STIX]{x1D711}}(\boldsymbol{k})U_{j}(r,\unicode[STIX]{x1D703})\,\text{d}\boldsymbol{k}}{\displaystyle \int U_{j}(r,\unicode[STIX]{x1D703})\,\text{d}\boldsymbol{k}},\quad -\frac{\unicode[STIX]{x03C0}}{2}\leqslant \unicode[STIX]{x0394}\unicode[STIX]{x1D703}\leqslant \frac{\unicode[STIX]{x03C0}}{2},\end{eqnarray}$$

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}l^{\prime }2^{-\lceil j/2\rceil }/2$ , $-2^{\lceil j/2\rceil }\leqslant l^{\prime }\leqslant 2^{\lceil j/2\rceil }$ is the discrete deviation angle away from the $x$ -axis in the physical space for scale $j$ . As illustrated in figure 2, we define the inclination angle $\unicode[STIX]{x1D6FC}$ between an inclined structure projected onto the streamwise and wall-normal ( $x$ – $y$ ) plane and the $x$ -direction. The averaged inclination angle away from $x$ -axis at each scale $j$ is calculated by

(2.22) $$\begin{eqnarray}\langle \unicode[STIX]{x1D6FC}({\mathcal{L}}_{j})\rangle =\frac{\displaystyle \mathop{\sum }_{l^{\prime }=0}^{l_{\max }^{\prime }}\unicode[STIX]{x1D6F7}_{j}(\unicode[STIX]{x0394}\unicode[STIX]{x1D703})\unicode[STIX]{x0394}\unicode[STIX]{x1D703}}{\displaystyle \mathop{\sum }_{l^{\prime }=0}^{l_{\max }^{\prime }}\unicode[STIX]{x1D6F7}_{j}(\unicode[STIX]{x0394}\unicode[STIX]{x1D703})},\end{eqnarray}$$

with $l_{max}^{\prime }=2^{\lceil j/2\rceil }$ .

4.1 Model equation

It is inspiring that the rapid growth of the averaged inclination angle in scalar images coincides with the increasing trend of the skin-friction coefficient in the boundary-layer transition. Next, we seek a quantitative connection between $\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle$ and $c_{f}(x)$ .

From the inclination-angle statistics of the filtered component fields on the $x$ – $y$ plane at different scales, we propose a simple image-based model of the skin-friction coefficient by blending $c_{fL}$ in (2.11) and $c_{fT}$ in (2.13) as

(4.1) $$\begin{eqnarray}c_{f,model}(x)=(1-\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D6FC}})\,c_{fL}(x_{L})+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D6FC}}\,c_{fT}(x_{T}),\end{eqnarray}$$

with

(4.2) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D6FC}}(x)=\left(\frac{\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle -\overline{\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle |L_{x1}\leqslant x\leqslant x_{L}}}{\overline{\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle |x_{T}\leqslant x\leqslant L_{x2}}-\overline{\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle |L_{x1}\leqslant x\leqslant x_{L}}}\right)^{\unicode[STIX]{x1D709}}.\end{eqnarray}$$

Here the subscript ‘ $L$ ’ denotes the quantities in the laminar region, and ‘ $T$ ’ the quantities in the fully developed turbulent region; $x_{L}$ and $x_{T}$ denote the reference locations of laminar and fully developed turbulent regions, respectively; $\overline{\langle \unicode[STIX]{x1D6FC}\rangle |L_{x1}\leqslant x\leqslant x_{L}}$ and $\overline{\langle \unicode[STIX]{x1D6FC}\rangle |x_{T}\leqslant x\leqslant L_{x2}}$ denote conditional averages of $\langle \unicode[STIX]{x1D6FC}\rangle$ over laminar and turbulent regions, respectively, for normalizing $\langle \unicode[STIX]{x1D6FC}\rangle$ in (4.2).

Figure 9 shows that the averaged inclination angle is close to $0^{\circ }$ at $x=x_{L}$ in the laminar state, then it increases sharply in the transitional region and finally it converges to a statistically steady value around $x=x_{T}$ in the fully developed turbulent state. The reference locations $x_{L}$ and $x_{T}$ are searched from the rightmost and leftmost of profiles of $\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle$ by criteria

(4.3) $$\begin{eqnarray}\langle \unicode[STIX]{x1D6FC}(x_{L},{\mathcal{L}}_{j})\rangle =(\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle _{max}-\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle _{min})\times 1\,\%+\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle _{min},\end{eqnarray}$$

and

(4.4) $$\begin{eqnarray}\langle \unicode[STIX]{x1D6FC}(x_{T},{\mathcal{L}}_{j})\rangle =(\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle _{max}-\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle _{min})\times 99\,\%+\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle _{min},\end{eqnarray}$$

respectively. Here, $\langle \unicode[STIX]{x1D6FC}\rangle _{min}$ and $\langle \unicode[STIX]{x1D6FC}\rangle _{max}$ are the minimum and maximum of $\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle$ , respectively, in the subdomain of interest for $\unicode[STIX]{x1D719}$ .

The modelled $c_{f}$ is bounded by the empirical formulae (2.11) for $c_{fL}$ and (2.13) for $c_{fT}$ in (4.1), and the transition from $c_{fL}$ to $c_{fT}$ depends on $\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D6FC}}$ generally ranging from 0 to 1. The imaged-based model is reminiscent of the intermittency-based model of Dhawan & Narasimha (Reference Dhawan and Narasimha1958), and $\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D6FC}}$ in terms of the normalized averaged inclination angle in (4.1) can be considered as an ‘intermittency factor’ to characterize the emergence of coherent structures or the ‘turbulent spot’. The exponent in (4.2) is generally related to the steepness of the transition from $c_{fL}$ to $c_{fT}$ , and it is set to be $\unicode[STIX]{x1D709}=1$ unless a variable $\unicode[STIX]{x1D709}$ is explicitly specified. In appendix B, we demonstrate that the optimal $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}({\mathcal{L}}_{j})$ is scale dependent.

4.2 Model validation

The model (4.1) of $c_{f}$ based on averaged inclination angles at different length scales is assessed by the relative error

(4.5) $$\begin{eqnarray}\unicode[STIX]{x1D700}_{r}=\frac{\displaystyle \int _{L}|c_{f,model}-c_{f,DNS}|\,\text{d}x/L}{\displaystyle \int _{L}c_{f,DNS}\,\text{d}x/L}\times 100\,\%,\end{eqnarray}$$

with the integral length $L=L_{x2}-L_{x1}$ . The relative errors for DNS at two $Ma_{\infty }$ are shown in figure 11. In general, the discrepancy between the model prediction $c_{f,model}$ and the DNS result $c_{f,DNS}$ decreases with decreasing the characteristic length scale in the scalar image; we achieve overall good predictions of $c_{f}$ with $\unicode[STIX]{x1D700}_{r}<15\,\%$ at intermediate and small scales as $0.1\unicode[STIX]{x1D6FF}$ or ${\mathcal{L}}_{j}<100\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ . Therefore, we suggest using $c_{f,model}$ at ${\mathcal{L}}_{j}\leqslant O(0.1\unicode[STIX]{x1D6FF})$ for the model prediction from experimental images of boundary-layer transition (see appendix A).

As discussed in appendix B, the model prediction can be further improved by using scale-dependent $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}({\mathcal{L}}_{j})$ instead of $\unicode[STIX]{x1D709}=1$ in (4.2). Moreover, the intermittency factor in (4.1) and (4.2) can be modelled by the form suggested by Dhawan & Narasimha (Reference Dhawan and Narasimha1958), which is discussed in appendix C.

Figure 11.

The relative error of $c_{f,model}$ calculated from the model (4.1) based on $\langle \unicode[STIX]{x1D6FC}(x,{\mathcal{L}}_{j})\rangle$ at various scales in terms of (a)  ${\mathcal{L}}_{j}/\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ and (b)  ${\mathcal{L}}_{j}/\unicode[STIX]{x1D6FF}$ .

Based on the error analysis, for example, we use the values of $\langle \unicode[STIX]{x1D6FC}\rangle$ at the length scale ${\mathcal{L}}_{j}\approx 0.08\unicode[STIX]{x1D6FF}$ in the model equation (4.1). As shown in figure 12, the trend of the increasing averaged inclination angle is qualitatively similar to the rise of $c_{f}$ , and the modelled $c_{f}$ in equation (4.1) agrees well with DNS results for both Mach numbers. This implies that it is possible to quantify $c_{f}$ , which is hard to measure in high- $Ma$ boundary-layer transition, from two-dimensional scalar images that are relatively easy to obtain from experimental visualization. In addition, the slight overshoot of $c_{f}$ at the end of the transitional region in DNS is not observed in most of the profiles of $\langle \unicode[STIX]{x1D6FC}\rangle$ in figure 9, resulting in the under-prediction of $c_{f}$ from (4.1).

Figure 12.

Comparison of the skin-friction coefficients from DNS and the model (4.1) with ${\mathcal{L}}_{j}\approx 0.08\unicode[STIX]{x1D6FF}$ or $O(10)\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ , along with the averaged inclination angle of scalar structures. (a)  $Ma_{\infty }=2.25$ , (b)  $Ma_{\infty }=6$ .

For the further application of the image-based model in experiments, the scalar field for experimental visualization generally has a finite diffusivity instead of the vanishing diffusivity in the Lagrangian field in the present DNS, so the real scalar field has a finite smallest scale $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D719}}$ . Since the effect of the diffusivity on the scalar field is similar to filtering, the discrepancy assessment of the Lagrangian field in figure 11 can be still useful for ${\mathcal{L}}_{j}>\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D719}}$ if the image-based model is applied to other scalar images from experiments. In appendix A, we demonstrate that the image-based modelling of $c_{f}$ is applicable to experimental images in compressible boundary-layer transition.

4.3 Inclined structures and drag production

The relatively accurate prediction of $c_{f}$ in the image-based model indicates that the generation of inclined small-scale flow structures is closely related to the drag production. As sketched in figure 13, one possible reason is that the lifts of material surfaces during the transition, which are good surrogates of vortex surfaces consisting of vortex lines, can generate strong inclined shear layers (see Zhao et al. Reference Zhao, Yang and Chen2016, Reference Zhao, Xiong, Yang and Chen2018) to increase $c_{f}$ .

Figure 13.

Diagram of the geometry of material surfaces and typical vortex lines near the surfaces, along with the rise of $c_{f}$ . Solid lines denote vortex lines, and solid vectors $\boldsymbol{n}_{\unicode[STIX]{x1D719}}$ denote the normal of material surfaces.

Assume the flow field is filled with wall-parallel material surfaces in the laminar state with all the surface normal $\boldsymbol{n}_{\unicode[STIX]{x1D719}}\equiv \unicode[STIX]{x1D735}\unicode[STIX]{x1D719}/|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|$ pointing to the wall-normal direction. In the transitional region, the near-wall material surfaces are lifted due to the growing streamwise vorticity. This elevation event is quantified by the wall-normal Lagrangian displacement

(4.6) $$\begin{eqnarray}\unicode[STIX]{x0394}Y(\boldsymbol{x}_{0},t_{0}|t)=Y(\boldsymbol{x}_{0},t_{0}|t)-y(\boldsymbol{x}_{0},t_{0}),\end{eqnarray}$$

where $Y$ is the wall-normal location of a fluid particle on a material surface. The displacement $\unicode[STIX]{x0394}Y$ quantifies the scalar transport in the wall-normal direction within a time interval of interest, and Zhao et al. (Reference Zhao, Yang and Chen2016) define the Lagrangian events ‘elevation’ with $\unicode[STIX]{x0394}Y>0$ and ‘descent’ with $\unicode[STIX]{x0394}Y<0$ . The contour of $\unicode[STIX]{x0394}Y$ and the contour lines of high shear $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y$ in the transitional region are shown in figure 14. In general, the inclined high shear layers cover the region with $\unicode[STIX]{x0394}Y>0$ , which is similar to the observation in an incompressible temporal transitional channel flow in Zhao et al. (Reference Zhao, Yang and Chen2016), because the strong shear layer can be generated between the elevated low-speed fluid and the surrounding high-speed fluid. Furthermore, the region with $\unicode[STIX]{x0394}Y>0$ , which also corresponds to the inclined scalar structure with $\boldsymbol{n}_{\unicode[STIX]{x1D719}}$ , deviates from the wall-normal direction as sketched in figure 13, which can be characterized as a finite $\langle \unicode[STIX]{x1D6FC}\rangle$ in the multi-directional analysis. Thus the inclined high shear layers accelerate the momentum transport and produce the large Reynolds shear stress (see Zhao et al. Reference Zhao, Yang and Chen2016), which can increase $c_{f}$ implied by the relation between the Reynolds shear stress and $c_{f}$ (see Fukagata et al. Reference Fukagata, Iwamoto and Kasagi2002; Gomez, Flutet & Sagaut Reference Gomez, Flutet and Sagaut2009).

Figure 14.

The contour of the Lagrangian wall-normal displacement $\unicode[STIX]{x0394}Y$ and contour lines for the strong shear layers on the $x$ – $y$ plane in the transitional region at $Ma_{\infty }=6$ .