3.1. Initial conditions

The computational domain is initialised using the self-similar boundary-layer profiles based on Lees–Dorodnitsyn variables (Ren et al. Reference Ren, Marxen and Pecnik2019; Boldini et al. Reference Boldini, Bugeat, Peeters, Kloker and Pecnik2024b ). The initial flow profiles for all cases listed in table 1 are plotted over the wall-normal coordinate $\mathrm{d}\eta =\rho ^* u^*_\infty /\sqrt {2\xi } \,\mathrm{d}y^*$ in figure 3. As temperature increases from liquid-like to vapour-like conditions (figure 3 a), the largest gradients in thermodynamic and transport properties, particularly in density (figure 3 c), occur at the pseudo-critical point (green dashed line). This gives rise to an inflectional base-flow profile in case Tw110, as defined by the GIP criterion, i.e. $\mathrm{d}(\bar {\rho } \ \mathrm{d}\bar {u}/\mathrm{d}y)/\mathrm{d}y=0$ . As shown in figure 3(d), the GIP coincides with the minimum of kinematic viscosity near the pseudo-critical point, and occurs for non-polar fluids at supercritical pressure and transcritical temperature (Bugeat et al. Reference Bugeat, Boldini and Pecnik2022, Reference Bugeat, Boldini, Hasan and Pecnik2024). The streamwise velocity (figure 3 b) highlights the more pronounced profile of case Tw110, which resembles that of a transitional boundary-layer profile. Appendix B presents a comparison between the unperturbed 2-D DNS solutions, which have reached a steady-state solution, and the initial self-similar solutions for cases Tw095 and Tw110.

Figure 3. Laminar profiles for the considered cases: (a) temperature $T^*/T^*_\infty$ , (b) streamwise velocity $u^*/u^*_\infty$ , (c) density $\rho ^*/\rho ^*_\infty$ , and (d) kinematic viscosity $\nu ^*/\nu ^*_\infty$ , as functions of the self-similar wall-normal coordinate $\eta$ . The line legend is in agreement with table 1 for cases TadIG, Tw095 and Tw110. The dashed green line indicates the pseudo-critical point, i.e. at the pseudo-critical temperature $T^*=T^*_{pc}$ . The location of the GIP for the transcritical case Tw110 is marked by the circle symbols in (b–d).

3.2. Disturbance strip

Disturbances are introduced via a blowing and suction disturbance strip on the flat-plate surface. Once the laminar base flow reaches a steady-state solution, the localised 3-D disturbance (Sayadi et al. Reference Sayadi, Hamman and Moin2013) is activated as

(3.1) \begin{equation} v(x,y=0,z,t)=f(x) \big [ A_{{2\hbox{-}{\rm D}}} \sin (\omega _{{2\hbox{-}{\rm D}}}t) + A_{{3\hbox{-}{\rm D}}} \sin (\omega _{{3\hbox{-}{\rm D}}} t) \cos ( \beta _0 z) \big ], \end{equation}

where $A_{{2\hbox{-}{\rm D}}}=A^*_{{2\hbox{-}{\rm D}}}/u^*_\infty$ and $A_{{3\hbox{-}{\rm D}}}=A^*_{{3\hbox{-}{\rm D}}}/u^*_\infty$ are the wave amplitudes for the primary 2-D and $z$ -symmetric 3-D waves, respectively, and $\omega _{{2\hbox{-}{\rm D}}}=\omega ^*_{{2\hbox{-}{\rm D}}}\delta ^*_{99,0}/u^*_\infty$ and $\omega _{{3\hbox{-}{\rm D}}}=\omega ^*_{{3\hbox{-}{\rm D}}}\delta ^*_{99,0}/u^*_\infty$ are the corresponding angular frequencies. The spanwise wavenumber $\beta _0=\beta ^*_0\delta ^*_{99,0}=2\pi /\lambda _{z}$ equals the spanwise size $z_{e}$ of the computational domain (see Table 3). The streamwise variation in (3.1) is governed by the function $f(x)=15.1875\xi ^5-35.4375\xi ^4+20.25\xi ^3$ , with $\xi =(x-x_1)/(x_{\textit{mid}}-x_1)$ for $x_1\lt x\lt x_{\textit{mid}}$ , and $\xi =(x_2-x)/(x_2-x_{\textit{mid}})$ for $x_{\textit{mid}}\lt x\lt x_2$ , where $x_{\textit{mid}}=(x_1+x_2)/2$ corresponds to $\textit{Re}_{x,\textit{mid}}$ (see table 2). The disturbance strip is positioned upstream of branch I (see figure 4). Since no experimental studies on controlled transition with supercritical fluids are currently available, the ideal-gas H-type breakdown scenario from Sayadi et al. (Reference Sayadi, Hamman and Moin2013), validated by Boldini et al. (Reference Boldini, Hirai, Costa, Peeters and Pecnik2025), is used as a reference, with $F_{{2\hbox{-}{\rm D}}}=124 \times 10^{-6}$ and $\lambda _{z}=9.63$ . In §§ 4 and 5, 2-D and 3-D simulations are performed, respectively. In the 2-D simulations, a 2-D wave ( $\beta _0=0$ ) with frequency $F_{{2\hbox{-}{\rm D}}}$ is forced, with amplitude $A_{{2\hbox{-}{\rm D}}}=1.0 \times 10^{-8}$ in the linear regime, and increased to either $7.5 \times 10^{-4}$ or $7.5 \times 10^{-3}$ in the nonlinear regime. The latter corresponds to the amplitude used in the ideal-gas H-type breakdown reference simulation from Boldini et al. (Reference Boldini, Hirai, Costa, Peeters and Pecnik2025). For the 3-D simulations, $A_{{2\hbox{-}{\rm D}}}$ is set to $7.5 \times 10^{-3}$ , while two amplitude levels for the $z$ -symmetric 3-D wave with $F_{{3\hbox{-}{\rm D}}}=62 \times 10^{-6}$ are considered: infinitesimally small (denoted as ‘IA’) or finite (denoted as ‘LA’). The differences in forcing parameters for all simulations are summarised in table 2. The non-dimensional angular frequency $\omega$ and the frequency parameter $F$ are related by $\omega =F\, Re_{0}$ , with $F=\omega ^*\mu ^*_\infty /(\rho ^*_\infty u^{*2}_\infty )$ and local Reynolds number $\textit{Re}_{0}=\sqrt {Re_{x,0}}$ .

Table 2. Forcing set-up: ‘LA’ and ‘IA’ denote finite 3-D amplitude forcing and infinitesimally small 3-D amplitude forcing, respectively. Others parameters are fixed: $A_{{2\hbox{-}{\rm D}}}=7.5 \times 10^{-3}$ at $F_{{2\hbox{-}{\rm D}}}=124 \times 10^{-6}$ , $z$ -symmetric 3-D wave at $F_{{3\hbox{-}{\rm D}}}=62 \times 10^{-6}$ , with $\textit{Re}_{x,\textit{mid}}=1.72 \times 10^5$ (cases TadIG and Tw095) or $\textit{Re}_{x,\textit{mid}}=9.61 \times 10^4$ (case Tw110).

Figure 4. Growth-rate ( $-\alpha _{i}$ ) contours in the $\textit{Re}$ – $F$ stability diagram: (a) TadIG, (b) Tw095, and (c) Tw110 (Modes I and II). The dotted blue lines in (b,c) represent the ideal-gas neutral stability at equal $T^*_{w}/T^*_\infty$ ratios. In the inset of (c), the wide frequency band of Modes I and II is displayed. The locations of the DNS domain and perturbation strip for subharmonic breakdown, i.e. $F_{{3\hbox{-}{\rm D}}}=0.5F_{{2\hbox{-}{\rm D}}}=62\times 10^{-6}$ , are marked by white and cyan bars, respectively, as described in § 3.2.

4.1. Linear evolution

Starting from the base-flow profiles in § 3.1, linear stability analysis is performed for the cases listed in table 1. Figure 4 displays the growth rate $-\alpha _{i}$ in the $\textit{Re}$ – $F$ stability diagram, where $\textit{Re}=\sqrt {Re_{x}}$ . In figure 4(b), increasing the wall temperature towards $T_{r,pc}$ stabilises the flow. This trend is reversed under ideal-gas conditions, where a higher $T^*_{w}/T^*_\infty$ ratio leads to destabilisation (dotted blue line). In figure 4(c), upon crossing the pseudo-critical point from the liquid-like free stream, Mode II appears (Ren et al. Reference Ren, Marxen and Pecnik2019), becoming highly unstable even at low Reynolds numbers. Simultaneously, Mode I, exhibiting growth rates an order of magnitude lower, becomes unstable over a much broader frequency band, extending up to $F\approx 2 \times 10^{-3}$ . This modal behaviour shifts the critical Reynolds number $\textit{Re}_{cr}$ , i.e. $\alpha _{i}=0$ , towards lower values. The DNS domain, indicated by rectangular bars in figure 4, is placed inside the linearly unstable region, with the disturbance strip (coloured in cyan) positioned somewhat upstream of, or close to, $\textit{Re}_{cr}$ at the selected frequency $F$ . As shown in figure 4(c) for case Tw110, the DNS domain spans nearly the entire linearly unstable region of Mode II at $F_{{2\hbox{-}{\rm D}}}=124 \times 10^{-6}$ , consistent with our aim of investigating transition to turbulence triggered by Mode-II instability.

In the context of the 2-D DNS, a steady laminar solution is first obtained (see § 3.1), and LST-DNS comparisons of growth rate and phase speed are provided in Appendix C. The eigenfunctions, normalised by $\max \{|\hat {u}|\}$ , are successfully compared in figure 5(a,c) at $\textit{Re}=500$ for case Tw095, and at $\textit{Re}=650$ for case Tw110. In the subcritical regime, $u^{\prime }$ shows a phase jump near $y/\delta _{99,0} \approx 1$ , and $\rho ^{\prime }$ is confined near the wall around the critical layer $y=y_{c}$ , defined by $\bar {u}(y_{c})=c_{r}$ , resembling the ideal-gas case TadIG (not shown here). For Tw110, Mode II is affected by the pseudo-critical point (dashed green line at $y=y_{pc}$ in figure 5 c), with $\max \{|\hat {u}|\}$ located in the vapour-like regime. Figure 5(b,d) provide an overview of the modal instability in Tw095 and Tw110 via contours of $\rho ^{\prime }$ . In the latter case, the density fluctuations form ‘rope-shaped’ patterns around the pseudo-critical point, where gradients in transport and thermodynamic properties are largest (see figure 3), near the GIP (horizontal grey line), in agreement with Ren et al. (Reference Ren, Marxen and Pecnik2019). While similar density patterns occur for the second-mode instability in a hypersonic boundary layer (Unnikrishnan & Gaitonde Reference Unnikrishnan and Gaitonde2020), two key distinctions apply here: (i) transcritical Mode II is not linked to Mack’s second mode (Ren et al. Reference Ren, Marxen and Pecnik2019), and (ii) at $M_\infty =0.2$ , pressure perturbations are predominantly of hydrodynamic nature. Note that, contrary to Ren et al. (Reference Ren, Marxen and Pecnik2019), the wall-normal distribution of the pressure eigenfunction $\hat {p}$ peaks at the wall rather than at the pseudo-critical point (see figure 5 c).

Figure 5. Cases (a,b) Tw095 and (c,d) Tw110: (a,c) wall-normal eigenfunctions (lines DNS data, circles LST) of $u^{\prime }$ , $v^{\prime }$ , $ T^{\prime }$ , $\rho ^{\prime }$ and $p^{\prime }$ normalised by $\max \{|\hat {u}|\}$ at $\textit{Re}=500$ (Tw095) and $\textit{Re}=650$ (Tw110); (b,d) contours $\rho ^{\prime }$ normalised by their respective maxima. The locations of the pseudo-critical point $y=y_{pc}$ , i.e. where $\bar {T}^*=T^*_{pc}$ , the GIP $y=y_{\textit{GIP}}$ , and the critical layer $y=y_{c}$ are indicated in dashed green, dashed grey and dashed brown, respectively.

A key question concerning Mode-II instability, as it manifests under transcritical heating conditions, concerns its physical linear mechanism. Bugeat et al. (Reference Bugeat, Boldini, Hasan and Pecnik2024) demonstrated in a plane Couette flow that shear and baroclinic effects interact to generate two vorticity waves around the central layer, coinciding with the location of the minimum of kinematic viscosity and the critical layer. Similarly, the present 2-D boundary-layer flow is analysed using the linearised vorticity equation, where $\bar {\varOmega }=\partial \bar {v}/\partial x - \partial \bar {u}/\partial y$ is the base-flow vorticity, for the disturbance vorticity $\xi =\partial v^{\prime }/\partial x-\partial u^{\prime }/\partial y$ as

(4.1) \begin{equation} \underbrace {\dfrac {{\rm D} \xi }{{\rm D} t }}_{\textit{LHS}} \approx \underbrace {v^{\prime } \dfrac {\partial ^2 \bar {u}}{\partial y^2}}_{S_{\xi }} + \underbrace {\dfrac {\partial \bar {u}}{\partial y} \left ( \dfrac {\partial u^{\prime }}{\partial x} + \dfrac {\partial v^{\prime }}{\partial y} \right )}_{C_{\xi }} \, \underbrace {{}-\dfrac {1}{\bar {\rho }^2} \dfrac {\partial \bar {\rho }}{\partial y}\dfrac {\partial p^{\prime }}{\partial x}}_{B_{\xi }} + O(\mu ). \end{equation}

Here, $S_{\xi }$ , $C_{\xi }$ and $B_{\xi }$ denote the shear, compressible stretching and baroclinic terms, respectively. The viscous term is represented by $O(\mu )$ . The term $(1/\bar {\rho }^2)(\partial \rho ^{\prime }/\partial x)(\partial \bar {p}/\partial y)$ , which belongs to $B_{\xi }$ , is negligible (see Appendix B). Figure 6 evaluates (4.1) for case Tw110. For case Tw095, the only significant term of (4.1) far from the wall is the shear contribution $S_{\xi }$ , where $|\partial ^2 \bar {u}/\partial y^2|$ is maximal. In the spectral domain at $\textit{Re}=650$ , figure 6(a) confirms that both $S_{\xi }$ and $B_{\xi }$ reach a maximum near the pseudo-critical point, where the minimum of $\bar {\nu }$ is located. At the GIP, where the density-weighted vorticity $\varPhi =\bar {\rho }\bar {\varOmega }$ is maximal (figure 6 b), the compressible stretching term $C_{\xi }$ also peaks and remains large in the vapour-like region. The viscous term $O(\mu )$ exhibits a local maximum at the pseudo-critical point due to its direct dependence on $|\partial ^2 \bar {u}/\partial y^2|$ , which is largest at $y_{pc}$ . In figure 6(c), the sum of $S_{\xi }$ and $B_{\xi }$ exhibits two out-of-phase waves with phase difference $\pi$ (not shown), located around $y_{pc}$ and exhibiting asymmetry due to the structure of $B_{\xi }$ . When the out-of-phase $C_{\xi }$ term is included, as shown in figure 6(e), the two vorticity waves shift around the critical layer at $y=y_{c}$ . This behaviour is consistent with Bugeat et al. (Reference Bugeat, Boldini, Hasan and Pecnik2024), where $C_{\xi }$ was absent and the critical layer coincided with the GIP, i.e. the centreline of the plane Couette flow under the parallel flow assumption. In other words, the misalignment between the critical layer and the GIP in the boundary layer results in a corresponding shift of the two vorticity waves. With the addition of the viscous $O(\mu )$ term, the vorticity perturbation $\xi$ in figure 6( f) is further amplified near the wall, tilting the near-wall vorticity wave in the upstream direction.

Figure 6. Case Tw110. Terms of the vorticity perturbation (4.1): (a) spectral domain at $\textit{Re}=650$ (all terms are normalised by $\max \{|{\rm D} \xi /{\rm D} t|\}$ ); (b) normalised density-weighted vorticity $|\varPhi |=|\bar {\rho }\bar {\varOmega }|$ ; (c) $S_{\xi }+B_{\xi }$ ; (d) $C_{\xi }$ ; (e) $S_{\xi }+B_{\xi }+C_{\xi }$ ; and ( f) $\xi$ . The locations of the pseudo-critical point $y=y_{pc}$ , i.e. where $\bar {T}^*=T^*_{pc}$ , the GIP $y=y_{\textit{GIP}}$ , and the critical layer $y=y_{c}$ are indicated in dashed green, dashed grey and dashed brown, respectively.

4.2. Nonlinear evolution

We now focus on the nonlinear response of the boundary layer to a finite-amplitude perturbation, with the blowing–suction amplitude increased to $A_{{2\hbox{-}{\rm D}}}=7.5 \times 10^{-3}$ , compared to the infinitesimal amplitude used in § 4.1. Fourier modes are denoted using the double-spectral notation $(\omega /\omega _{{2\hbox{-}{\rm D}}},0)$ , where $\omega _{{2\hbox{-}{\rm D}}}$ is the fundamental frequency.

Figure 7(a) shows the downstream modal evolution for case Tw095. Once the $1\,\%$ threshold is crossed at $\textit{Re} \approx 600$ , the primary (1, 0) mode decays. Higher harmonics (2, 0), (3, 0) and beyond are slaved to the primary wave through weakly nonlinear effects in the receptivity region. Despite being linearly stable, they follow the streamwise growth, and subsequent stabilisation, of the forced (1, 0), with amplitudes approximately two orders of magnitude lower. Wall-normal profiles of streamwise velocity and density at $\textit{Re}=500$ and $700$ (figure 7 b,c) reveal only minor nonlinear effects, with the boundary-layer receptivity resembling the incompressible TadIG case.

Figure 7. Case Tw095 with $A^{(1,0)}_{{2\hbox{-}{\rm D}}}=7.5 \times 10^{-3}$ : (a) maximum wall-normal mass-flux amplitude for mode (1, 0) (solid line), (2, 0) (dash-dotted line), and (3, 0) (dashed line); (b) streamwise velocity and (c) density perturbations as functions of the wall-normal coordinate $y/\delta _{99,0}$ at $\textit{Re}=500$ (solid line) and $\textit{Re}=700$ (dash-dotted line), normalised by their respective $\max \{|\hat {u}^{(1,0)}|\}$ . In (b,c), the scaled LST solution is represented with circle symbols at $\textit{Re}=500$ , and with square symbols at $\textit{Re}=700$ .

Figure 8. Case Tw110, with (a,c) $A^{(1,0)}_{{2\hbox{-}{\rm D}}}=7.5 \times 10^{-4}$ and (b,d) $A^{(1,0)}_{{2\hbox{-}{\rm D}}}=7.5 \times 10^{-3}$ , for: (a,b) maximum wall-normal mass-flux amplitude for mode $(1,0)$ (solid line), $(2,0)$ (dash-dotted line), $(3,0)$ (dashed line), $(4,0)$ (dotted line), $(5,0)$ (solid line with triangles), and $(6,0)$ (solid line with diamonds); (c,d) phase speed $c_{r}$ for modes $(1,0)$ (solid line) and $(2,0)$ (dash-dotted line). In (b), the mean-flow distortion (MFD) $(0,0)$ is indicated with a black solid line. The LST solution is represented with circle symbols for mode $(1,0)$ , square symbols for mode $(2,0)$ , and asterisk symbols for mode $(3,0)$ .

In contrast to the subcritical regime, nonlinearity significantly influences the 2-D modal evolution in case Tw110. To highlight this, $A_{{2\hbox{-}{\rm D}}}$ is increased from $7.5 \times 10^{-4}$ in figure 8(a) to $7.5 \times 10^{-3}$ in figure 8(b). In the low-amplitude case, modes $(1,0)$ and $(2,0)$ follow the LST prediction up to $\textit{Re} \approx 700$ , and as shown in figure 4(c), linear instability also arises at $2 F_{{2\hbox{-}{\rm D}}}$ , with a larger growth rate than that of $(1,0)$ . Higher harmonics ( $\omega /\omega _{{2\hbox{-}{\rm D}}}\geq 3$ ) deviate from their linear evolution, growing rapidly to high amplitudes – unlike in the subcritical regime (figure 7 a). This behaviour resembles the Kelvin–Helmholtz (KH) instability in a mixing layer (Babucke, Kloker & Rist Reference Babucke, Kloker and Rist2008). In the high-amplitude case in figure 8(b), higher harmonics exceed the $0.1\,\%$ amplitude threshold at lower values of $\textit{Re}$ . Mode $(2,0)$ , which is more unstable than mode $(1,0)$ in the linear regime, surpasses $(1,0)$ at $\textit{Re} \approx 600$ , before nonlinearly saturating after reaching the $2\,\%$ amplitude threshold. As it peaks at $\textit{Re} \approx 700$ , a subharmonic resonance mechanism with its subharmonic mode $(1,0)$ , typically observed as vortex pairing in mixing layers (Monkewitz Reference Monkewitz1988), emerges, destabilising the latter. Thus mode $(1,0)$ undergoes strong destabilisation and becomes the most dominant mode again further downstream. At this streamwise location, the MFD, i.e. mode $(0,0)$ , reaches up to $5\,\%$ , and all higher harmonics are fully nonlinear. It is worth noting that, contrary to the conclusions of Kachanov (Reference Kachanov1994) – who reported that a threshold amplitude of the 2-D fundamental TS wave of the order of $29\,\%$ of the mean flow was required to amplify 2-D subharmonics – the amplitude of mode $(2,0)$ reaches only $5\,\%$ in the present case under transcritical conditions. The resonant interaction between modes $(1,0)$ and $(2,0)$ is further evidenced by their phase speeds $c_{r}$ in figure 8(c,d). For the low-amplitude case ( $A^{(1,0)}_{{2\hbox{-}{\rm D}}}=7.5 \times 10^{-4}$ ), the phase speeds are sufficiently close near $\textit{Re} \approx 800$ . In contrast, for the high-amplitude case, $c^{(1,0)}_{r}$ and $c^{(2,0)}_{r}$ remain nearly identical across the entire computational domain, triggering subharmonic resonance at $\textit{Re}\approx 700$ . Notably, at $\textit{Re} \approx 730$ , $c^{(2,0)}_{r}$ increases sharply, deviating from its linear evolution (in grey) and indicating strong nonlinear stabilisation (figure 8 b).

To qualitatively assess the nonlinear evolution of Mode II, selected perturbation contours for $A^{(1,0)}_{{2\hbox{-}{\rm D}}}=7.5 \times 10^{-3}$ are shown in figure 9(a,b). The region corresponding to highest specific heat at constant pressure, between $98\,\% \max \{c_{p}\}$ and $\max \{c_{p}\}$ , is shaded in green. At this high disturbance level, the constant-pressure assumption for a laminar boundary layer breaks down, and $\max \{c_{p}\}$ depends on both reduced pressure and temperature. It is estimated using the analytical Widom-line relation $T_{r,pc}=1/A_{\textit{VdW}} \ln (p_{r}) + 1$ , where the Widom line, defined as the locus of pseudo-critical points at supercritical pressure (figure 2), is characterised by the critical slope $A_{\textit{VdW}}=T_{c}/p_{c} (\mathrm{d}p/\mathrm{d}T)_{c}=4$ (Banuti Reference Banuti2015). Unlike in the linear regime (figure 5 c,d), where the pseudo-critical point height grows with $\sqrt {x}$ , the nonlinear disturbance wave propagation leads to increasing distortion of the Widom line, proportional to the perturbation magnitude (figure 9 a,b). The large harmonic mode $(2,0)$ , prominent between $\textit{Re}\approx 600$ and $\textit{Re}\approx720$ , halves the streamwise wavelength of the perturbation wave ( $(1,0)$ as the fundamental), thereby doubling the oscillation frequency of the Widom line. As shown in the inset of figure 9(a), the crests and troughs of the Widom line align with regions of positive and negative wall-normal velocity perturbation $v^{\prime }$ . Since $v^{\prime }$ and the pressure perturbation $p^{\prime }$ are out of phase (Luhar, Sharma & McKeon Reference Luhar, Sharma and McKeon2014; Bugeat et al. Reference Bugeat, Boldini, Hasan and Pecnik2024), the upward and downward displacements of the Widom line correspond to local pressure decreases and increases, respectively. Density perturbations in figure 9(b), initially confined near the Widom line in the linear regime (figure 5 d), now reach up to $20\,\%$ of the free-stream value in regions where $p^{\prime }_{r}$ is negative, i.e. closer to the critical point. These perturbations mirror the billowing behaviour of the Widom line, with $\rho ^{\prime }\lt 0$ around the Widom-line crests, and $\rho ^{\prime }\gt 0$ around the troughs.

Figure 9. Case Tw110. Instantaneous contours at $T/T_0=0$ , where $T_0=2\pi /\omega _0$ (fundamental frequency $\omega _0$ ), with $A^{(1,0)}_{{2\hbox{-}{\rm D}}}=7.5 \times 10^{-3}$ : (a) reduced pressure fluctuation $p^\prime _{r}=p^{*\prime }/p^*_{c}$ , (b) density fluctuation $\rho ^\prime$ , (c) vorticity $\varOmega$ , (d) streamwise velocity $u$ , with boundary-layer thickness $\delta _{99}$ and displacement thickness $\delta _1$ indicated by dotted and dashed lines, respectively, and (e) Mach number $M=u/a$ . The Widom line $y=y_{\textit{WL}}$ lies within the green region, i.e. between $98\,\% \max \{c_{p}\}$ and $\max \{c_{p}\}$ . Note that the Widom line is used here as a spatial reference for the local pseudo-critical point at supercritical pressure. Insets in (a–c) show the wall-normal velocity fluctuation $v^\prime$ , density $\rho$ , and vorticity fluctuation $\xi$ , respectively. The inset in (d) highlights separation zones ( $u \lt 0$ ) in blue and includes velocity vectors $|\boldsymbol{V}|=\sqrt {u^2+v^2}$ . A supplementary movie of the billow roll-ups is available in the supplementary material is available at https://doi.org/10.1017/jfm.2025.10993.

The term ‘billowing’ is used by analogy with classical shear-layer flows, where a periodic train of large, rolling wave-like structures forms as the KH instability evolves nonlinearly (Klaassen & Peltier Reference Klaassen and Peltier1985; Liu, Kaminski & Smyth Reference Liu, Kaminski and Smyth2023). Here, a train of streamwise-growing density billows arises (see inset of figure 9 b). In both shear layers and the present case, billowing arises from an excess of vorticity concentrated in a localised flow region (see figure 6 b). Although the train of billowing flow patterns is caused by the transcritical Mode-II instability rather than by the KH instability, Bugeat et al. (Reference Bugeat, Boldini, Hasan and Pecnik2024) proved that the resulting vorticity fields are identical in the linear regime. Thus it is not surprising that in the nonlinear regime, the billowing patterns observed here resemble those seen in classical KH-type roll-ups, despite the absence of a true shear-layer mechanism.

Figure 9(b,c) highlight the nonlinear evolution of these billows via contours of $\rho ^{\prime }$ and vorticity $\varOmega =\partial v/\partial x- \partial u/\partial y$ . Their deformation and streamwise growth progressively narrow the near-wall vapour-like region, which shifts significantly closer to the wall at $\textit{Re} \approx 800$ compared to the linear regime. Simultaneously, the upper vapour-like fluid is lifted up by the rolling billows. Peaks in $\varOmega$ , which were previously aligned with the Widom line in the laminar regime, now appear below the Widom-line troughs in the near-wall region and at the wall. Starting from $\textit{Re} \approx 780$ , and recurring periodically downstream, the $\varOmega$ -peaks strengthen (high-shear regions), and their magnitude scales with their fluctuation $\xi$ (inset of figure 9 c). The contours of streamwise velocity in figure 9(d) reveal the emergence of near-wall flow reversal regions ( $u\lt 0$ ), located just beneath the near-wall $\varOmega$ -peaks at the Widom-line troughs, indicating the formation of a developing shear layer near the wall. In the vicinity of these flow reversal zones, a region of low streamwise velocity forms, leading to the reduction of the boundary-layer thickness $\delta _{99}$ by approximately $10\,\%$ relative to the unperturbed profile. At $\textit{Re} \approx 800$ , the near-wall separation zone reaches a height of approximately $2\,\%$ of $\delta _{99}$ , while the billow crest extends beyond $y/\delta _{99,0}=1$ . Here, the local Mach number reaches its maximum, as shown in figure 9(e). This behaviour results from both the boundary-layer displacement effect and the reduction in the local speed of sound at the Widom line. Downstream of the crest, the Mach number rapidly returns to its laminar value.

Figure 10. Case Tw110 for $A^{(1,0)}_{{2\hbox{-}{\rm D}}}=7.5 \times 10^{-3}$ . Wall-normal slice at $\textit{Re}=802$ showing: (a,b) instantaneous streamwise velocity $u$ and reduced specific heat at constant pressure $c_{p,r}/c_{p,r,\infty }$ at time periods $t/T_0=0,0.25,0.5,0.95$ , where $T_0=2\pi /\omega _0$ is the fundamental forcing period; (c,d) time-averaged streamwise velocity $\langle u \rangle$ and density $\langle \rho \rangle$ profiles. In (c,d), the r.m.s. of $u^{\prime }$ and $\rho ^{\prime }$ , respectively, and higher harmonics (modes $(1,0)$ , $(2,0)$ and $(3,0)$ ) are shown. The locations of the GIP and inflection point at $t/T_0=0$ are marked in (a,b) by grey circle and purple star symbols, respectively.

To further investigate the near-wall high-shear region during billow roll-up, wall-normal profiles are extracted at $\textit{Re}=802$ , corresponding to the first flow reversal found in figure 9(d). Figure 10(a,b) show the streamwise velocity $u$ and the reduced specific heat at constant pressure $c_{p,r}/c_{p,r,\infty }$ , respectively, at four time instants within one forcing period $T_0=2\pi /\omega _0$ . Between $t/T_0=0.95$ and $t/T_0=0$ , the billow roll-up induces two distinct kinks in the $u$ profile: one at the billow crest at $y/\delta _{99,0}\approx 1.0$ , and another at the trough near $y/\delta _{99,0}\approx 0.3$ . These locations coincide with sharp peaks in $c_{p}$ (figure 10 b), with the highest values occurring at $t/T_0=0$ , linked to a local pressure drop $p^{\prime }_r\lt 0$ (see figure 9 a). Beneath the near-wall kink in the vapour-like region, flow reversal at the wall is observed (inset of figure 10 a), along with a GIP, where the density-weighted vorticity $\varPhi =\rho \varOmega$ is maximal (grey circle). Here, the vorticity experiences a maximum, with an inviscid instability found according to the generalised Fjørtoft’s criterion (Bugeat et al. Reference Bugeat, Boldini, Hasan and Pecnik2024). By $t/T_0=0.25$ , no further GIP or flow reversal is observed, confirming the periodic nature of the billows generation. At $t/T_0=0.50$ , the profiles resemble those of the unperturbed boundary layer (figure 3). The large near-wall shift in the instantaneous  $u$ observed over one $T_0$ is caused by strong near-wall velocity fluctuations, reaching up to $5\,\%$ of the root mean square (r.m.s.) value, (figure 10 c), sustained by both the primary mode $(1,0)$ and its first harmonic $(2,0)$ . A secondary peak in $u^{\prime }_{\textit{rms}}$ at $y/\delta _{99,0} \approx 0.7$ is primarily attributed to mode $(1,0)$ . Density fluctuations (figure 10 d) are dominant in the high- $c_{p}$ region, with the $\rho ^{\prime }_{\textit{rms}}$ -peak located at the height of the far-wall $u$ kink.

In conclusion, under transcritical conditions, we observe localised flow reversal beneath the billow roll-up at $\textit{Re}=802$ , along with the formation of a shear layer just above it. This behaviour is linked to the sharp near-wall increase in $c_{p}$ , which amplifies near-wall $u$ disturbances and destabilises the local mean flow. In other words, this effect is tied to the upward and downward displacements of the Widom line, which correspond to decreases and increases in reduced pressure $p_{r}$ , respectively. As $p_{r}$ decreases, steep gradients in thermodynamic properties intensify. Ultimately, the initial assumption of a ZPG laminar boundary layer breaks down. When the billow rolls up, a net $\mathrm{d}p/\mathrm{d}x\lt 0$ (APG) forms. Simultaneously, the MFD is insufficient to counteract the large near-wall (Mode II) $u$ fluctuations, leading to a travelling region of flow reversal. While Taylor (Reference Taylor1936) attributed transition in incompressible flows to unsteady pressure gradients caused by ‘external’ disturbances leading to separation (see note by Kloker Reference Kloker2024), similar localised separation zones, such as those in figure 9(d), were found under strong APG in the K-type breakdown of an ideal-gas boundary layer (Kloker Reference Kloker1993; Kloker & Fasel Reference Kloker and Fasel1995) and in the laminar-to-turbulent transition experiments of Kosorygin (Reference Kosorygin1994). In the current study, however, the disturbance is ‘internal’, originating from the displacement of the Widom line, and is absent under weakly non-ideal-gas conditions at the same forcing amplitude.

5.1. Modal analysis

Figure 11 presents the relevant modes for the subcritical cases with infinitesimal and finite amplitude, Tw095-IA and Tw095-LA, and the transcritical cases with infinitesimal and finite amplitude, Tw110-IA and Tw110-LA. All share the same fundamental frequency and spanwise wavenumber; the ‘LA’ cases in figure 11(b,c) use the same 2-D and 3-D forcing amplitudes, while the ‘IA’ cases in figure 11(a,d) have infinitesimal 3-D forcing amplitude $A_{{3\hbox{-}{\rm D}}}=1.33 \times 10^{-6} A_{{2\hbox{-}{\rm D}}}$ (see table 2).

Figure 11. Streamwise evolution of the $y$ -maximum of $(\rho u)^\prime$ for the most relevant modes $( \omega / \omega _{{2\hbox{-}{\rm D}}}, \beta / \beta _0)$ for cases (a) Tw095-IA, (b) Tw095-LA, (c) Tw110-LA, and (d) Tw110-IA. The minimum and maximum values of the time- and spanwise-averaged skin-friction coefficient are indicated as $C_{\!\textit{f},\textit{min}}=\min \{C_{\!{f}}\}$ and $C_{\!\textit{f},\textit{max}}=\max \{C_{\!{f}}\}$ , respectively. Insets in (c) and (d) highlight the relevant modes in the breakdown region. Note the different $\textit{Re}_{x}$ -axis limits between the subcritical (a,b) and (c,d) transcritical cases.

In case Tw095-IA (figure 11 a), no transition occurs. Although the oblique mode $(1/2,1)$ undergoes secondary subharmonic growth starting at $\textit{Re}_{x}/10^5 \approx 3$ , its amplitude remains insufficient to trigger nonlinear interactions. The primary 2-D mode $(1,0)$ decays in line with figure 7(a), while mode $(0,2)$ continues to grow due to transient growth (Boldini et al. Reference Boldini, Bugeat, Peeters, Kloker and Pecnik2024b ).

In contrast, case Tw095-LA (figure 11 b) transitions to turbulence via the subharmonic instability, which rapidly induces three-dimensionality. Its modal evolution closely resembles that of case TadIG (not shown). At the disturbance strip, higher harmonics such as $(2,0)$ and $(3/2,1)$ are also triggered but remain only moderately amplified. When the primary wave $(1,0)$ reaches ${\sim} 1\,\%$ at $\textit{Re}_{x}/10^5 \approx 2.9$ , subharmonic resonance sets in, and the forced subharmonic oblique mode $(1/2,1)$ , absent in the 2-D DNS set-up in figure 7(a), grows rapidly, surpassing $(1,0)$ near $\textit{Re}_{x}/10^5 \approx 4$ . Compared to TadIG, this process is delayed due to the slower destabilisation of $(1,0)$ as the wall temperature is increased towards the Widom line. During the secondary instability growth, higher modes experience strong nonlinear amplification, contributing to the increase of the MFD $(0,0)$ near $C_{\!\textit{f},\textit{min}}=\min \{C_{\!{f}}\}$ . Later, all modes saturate except for $(1/2,3)$ , which matches the amplitude of $(1/2,1)$ at $\textit{Re}_{x}/10^5 \approx 5.3$ . Hereafter, a wave–vortex triad typical of the oblique breakdown (Chang & Malik Reference Chang and Malik1994) forms between the strongly nonlinear $(1/2,3)$ and steady vortex mode $(0,2)$ , further destabilising $(1/2,1)$ (see inset of figure 11 b). Further downstream where nonlinear saturation sets in, mode $(0, 0)$ reaches its maximum, marking the onset of turbulence at $C_{\!\textit{f},\textit{max}}=\max \{C_{\!{f}}\}$ .

For case Tw110-LA in figure 11(c), although the evolution of the fundamental wave and its higher harmonics closely follows the 2-D development in figure 8(b) up to $\textit{Re}_{x}/10^5 \approx 6.5$ , the subharmonic resonance is significantly delayed due to the strong 2-D nonlinearity of the higher harmonics. The oblique mode $(1/2,1)$ is initially strongly damped between $\textit{Re}_{x}/10^5 \approx 1.0$ and $2.0$ , and maintains a low amplitude level, along with $(3/2,1)$ , up to $\textit{Re}_{x}/10^5 \approx 4.7$ . Before mode $(1,0)$ undergoes subharmonic resonance with $(2,0)$ (see § 4.2), it saturates and a phase-speed-locking process (cf. Hader & Fasel Reference Hader and Fasel2019) occurs between mode $(1,0)$ and its secondary wave $(1/2,1)$ , triggering the growth of $(1/2,1)$ via the secondary-instability mechanism (Herbert Reference Herbert1988). Mode $(1/2,1)$ and its higher harmonics grow to nonlinear amplitude levels and saturate at $\textit{Re}_{x}/10^5 \approx 6.0$ . On the contrary, modes with spanwise-wavenumber parameter $ \beta / \beta _0=1$ briefly grow near $C_{\!\textit{f},\textit{min}}$ due to nonlinear interaction, but do not contribute to the later breakdown stage. When the amplitudes of subharmonic modes ( $\omega / \omega _0=1/2$ ) approach those of the primary wave and higher harmonics, additional higher modes with $ \beta / \beta _0 \geq 2$ are nonlinearly amplified, rapidly growing downstream and enhancing both the destabilisation of $(1,0)$ and the MFD. In the early breakdown stage, following the saturation of mode $(1,0)$ at $\textit{Re}_{x}/10^5 \approx 7.1$ , mode $(1/2,1)$ grows nonlinearly, surpassing mode $(1,0)$ at $\textit{Re}_{x}/10^5 \approx 7.35$ , while the MFD $(0,0)$ exceeds the $20\,\%$ amplitude threshold. In the late transitional stage, before $C_{\!\textit{f},\textit{max}}$ at $\textit{Re}_{x}/10^5 \approx 8.5$ , mode $(1/2,1)$ remains dominant, eventually followed by the abrupt growth of $(1/2,3)$ at $\textit{Re}_{x}/10^5 \approx 9.0$ .

Figure 11(d) illustrates the modal evolution of case Tw110 with infinitesimal 3-D forcing (‘IA’). As expected, up to $\min \{C_{\!{f}}\}$ , the evolution of modes $(h,0)$ , with $h \geq 0$ , matches that of case Tw110-LA in figure 11(c). No subharmonic resonance of $(1/2,1)$ occurs further downstream in this case, as no phase-speed locking mechanism between the primary wave and its subharmonic is present. In contrast, fundamental resonance with $(1,1)$ sets in at $\textit{Re}_{x}/10^5 \approx 6.2$ , along with the amplification of 3-D harmonic modes with spanwise-wavenumber parameter $\beta /\beta _0 = 2$ , before mode $(1/2,1)$ grows nonlinearly only from $\textit{Re}_{x}/10^5 \approx 6.5$ onwards. Around $\textit{Re}_{x}/10^5 \approx 7.2$ , the spectrum rapidly fills up with all other nonlinearly amplified 3-D modes, contributing to a significant increase in the MFD $(0,0)$ . Note that mode $(1,1)$ arises initially due to numerical background noise, and the forcing of $(1/2,1)$ remains likewise at the noise level; therefore, both are insignificant and not representative of classical boundary-layer modes in the early nonlinear stage. Further downstream around $\textit{Re}_{x}/10^5 \approx 6.0$ , both modes $(1,1)$ and $(0,1)$ (forming a wave–vortex triad with $(1,0)$ ) begin to behave as discrete boundary-layer wave modes, whereas $(1/2,1)$ continues to be a continuous, non-resonant mode. Under such numerical background noise conditions, and in the absence of primary instability of the investigated 3-D modes, the fundamental-resonance mechanism proves to be more robust than the subharmonic one. Eventually, the steady mode $(0,1)$ , which dominates the transitional stage in the K-type breakdown (Boldini et al. Reference Boldini, Bugeat, Peeters, Kloker and Pecnik2024a ), reaches the amplitude level of $(1,0)$ . In the final breakdown stage, $(1,1)$ and $(1/2,3)$ alternate as the dominant modes before $C_{\!\textit{f},\textit{max}}$ is reached. Compared to Tw110-LA, the transition is more gradual, but follows a K-type breakdown scenario. Notably, despite the infinitesimal 3-D forcing, breakdown to turbulence is triggered solely by the 2-D fundamental wave and numerical background noise – a behaviour not observed in case Tw095-IA under weakly non-ideal-gas conditions.

5.2. Flow structures

To investigate the mechanisms driving the breakdown in both thermodynamic regimes, figure 12 displays the streamwise evolution of the relevant transitional flow structures using isocontours of the $Q$ -criterion.

Figure 12. Instantaneous isosurfaces of the $Q$ -criterion, coloured by the streamwise velocity magnitude: (a) case Tw095-LA ( $Q = 0.015$ ) at $t/T_0=0$ , (b) case Tw110-LA ( $Q = 0.020$ ) at $t/T_0=0.5$ , and (c) case Tw110-IA ( $Q = 0.020$ ) at $t/T_0=0.5$ . Here, $T_0$ is the period of the fundamental wave. The side $x{-}y$ plane shows the instantaneous spanwise vorticity $\omega _{z}$ . Isosurfaces of the separation zones, i.e. regions with $u\lt 0$ , are coloured in cyan. For better visualisation, the domain is copied twice in the spanwise direction. Supplementary movies are available in the supplementary material is available at https://doi.org/10.1017/jfm.2025.10993.

The subcritical case Tw095-LA (figure 12 a) exhibits the classical H-type breakdown, similar to the ideal-gas case TadIG, with staggered $\Lambda$ -vortices and high-shear layers (Bake et al. Reference Bake, Meyer and Rist2002) above the vortices’ tips, characterised by peaks in $\omega _{z}\propto\partial u/\partial y$ . Conversely, no vortices are present at $\lambda _{z}/2$ , i.e. half a spanwise wavelength apart. While the early breakdown stage closely resembles that of the ideal-gas case TadIG, the $\Lambda$ -vortices become more elongated in the streamwise direction from $\textit{Re}_{x}/10^5 \approx 5.2$ onwards, coinciding with the large amplitude of 3-D modes $(3/2,1)$ and $(0,2)$ in figure 11(b). This elongation of the $\Lambda$ -vortices, along with longitudinal structures on their sides, is particularly visible for $5.2 \leq Re_{x}/10^5 \leq 6$ in figure 13(a), which shows a horizontal flow snapshot ( $x{-}z$ plane) of the streamwise velocity at $y/\delta _{99,0}=0.49$ . In contrast, the staggered $\Lambda$ -vortex arrangement in case TadIG (figure 13 b) breaks down more rapidly during the transitional stage.

Figure 13. Contours of instantaneous streamwise velocity ( $x{-}z$ plane at $y/\delta _{99,0}=0.49$ ): (a) Tw095-LA and (b) TadIG. The dashed and dash-dotted white vertical lines represent the locations of $C_{\!\textit{f},\textit{min}}=\min \{C_{\!{f}}\}$ and $C_{\!\textit{f},\textit{max}}=\max \{C_{\!{f}}\}$ , respectively. For better visualisation, the domain is copied once in the spanwise direction.

Compared to the subcritical regime, figure 12(b,c) reveal distinct vortical structures in the transcritical heating regime. During the early stage of breakdown, dominated by mode $(1,0)$ , 2-D billow structures with near-wall separation zones (highlighted in cyan) are convected downstream. At later stages, cases Tw110-LA and Tw110-IA exhibit fundamentally different breakdown dynamics. In Tw110-LA, the forced 3-D mode $(1/2,1)$ rapidly grows from $0.01\,\%$ to $13\,\%$ between $\textit{Re}_{x}/10^5 \approx 6$ and $7.3$ , leading to the sudden appearance of staggered $\Lambda$ -vortices lifting off from the billow roll-ups. Unlike the elongated staggered alignment observed in Tw095-LA, $\Lambda$ -vortices here emerge more abruptly and locally, with secondary structures forming at half a spanwise wavelength apart (see § 5.2.1). By $\textit{Re}_{x}/10^5 \approx 7.5$ , where mode $(0,0)$ reaches approximately $30\,\%$ , the first row of hairpin vortices becomes visible, with small-scale structures at their legs. Farther downstream, the complete breakdown of the preceding 2-D billow roll-up becomes evident, characterised by near-wall longitudinal structures between $\textit{Re}_{x}/10^5 \approx 7.5$ and $7.7$ .

Figure 14. Case Tw110-LA. Instantaneous isosurfaces of (a) spanwise vorticity $|\omega _{z}|=0.45$ and (b) streamwise vorticity $|\omega _{x}|=0.45$ , coloured by the streamwise momentum $\rho u$ magnitude at $t/T_0=0.5$ . Here, $T_0$ is the period of the fundamental wave. For better visualisation, the domain is copied once in the spanwise direction.

In contrast, case Tw110-IA exhibits a delayed breakdown of each 2-D billow, as indicated by the extended near-wall separation zones. The onset of three-dimensionality is significantly delayed, consistent with the more gradual growth of 3-D modes in figure 11(d). The $\Lambda$ -like vortices emerge within the 2-D billow roll-ups, initially displaying an aligned peak–valley splitting, which is characteristic of the fundamental K-type breakdown (see § 5.2.2). The $\Lambda$ -vortices exhibit a shorter streamwise than spanwise wavelength, consistent with the large amplitude of mode $(2,1)$ in the early breakdown stage in figure 11(d). No strong lift-up of hairpin-like vortices is observed; however, near-wall longitudinal structures, resembling those in Tw110-LA, gradually emerge. Ultimately, each billow independently breaks into small-scale turbulent structures, eventually merging with the turbulent structures originating from the previous billow. From $\textit{Re}_{x}/10^5 \approx 8.0$ onwards, as mode $(0,0)$ reaches amplitude ${\sim} 30\,\%$ , the resulting turbulent flow pattern resembles that of Tw110-LA between $\textit{Re}_{x}/10^5 \approx 7.5$ and $7.7$ , as shown in figure 12(b).

5.2.1. Case Tw110-LA: detailed breakdown analysis

Figure 14 shows close-up views of the initial breakdown region in case Tw110-LA ( $6.90 \leq Re_{x}/10^5 \leq 7.52$ ), with isosurfaces of spanwise $\omega _{z}$ and streamwise $\omega _{x}$ vorticity in figure 14(a,b), respectively, at the same time instant as in figure 12(b). The visualisation spans two spanwise wavelengths to capture the full structure. The $\omega _{z}$ distribution reveals a delta-wing-like shear layer at $z=0.5\lambda _{z},1.5\lambda _{z},\ldots$ (dashed black lines), with peak spanwise vorticity at its tip. These structures, featuring round, slightly inclined heads, are aligned above the oppositely signed $\omega _{x}$ legs of the $\Lambda$ -shaped vortex (figure 14 b). At these locations, hereafter referred to as ‘peak’ locations, the strongest $u^{\prime }$ and $\rho ^{\prime }$ perturbations occur with $u^{\prime }_{\textit{rms}}\approx 0.2$ and $\rho ^{\prime }_{\textit{rms}}\approx 0.18$ . This layer is hereafter referred to as the ‘upper’ high-shear layer and is denoted as ‘UL’. Half a spanwise length apart, at $z=0,\lambda _{z},\ldots$ (dashed red line), $u^{\prime }$ and $\rho ^{\prime }$ , unlike in case Tw095-LA, do not decay but exhibit the second-largest r.m.s.amplitudes in the spanwise direction. This location is thus hereafter referred to as the ‘co-peak’; cf. the ideal-gas APG case (Kloker Reference Kloker1993; Kosorygin Reference Kosorygin1994; Kloker & Fasel Reference Kloker and Fasel1995). Here, the adjacent $\Lambda$ -vortices are more closely spaced than in cases TadIG and Tw095-LA, consistent with a delta-wing-like topology (Kloker & Fasel Reference Kloker and Fasel1995), and a secondary streamwise vortical structure develops between them.

Figure 15. Case Tw110-LA. Instantaneous contours of spanwise vorticity $\omega _{z}$ in the $\textit{Re}_{x}$ – $y/\delta _{99,0}$ plane at (a,b) $t/T_0=0$ , (c,d) $t/T_0=0.25$ , (e, f) $t/T_0=0.5$ , (g,h) $t/T_0=0.75$ , and (i, j) $t/T_0=1.0$ . Here, $T_0$ is the period of the fundamental wave. The first column (a,c,e,g,i) corresponds to the spanwise ‘co-peak’ location at $z/\lambda _{z}=0$ (see figure 14), while the second column (b,d, f,h, j) corresponds to the spanwise ‘peak’ location at $z/\lambda _{z}=0.5$ (see figure 14). The ‘upper’ and ‘inverted lower’ high-shear layers are labelled as ‘UL’ and ‘IL’, respectively. The near-wall region for which $u\lt 0$ is coloured in cyan. The Widom line $y=y_{\textit{WL}}$ lies within the green region, i.e. between $95\,\% \max \{c_{p}\}$ and $\max \{c_{p}\}$ .

The origin and evolution of this secondary vortex system is illustrated in figure 15, which shows instantaneous contours of $\omega _{z}$ in longitudinal planes (between $\textit{Re}_{x}/10^5 = 6.90$ and $7.52$ ) at five time steps ( $t/T_0=0,0.25,0.5,0.75,1.0$ ). At the ‘co-peak’ plane ( $z/\lambda _{z}=0$ , figure 15 a,c,e,g,i), a local maximum of $\omega _{z}$ forms near the wall, hereafter referred to as the ‘inverted lower’ high-shear layer and labelled as ‘IL’, resembling that of the ideal-gas APG case. Its origin is linked to the downstream-convected near-wall separation zone (highlighted in cyan) below the 2-D billow trough, which induces a local high-shear velocity profile (see figure 10 a). This layer is unstable to 3-D perturbations (cf. the stationary laminar separation bubble; Maucher et al. Reference Maucher, Rist, Kloker and Wagner2000) due to strong near-wall  $u^{\prime }$ disturbances (see inset in figure 15 c), which are transported downstream with the shear layer. Similar to the 2-D investigation (figure 10 c), they correspond to a strong near-wall amplitude of mode $(1,0)$ and its higher harmonics. Mode $(1/2,1)$ , which was absent in § 4.2, reaches its peak at $y/\delta _{99,0}\approx 1$ (not shown). Over time, the ‘inverted lower’ layer evolves: its head is pushed towards the wall, while its leg is lifted upwards (figure 15 e). The secondary vortex system shown in figure 14(b) then develops above it, again reminiscent of the ideal-gas APG case. The two oppositely signed vortical structures transport low-velocity, low-density fluid upwards, as seen in the $\rho u$ contours in figure 14(b). At later times (figure 15 g,i), the ‘inverted lower’ layer merges with a second high-shear layer, crossing the entire boundary layer.

Figure 16. Case Tw110-IA. Instantaneous isosurfaces of (a) spanwise vorticity $|\omega _{z}|=0.45$ and (b) streamwise vorticity $|\omega _{x}|=0.45$ , coloured by the streamwise momentum $\rho u$ magnitude at $t/T_0=0.5$ . Here, $T_0$ is the period of the fundamental wave. For better visualisation, the domain is copied once in the spanwise direction.

At $t/T_0=1$ , the flow structures at the ‘co-peak’ location match those at the ‘peak’ location at $t/T_0=0$ in figure 15(b), indicating the staggered, periodic breakdown pattern of case Tw110-LA. The newly formed shear layer initially consists of two adjacent, streamwise-aligned, counter-rotating vortex filaments. As it becomes increasingly unstable in figure 15(b,d, f), it rolls up and breaks into small-scale structures (see the secondary vortex system at ‘peak’ locations near $\textit{Re}_{x}/10^5 \approx 7.35$ in figure 14). Simultaneously, the ‘upper’ high-shear layer (‘UL’ in figure 15 b,d, f) also destabilises and breaks up into numerous small vortical structures (figure 15 h, j), generating an isolated hairpin-like vortex, captured by the $\omega _{z}$ isocontours at $\textit{Re}_{x}/10^5 \approx 7.5$ in figure 14(a). Note that the breakdown to turbulence is initiated by the earlier break-up of the ‘inverted lower’ layer at the ‘co-peak’ positions, which precedes the break-up of the ‘upper’ layer at the ‘peak’ positions. This sequence is also evident in figure 12(b), where near-wall turbulence triggered by the ‘inverted lower’ layer emerges significantly farther upstream than the small-scale structures in the boundary-layer core, which are caused by the break-up of the ‘upper’ layer.

5.2.2. Case Tw110-IA: detailed breakdown analysis

Figure 16(a,b) show isocontours of $\omega _{z}$ and $\omega _{x}$ for case Tw110-IA, respectively, at the same time step as in figure 12(c). The visualisation spans two spanwise wavelengths to capture the full structure. Unlike case Tw110-LA, mode $(1,0)$ retains a high amplitude over a longer streamwise distance, allowing the 2-D billow roll-ups to grow farther downstream, with intensified billow crests and larger near-wall separation zones.

Figure 17. Case Tw110-IA. Instantaneous contours of streamwise vorticity $\omega _{x}$ and velocity $u^*/u^*_{\infty }$ in the $z/\delta _{99,0}$ – $y/\delta _{99,0}$ plane at (a,b) $t/T_0=0.5$ ( $\textit{Re}_{x}/10^5 = 7.62$ ), (c,d) $t/T_0=0.75$ ( $\textit{Re}_{x}/10^5 = 7.70$ ), and (e–h) $t/T_0=1.0$ . Here, $T_0$ is the period of the fundamental wave. The dashed black lines in (a,c,e,g) indicate contours of $|\omega _{z}|=0.45$ , while the black lines in (b,d, f,h) correspond to $\delta _{99}$ . The ‘upper’ and ‘inverted lower’ high-shear layers are labelled as ‘UL’ and ‘IL’, respectively. The near-wall region for which $u\lt 0$ is coloured in cyan. The Widom line $y=y_{\textit{WL}}$ lies within the green region, i.e. between $95\,\% \max \{c_{p}\}$ and $\max \{c_{p}\}$ .

A key difference between Tw110-IA and Tw110-LA lies in the 3-D disturbance evolution during the breakdown. In the former, following the nonlinear amplification of 3-D modes, $\Lambda$ -like vortices (see figure 16) form gradually inside each streamwise-convected 2-D billow roll-up. Unlike the abrupt, staggered onset in case Tw110-LA, these vortices emerge naturally in an aligned formation without the external forcing of modes $(0,1)$ or $(1,\pm 1)$ ; cf. Rist & Fasel (Reference Rist and Fasel1995). This development inside the single billow roll-up is highlighted in figure 17 between $t/T_0=0.5$ and $t/T_0=1.0$ ( $7.42 \leq Re_{x}/10^5 \leq 7.79$ ). As the billow roll-up travels downstream, the separation zone at the spanwise ‘peak’ location weakens (see inset of figure 17 b) due to the recirculation of high-velocity fluid toward the wall by the $\Lambda$ -vortex legs. The rapid growth of the 3-D modes around $\textit{Re}_{x}/10^5 \approx 7.6$ intensifies the ‘upper’ layers at the ‘peak’ locations ( $z=0.5\lambda _{z},1.5\lambda _{z},\ldots$ ), which emerge from the billow crest. Meanwhile, the near-wall separation zones generate small, oppositely signed longitudinal vortices (inset of figure 17 c). At $t/T_0=1$ , the spanwise ‘peak’ features an ‘upper’ high-shear $\Lambda$ -shaped layer (see $|\omega _{z}|=\mathrm{const.}$ contours in figure 17 e), which enhances the spanwise flow modulation, while its legs are pushed closer to the wall, giving rise to regions of high $\omega _{z}$ and additional secondary vortical structures; cf. Rist & Fasel (Reference Rist and Fasel1995). This ‘upper’ layer is analogous to that of case Tw110-LA in figure 15. Upstream, the ‘inverted lower’ layer develops inside the boundary layer in a manner analogous to that observed in case Tw110-LA, transforming the initial ‘peak’–‘valley’ splitting into a ‘peak’–‘co-peak’ breakdown scenario. This high-shear layer originates at the trough of the billow roll-up and generates strong secondary low-density co-rotating vortices at the intermediate spanwise ‘co-peak’ locations (figure 17 g,h), which travel faster than the flow in the near-wall region of the spanwise ‘peak’. These longitudinal structures, visible in figure 16 around $\textit{Re}_{x}/10^5 \approx 7.8$ , break up earlier than the far-wall hairpin-like vortices, and contribute to the formation of turbulent structures also at the spanwise ‘co-peak’ locations.

Figure 18. Contours of the time- and spanwise-averaged (a,b) streamwise velocity $\langle u\rangle _{t,z}$ and (c,d) density $\langle \rho \rangle _{t,z}$ for (a,c) case Tw095-LA and (b,d) case Tw110-LA. In (a,b), the displacement thickness $\delta _{1}$ and momentum thickness $\theta$ are indicated by white dotted and dashed lines, respectively. Insets in (a,b) show selected streamwise velocity profiles normalised by the corresponding $u(\delta _{99})=0.99$ . The insets in (c,d) depict selected density profiles normalised by the corresponding $\rho (\delta _{99})$ . In (b) and (d), the Widom line $y=y_{\textit{WL}}$ lies within the green-shaded region, i.e. between $95\,\% \max \{c_{p}\}$ and $\max \{c_{p}\}$ .

5.3. Integral quantities

To illustrate the transition process, figure 18 presents isocontours of time- and spanwise-averaged streamwise velocity and density. The velocity profiles in figure 18(a,b) illustrate the transition to turbulence for cases Tw095-LA and Tw110-LA, respectively. In the latter, the near-wall separation zones are absent, as they are convected downstream and disappear on average. In the turbulent regime, the velocity and density profiles of both cases become qualitatively similar, regardless of whether the fluid is weakly or strongly non-ideal. In the strongly non-ideal case Tw110-LA, after transition to turbulence, the Widom line (coloured in green in figure 18 d) shifts significantly closer to the wall, considerably reducing the height of the vapour-like region compared with the laminar regime. The implications of this shift on the turbulent boundary layer are discussed in § 6.

To quantitatively characterise the transitional boundary layers following the qualitative analysis in § 5.2, the streamwise evolution of the time- and spanwise-averaged skin-friction coefficient $C_{\!{f}}$ and the Stanton number $St$ is examined. These quantities are defined in dimensionless form as

(5.1a,b) \begin{align} C_{\!{f}}=\dfrac {\overline {\tau }^*_{w}}{0.5\rho ^*_\infty u^{*2}_\infty }=\dfrac {2\overline {\tau }_{w}}{\textit{Re}}, \quad St=\dfrac {\overline {q}^*_{w}}{\rho ^*_\infty u^*_\infty \big(h^*_{aw}-h^*_{w}\big)}=\dfrac {\overline {q}_{w}}{\textit{Re}\, \textit{Pr}_\infty\, \textit{Ec}_\infty\, (h_{aw}-h_{w})},\\[-18pt]\nonumber \end{align}

where the non-dimensional characteristic parameters are defined in (2.3a – c ). The adiabatic wall enthalpy $h^*_{aw}=h^*_\infty +ru^{*2}_\infty /2$ , where the recovery factor $r=(h_{aw}-h_\infty )/(h_0-h_\infty )$ ( $h_0$ is the total enthalpy) is given as $\textit{Pr}^{1/3}_\infty =1$ (White Reference White2006), can be expressed in reduced form as $h_{r,aw}=h_{r,\infty }+rM^2_{\infty } a^2_{r,\infty }/2$ , where the free-stream reduced enthalpy is defined as $h_{r,\infty }=e_{r,\infty }+p_{r,\infty }/\rho _{r,\infty }$ . It is important to note that the assumption $r \approx 1$ has been verified in the laminar boundary layer under transcritical conditions (not shown). Figure 19(a,b) show the streamwise evolutions of $C_{\!{f}}$ and $St$ , respectively, for all cases in table 2. The theoretical laminar curves $C_{\textit{f},\textit{lam.}}$ and $St_{\textit{lam.}}$ for a non-ideal fluid are derived from (5.1a , b ) using the self-similar Lees–Dorodnitsyn variables as

(5.2a,b) \begin{align} C_{\textit{f},\textit{lam.}}=\dfrac {\sqrt {2}C_{w}}{\textit{Re}}\dfrac {\partial u}{\partial \eta }\bigg |_{w}, \quad St_{\textit{lam.}}=\dfrac {1}{\sqrt {2}\, \textit{Re}\, \textit{Ec}_\infty\, (h_{aw}-h_{w})} \dfrac {C_{w} \, c_{p,r,w}}{\textit{Pr}_{w} \, c_{p,r,\infty }}\dfrac {\partial T}{\partial \eta }\bigg |_{w} , \end{align}

where $\partial u/\partial \eta |_{w}$ and $\partial T/\partial \eta |_{w}$ are the wall-normal gradients of the laminar streamwise velocity and temperature in figure 3, respectively, and $C_{w}=\rho _{w}\mu _{w}$ is the Chapman–Rubesin parameter at the wall. In contrast, turbulent correlations for non-ideal wall-bounded flows are not yet available. An estimation is presented in § 6, where we introduce an accurate estimation solver.

Figure 19. Time- and spanwise-averaged (a) skin-friction coefficient $C_{\!{f}}$ , (b) Stanton number $St$ , and (c) shape factor $H_{12}$ , as functions of $\textit{Re}_{x}$ . Solid lines denote DNS results, while circle symbols represent the self-similar laminar correlations from (5.2a , b ) with initial conditions as in § 3.1. In (a,b), the theoretical incompressible skin-friction coefficient and the theoretical incompressible Stanton number using the Reynolds analogy $St=0.5C_{\!{f}}\textit{Pr}_\infty ^{-2/3}$ are represented by square and diamond black symbols, respectively (White Reference White2006). Note that for TadIG, $St=0$ due to adiabatic wall conditions. In (c), the theoretical incompressible shape factor is indicated with triangle black symbols (White Reference White2006).

The skin-friction and Stanton number curves initially follow the laminar trend in all cases. In Tw095-LA, both $\partial u/\partial \eta |_{w}$ and $C_{w}$ in (5.2a , b ) resemble the ideal-gas case, resulting in a $C_{\textit{f},\textit{lam.}}$ distribution that closely matches the ideal-gas case. In contrast, for both Tw110-LA and Tw110-IA, despite a fuller laminar velocity profile, $C_{w}$ is significantly below unity due to density and viscosity stratification. Thus $C^{{Tw110}}_{\textit{f},\textit{lam.}}\lt C^{{Tw095}}_{\textit{f},\textit{lam.}}$ . A similar trend holds for $St$ in figure 19(b). In Tw095-LA, $St$ closely matches the incompressible limit, i.e. $St=0.332\,\textit{Pr}^{-2/3}/Re$ . In Tw110-LA and Tw110-IA, the enthalpy difference $(h_{aw}-h_{w})$ is significantly larger, yielding $St^{{Tw110}}_{\textit{lam.}}\lt St^{{Tw095}}_{\textit{lam.}}$ .

In Tw095-LA, transition begins as mode $(1/2,1)$ overtakes $(1,0)$ and its amplitude exceeds $1\,\%$ . Between $\textit{Re}_{x}/10^5 \approx 5.3$ and $6.0$ , coinciding with the saturation of $(1/2,1)$ and other 3-D modes, both $C_{\!{f}}$ and $St$ level off as elongated $\Lambda$ -vortices become the predominant flow structures. At $\textit{Re}_{x}/10^5 \approx 6.0$ , a sharp rise in $C_{\!{f}}$ and $St$ is observed, as the amplification of modes $(1/2,1)$ , $(0,2)$ and higher 3-D modes leads to strong MFD and the roll-up of the $\Lambda$ -vortices into hairpin vortices that evolve into ring-like structures. Consequently, both $C_{\!{f}}$ and $St$ overshoot, consistent with the ideal-gas case TadIG.

In Tw110-LA and Tw110-IA, transition is delayed. The $C_{\!{f}}$ and $St$ curves begin to deviate from the laminar predictions around $\textit{Re}_{x}/10^5\approx 5.7$ , following the subharmonic resonance of modes $(1,0)$ and $(2,0)$ , and the appearance of near-wall separation zones. However, these separation zones are convected downstream such that the average $C_{\!{f}}$ remains positive, unlike a classic laminar separation bubble, which exhibits negative $C_{\!{f}}$ due to a larger, quasi-steady separation region (Alam & Sandham Reference Alam and Sandham2000). The $C_{\!{f}}$ and $St$ curves for cases Tw110-LA and Tw110-IA remain identical up to $\textit{Re}_{x}/10^5\approx 6.9$ , beyond which Tw110-LA exhibits a rapid rise in $C_{\!{f}}$ due to the abrupt amplification of 3-D disturbances. In contrast, Tw110-IA shows a sharp increase in both $C_{\!{f}}$ and $St$ only at $\textit{Re}_{x}/10^5\approx 7.65$ , followed by a kink in their evolutions around $\textit{Re}_{x}/10^5\approx 8.2$ , associated with the strong saturation of the steady modes $(0,1)$ and $(0,2)$ . Qualitatively, the minor decrease in $C_{\!{f}}$ and $St$ corresponds to the region between the break-up of two spanwise rollers, where no prominent vortical structures are observed (see figure 12(c) between $\textit{Re}_{x}/10^5\approx 8.0$ and $8.1$ ). Further downstream, the distributions of both transcritical cases gradually level off without a pronounced overshoot. Although mode $(0,0)$ , representing the wall-normal amplitude maximum, is approximately as large as in the subcritical case, the longitudinal vortex mode $(0,2)$ is significantly smaller, and no distinct overshoot appears. Accordingly, the reduction in $C_{\!{f}}$ and $St$ is more pronounced in the turbulent regime than in the laminar region, when compared to the subcritical case. The contribution of $\overline {\mu ^{\prime }\,\partial u^{\prime }/\partial y|_{w}}$ to the wall shear stress $\overline {\tau }_{w}$ is found to be negligible. Despite different initial forcings, the $C_{\!{f}}$ and $St$ curves of Tw110-LA and Tw110-IA converge as expected towards the same turbulent values.

The streamwise development of the shape factor $H_{12}=\delta _{1}/\theta$ , where $\delta _{1}=\int _0^{\infty } (1 - \rho u) \, \mathrm{d}y$ is the displacement thickness, and $\theta =\int _0^{\infty } (\rho u)(1-u) \, \mathrm{d}y$ is the momentum thickness, is shown in figure 19(c). In the transcritical cases, the increase of $H_{12}$ ( $H_{12} \approx 3.76$ ) in the laminar boundary layer is driven by a significant rise in $\delta _1$ , caused by the reduction in streamwise momentum $\rho u$ in the near-wall vapour-like regime, and a decrease in $\theta$ resulting from the fuller velocity profile compared to the subcritical case. Before $H_{12}$ drops to approximately $2.0$ , primarily due to the sharp rise in the turbulent $\theta$ value (see figure 18 b), it exhibits minor oscillations between $\textit{Re}_{x}/10^5 \approx 5.0$ and $6.0$ , attributed to the localised increase and decrease in $\delta _1$ over the near-wall separation zones.

6.1. Mean velocity and enthalpy–velocity transformations

We consider the transformation proposed by Patel et al. (Reference Patel, Boersma and Pecnik2016) and Trettel & Larsson (Reference Trettel and Larsson2016), which extends the van Driest (Reference van Driest1951) velocity transformation (subscript ${vD}$ ) defined as $\overline {u}^+_{vD} = \int ^{\overline {u}^+}_0 \sqrt {\overline {\rho }/\overline {\rho }_{w}} \, \mathrm{d}\overline {u}^+$ , by accounting for spatial variations in the viscous length scale $\delta _{v}^\star$ in the near-wall region. It is formulated as

(6.1) \begin{equation} \overline {u}^{\kern0.75pt\star} = \int ^{\overline {u}^+}_0 \left ( 1 - \dfrac {y}{\delta _{v}^\star } \dfrac {\mathrm{d}\delta _{v}^\star }{\mathrm{d}y} \right ) \underbrace {\frac {u_\tau }{u^\star _\tau }\, \mathrm{d}\overline {u}^+}_{\mathrm{d}\overline {u}^+_{vD}}, \end{equation}

where the superscripts ${+}$ and ${\star }$ indicate non-dimensionalisation by the viscous length scale $\delta _{v}$ and friction velocity $u_\tau$ (as defined in Appendix A), and by the semi-local viscous length scale $\delta ^{\star }_{v}=\overline {\mu }/(\overline {\rho } {u^{\star }_{\tau }})$ and semi-local friction velocity ${u^{\star }_{\tau }}=\sqrt {{\overline \tau _{w}}/\overline {\rho }}$ , respectively. The Reynolds (time and spanwise) average of a given variable $\chi$ is expressed as $\overline {\chi }=\chi -\chi ^{\prime }$ , with $\chi ^{\prime }$ denoting the corresponding fluctuation. It is worth noting that the semi-local Mach number $M_\tau =u_\tau /\overline a_{w}$ remains much smaller than unity across the entire domain (with maximum value $0.025$ ). Therefore, the transformation recently proposed by Hasan et al. (Reference Hasan, Larsson, Pirozzoli and Pecnik2023), which extends (6.1) to account for intrinsic compressibility effects, would yield results equivalent to $\overline {u}^{\kern0.75pt\star}$ .

Figure 20. Wall-normal profiles of the transformed streamwise velocity using (a) van Driest (Reference van Driest1951) and (b) Patel et al. (Reference Patel, Boersma and Pecnik2016). Cases Tw095-LA and Tw110-LA are shown in orange at $\textit{Re}_{\theta }=1387$ ( $\textit{Re}_{x}=1.0 \times 10^6$ ) and in red at $\textit{Re}_{\theta }=881$ ( $\textit{Re}_{x}=1.06 \times 10^6$ ), respectively. Grey dash-dotted lines denote the linear and logarithmic laws (a) $(1/\kappa )\log y^+ + C$ , (b) $(1/\kappa )\log y^{\star } + C$ , with $\kappa =0.41$ and $C=5.2$ . Blue shows the ideal-gas case TadIG at $\textit{Re}_{\theta }=1190$ ( $\textit{Re}_{x}=0.8 \times 10^6$ ).

Figure 20(a,b) show the transformed Reynolds-averaged mean velocity profiles for $\overline {u}^+_{vD}$ and $\overline u^\star$ on logarithmic scales of $y^+$ and $y^\star$ , respectively. The turbulent profiles are extracted near the end of the computational domains to ensure well-developed turbulent boundary layers: for Tw095-LA at $\textit{Re}_{\theta }=1387$ ( $\textit{Re}_{x}=1.0 \times 10^6$ ), for Tw110-LA at $\textit{Re}_{\theta }=881$ ( $\textit{Re}_{x}=1.06 \times 10^6$ ), and for TadIG at $\textit{Re}_{\theta }=1190$ ( $\textit{Re}_{x}=0.8 \times 10^6$ ). As expected, both velocity transformations of the subcritical case Tw095-LA show a good collapse with the velocity of the constant-property turbulent boundary layer (case TadIG), due to the low variation of thermophysical properties across the boundary layer. In contrast, the van Driest (Reference van Driest1951) scaling shows a large offset in the logarithmic region for the transcritical case Tw110-LA. For this case, the large near-wall variation of the viscous length scale must be taken into account. With this correction, (6.1) clearly improves the collapse with the constant-property boundary layer (see figure 20 b).

Contrary to the findings of Kawai (Reference Kawai2019), this collapse is achieved for the following reasons: in case Tw110-LA, the wall-normal velocity in the boundary layer remains relatively small, with $\overline {v}/u_\infty \approx 0.2 \,\%$ . Furthermore, the density and viscosity fluctuations, which peak in the buffer layer, are moderate, of the order of $\sqrt {\overline {\rho ^{\prime }\rho ^{\prime }}}/\overline {\rho }$ and $\sqrt {\overline {\mu ^{\prime }\mu ^{\prime }}}/\overline {\mu } \approx 0.3$ , respectively. By contrast, in the transcritical case investigated by Kawai (Reference Kawai2019) at a higher reduced pressure ( $p_{r} \approx 1.28$ ) but significantly larger wall-to-free-stream temperature ratios ( $T^*_{w}/T^*_\infty \approx 4$ – $8$ ), the wall-normal velocity reached approximately $1.5\,\%$ , while the density and viscosity fluctuations reached $\sqrt {\overline {\rho ^{\prime }\rho ^{\prime }}}/\overline {\rho } \approx 1.0$ and $\sqrt {\overline {\mu ^{\prime }\mu ^{\prime }}}/\overline {\mu } \approx 0.4$ , particularly in the logarithmic layer. These strong fluctuations led to large near-wall convective flux and an overshoot of the Reynolds shear stress surpassing unity in the logarithmic region, ultimately causing the failure of the velocity transformation. In our case, however, the total stress balance, on which the Patel et al. (Reference Patel, Boersma and Pecnik2016) transformation relies, is still mainly governed by the viscous and turbulent stresses under transcritical conditions. Overall, the transformed velocity profiles indicate that the flow has transitioned into the fully turbulent regime by the end of the computational domain.

Figure 21. Case Tw095-LA ( $\textit{Re}_{\theta }=1387$ ) in orange and Tw110-LA ( $\textit{Re}_{\theta }=881$ ) in red: (a) mixed Prandtl number $\textit{Pr}_{m}$ , (b) turbulent Prandtl number $\textit{Pr}_{t}$ , and (c) mean molecular Prandtl number $\overline {\textit{Pr}}$ and $\overline {c}_{p}/c_{p,\infty }$ (red dash-dotted line) for case Tw110-LA.

Next, we focus on estimating the mean thermodynamic properties to develop a predictive model for $C_{\!{f}}$ and $St$ , which remains unknown for turbulent boundary layers with highly non-ideal fluids. One possible approach is to directly integrate the total heat flux from the wall to the free stream to reconstruct the temperature field. However, this method is cumbersome, as it requires an iterative procedure on the heat flux to match the prescribed wall-to-free-stream temperature ratio, combined with a turbulence model that itself depends on the velocity field as well as the density- and temperature-dependent thermodynamic properties. As a more practical alternative, we examine the mean temperature–velocity correlations, which offer a simpler route. When non-ideal-gas effects become significant, enthalpy-based relations must be employed, as enthalpy is no longer a linear function of $T$ but instead depends on another thermodynamic property. Furthermore, classical enthalpy-based relations such as Crocco–Busemann (Smits & Dussauge Reference Smits and Dussauge2006), or the relations of Walz (Reference Walz1969) and Duan & Martín (Reference Duan and Martín2011), assume a constant mixed Prandtl number $\textit{Pr}_{m}$ ( $\textit{Pr}_{m}=1$ in Crocco–Busemann), defined as

(6.2) \begin{align} \textit{Pr}_{m}=\dfrac {(\overline {\mu }+\mu _t)\overline {c}_{p}}{\overline {\kappa }+\kappa _t}, \quad \text{with}\ \mu _{t}=\dfrac {-\overline {\rho }\,\widetilde {u^{\prime \prime }v^{\prime \prime }}}{\partial \widetilde {u}/\partial y}\ \text{and}\ \kappa _{t} = \dfrac {-\overline {c}_{p} \overline {\rho }\, \widetilde {v^{\prime \prime }T^{\prime \prime }} }{\partial \widetilde {T}/\partial y},\end{align}

where $\mu _{t}$ and $\kappa _{t}$ are the eddy viscosity and eddy conductivity, respectively. Note that $\widetilde {\chi }=\overline {\rho \chi }/\overline {\rho }$ denotes the Favre average, with $\chi ^{\prime \prime }$ as the Favre fluctuation.

The mixed Prandtl number $\textit{Pr}_{m}$ , turbulent Prandtl number $\textit{Pr}_{t}=\mu _{t}\overline {c}_{p}/\kappa _{t}$ , and molecular Prandtl number $\overline {Pr}$ are plotted in figure 21(a–c), respectively, for the subcritical and transcritical cases. In case Tw095-LA, the mixed Prandtl number remains approximately constant and close to unity across the boundary layer ( $\textit{Pr}_{m,w} \approx 1.06$ ), despite the turbulent Prandtl number decreasing to approximately $0.7$ in the outer region (see figure 21 b). The evolution of $\textit{Pr}_{m}$ confirms the assumption underlying the aforementioned enthalpy-based relations with a constant Prandtl number. Under transcritical conditions, however, $\textit{Pr}_{m}$ deviates from unity in the logarithmic region, reaching $1.84$ in the buffer layer before dropping below unity in the viscous sub-layer ( $\textit{Pr}_{m,w} \approx 0.51$ ). Observing the mean molecular Prandtl number $\overline {Pr}$ in figure 21(c), its behaviour underlines that the evolution of the $\textit{Pr}_{m}$ is dominated by molecular diffusion. The strong increase and subsequent decrease of $\textit{Pr}_{m}$ towards the wall are proportional to the variation of $\overline {c}_{p}$ , indicating that the Widom line is predominantly located in the buffer layer in the turbulent regime. In contrast, the turbulent Prandtl number moderately decreases from unity in the near-wall region to approximately $0.6$ .

To address the variable $\textit{Pr}_{m}$ in the highly non-ideal case Tw110-LA, we consider the variable-Prandtl theory of van Driest (Reference van Driest1955). For a turbulent boundary layer in mean steady state over a ZPG flat plate, it holds that

(6.3) \begin{equation} \left (\frac {\overline {i}^{\prime }}{\textit{Pr}_{m}}\right )^{\prime }+\left (1-\textit{Pr}_{m}\right ) \frac {\tau ^{\prime }}{\tau }\left (\frac {\overline {i}^{\prime }}{\textit{Pr}_{m}}\right )=-\frac {u_{\infty }^{2}}{h_{\infty }}, \end{equation}

where $\overline {i}=\overline {h}/h_\infty$ is the mean enthalpy ratio, $\tau$ is the total shear stress, and $({\cdot })^{\prime }$ indicates differentiation with respect to $\overline {u}_{r}=\overline {u}/u_\infty$ . Integrating (6.3) yields

(6.4) \begin{equation} \overline {i} = \dfrac {h_{w}}{h_\infty } - \dfrac {\mathcal{S}}{\mathcal{S}_\infty } \left ( \dfrac {h_{w}}{h_\infty } - 1 \right ) + \dfrac {u_\infty ^2}{h_\infty } \left [ \dfrac {\mathcal{S}}{\mathcal{S}_\infty } \mathcal{R}_\infty - \mathcal{R} \right ]\!, \end{equation}

where the functions $\mathcal{S}$ and $\mathcal{R}$ are defined as

(6.5a) \begin{align}\mathcal{S}=\int _{0}^{\overline {u}_{r}} \textit{Pr}_{m} \exp \left [-\int _{\tau _{w}}^{\tau } \frac {1-\textit{Pr}_{m}}{\tau }\, \mathrm{d} \tau \right ] \mathrm{d} \overline {u}_{r},\end{align}

(6.5b) \begin{align}\mathcal{R}=\int _{0}^{\overline {u}_{r}} \textit{Pr}_{m} \exp \left [-\int _{\tau _{w}}^{\tau } \frac {1-\textit{Pr}_{m}}{\tau }\, \mathrm{d} \tau \right ]\left \{\int _{0}^{\overline {u}_{r}} \exp \left (\int _{\tau _{w}}^{\tau } \frac {1-\textit{Pr}_{m}}{\tau }\, \mathrm{d} \tau \right ) \mathrm{d} \overline {u}_{r}\right \} \mathrm{d} \overline {u}_{r}.\end{align}

Note that for $\textit{Pr}_{m}=1$ and $\textit{Pr}_{m}=\mathrm{const.}$ , the classical relations of Crocco–Busemann (Smits & Dussauge Reference Smits and Dussauge2006) and Walz (Reference Walz1969) are recovered, respectively.

Figure 22(a,b) compare (6.4) against DNS data for the heated case Tw095-LA ( $T^*_{w}/T^*_\infty =1.056$ ) and Tw110-LA ( $T^*_{w}/T^*_\infty =1.222$ ), respectively. For the subcritical case, van Driest’s relation accurately predicts $\overline {h}/h_\infty$ and aligns with Walz’s relation as $\textit{Pr}_{m} \approx \mathrm{const.}$ However, for the highly non-ideal fluid with strongly varying mixed Prandtl number, as seen in figure 21(a), Walz’s relation becomes particularly inaccurate in the inner layer. From $\overline {u}/u_\infty \approx 0.6$ , or $y^\star \approx 20$ , the enthalpy undergoes a significant increase, which only van Driest’s relation is able to capture in good agreement with the DNS data.

Figure 22. Reynolds-averaged mean enthalpy from the variable-Prandtl theory of van Driest (Reference van Driest1955) in black for (a) Tw095-LA at $\textit{Re}_{\theta }=1387$ (DNS profile in orange) and (b) Tw110-LA at $\textit{Re}_{\theta }=881$ (DNS profile in red). In grey is shown the relation of Walz (Reference Walz1969) as $\overline {h}/h_\infty =h_{w}/h_\infty +(h_{aw}-h_{w})(\overline {u}/u_\infty )/h_\infty -ru^2_\infty (\overline {u}/u_\infty )^2/(2h_\infty )$ .

6.2. Estimating mean profiles and fluxes

We now follow the approach of Hasan et al. (Reference Hasan, Larsson, Pirozzoli and Pecnik2024) to predict the drag and heat transfer, and compare the results to the available DNS data (see figure 19). In this approach, the mean shear is integrated from the wall to the free stream to obtain the streamwise velocity, using a combination of inner- and outer-layer modelling approximations. The inner-layer eddy viscosity is modelled using the Johnson–King mixing-length formulation (Johnson & King Reference Johnson and King1985), accounting for variations in viscous length and velocity scales, while the outer layer is described by Coles’ law of the wake (Coles Reference Coles1956), modified to include mean density variations. The corresponding thermodynamic properties, e.g. $\rho$ , $T$ , $\mu$ , $\kappa$ and $\textit{Pr}$ , are computed using the VdW EoS, under the assumption of constant thermodynamic pressure, and the JST model, with both expressed as functions of the enthalpy obtained via van Driest’s variable-Prandtl theory. Note that two additional assumptions are made: (i) $\textit{Pr}_{m} \approx \overline {Pr}$ (compare figure 21 a,c), and (ii) the shear distribution follows van Driest (Reference van Driest1955), i.e. $\tau /\tau _{w} \approx 1 - \exp [-\kappa /\sqrt {C_{\!{f}}/2}(1-\overline {u}_{r})]$ , with $\kappa =0.41$ as the von Kármán constant. For further details on the analytical solver, refer to the source code available on GitHub (https://github.com/pcboldini/DragAndHeatTransferEstimation_NonIdealFluids).

Figure 23 presents the predicted temperature, density, viscosity and velocity profiles from the analytical solver (‘Estimation’) compared with the DNS profiles for cases Tw095-LA and Tw110-LA. In the subcritical case, the estimated profiles are in very good agreement with the DNS solution. In the transcritical case, the temperature profile in figure 23(b) shows minimal deviations in the buffer layer but successfully approximates the wall heat flux $\overline {q}_{w}$ . The estimated profiles for density (figure 23 d) and viscosity (figure 23 f) differ significantly only in the buffer layer, where the largest discrepancies in $\overline {h}$ are observed. Note that Jensen’s inequality, i.e. $\overline {\chi } \neq \chi (\overline {h},\overline {p})$ , applies to these thermophysical quantities (Nemati et al. Reference Nemati, Patel, Boersma and Pecnik2015), due to their sharp curvature around the pseudo-critical point, i.e. $\max \{c_{p}\}$ , and their large fluctuations present specifically in the buffer layer. Conversely, the estimated velocity profile agrees very well with the DNS solution, as it results from integration across both inner and outer layers.

Figure 23. Estimated mean (a,b) temperature $\overline {T}/T_\infty$ , (c,d) density $\overline {\rho }/\rho _\infty$ , (e, f) viscosity $\overline {\mu }/\mu _\infty$ , and (g,h) streamwise velocity $\overline {u}^{{\kern0.5pt}+}$ profiles (dashed grey lines) compared to DNS results (solid lines, case Tw095-LA ( $\textit{Re}_{\theta }=1387$ ) in orange, and case Tw110-LA ( $\textit{Re}_{\theta }=881$ ) in red).

Figure 24. Reynolds analogy factor $s=2\,St/C_{\!{f}}$ as a function of the momentum-thickness Reynolds number $\textit{Re}_{\theta }$ : case Tw095-LA (orange) and Tw110-LA (red). Solid lines correspond to the DNS results, while dashed lines represent the turbulent Reynolds analogy factor $s=\mathcal{S}_\infty$ according to the variable-Prandtl theory of van Driest (Reference van Driest1955). The Reynolds analogy is indicated by a black dotted line.

Figure 25. Case Tw095-LA (orange) and Tw110-LA (red): (a) skin-friction coefficient $C_{\!{f}}$ and (b) Stanton number $St$ as functions of the momentum-thickness Reynolds number $\textit{Re}_{\theta }$ . In (a,b), solid lines correspond to the DNS results, while dashed lines denote the analytical estimations. In (b), dotted lines represent $St=C^{\prime }_{f}/2$ ( $\textit{Pr}_\infty =1$ ) according to the Reynolds analogy, where $C^{\prime }_{f}$ is the analytical prediction from (a).

Before addressing the prediction of the skin-friction coefficient and Stanton number, we first examine the Reynolds analogy, which relates skin friction and heat transfer. Using the Reynolds analogy factor $s=2St/C_{\!{f}}=\overline {q}^*_{w}u^*_\infty /(\overline {\tau }^*_{w}(h^*_{aw}-h^*_{w}))$ , we investigate how non-ideal-gas effects may disrupt the classical coupling between momentum and thermal transport. The impact of the Reynolds analogy is shown in figure 24. In the laminar boundary layer, strong property variations in the transcritical regime reduce the Reynolds analogy factor below unity ( $s \approx 0.89$ ). In the transitional region, the value of $s$ rises above unity in both thermodynamic regimes due to the gradual development of turbulent structures. Interestingly, $s$ appears largely insensitive to the degree of fluid non-ideality in the turbulent boundary layer. Here, the Reynolds analogy holds in both thermodynamic regimes, with $s \approx 1$ . The behaviour aligns with trends observed in ideal-gas turbulent boundary layers, where $s$ was found to be unaffected by variations in wall temperature and Mach number (Wenzel et al. Reference Wenzel, Gibis, Kloker and Rist2021). Nevertheless, a more comprehensive investigation over a broader Reynolds number range would be needed to extend the applicability of these results to other non-ideal fluid flows. Note that the turbulent Reynolds analogy factor, calculated according to the variable-Prandtl theory of van Driest (Reference van Driest1955) in (6.3), with $s=\mathcal{S}_{\infty }$ , follows the DNS results only in the subcritical regime. In the transcritical regime, larger deviations in the estimation of $\tau$ are observed.

In figure 25(a,b), the DNS and estimated values of the skin-friction coefficient $C_{\!{f}}$ and Stanton number $St$ are plotted as functions of $\textit{Re}_{\theta }$ . The predictions of $C_{\!{f}}$ show good agreement with the DNS curves beyond the respective overshoots. Moreover, the less pronounced overshoot in $C_{\!{f}}$ for case Tw110-LA is confirmed. In contrast, the $St$ prediction reveals good accuracy only for the subcritical case Tw095-LA. For the highly non-ideal case Tw110-LA, the $St$ number is over-predicted by ${\sim} 30\,\%$ at $\textit{Re}_{\theta }=881$ due to a slightly higher estimated wall heat flux in figure 23(b). Note that the error in the Stanton number prediction using the relation of Walz (Reference Walz1969) reaches ${\sim} 200\,\%$ at $\textit{Re}_{\theta }=881$ . Interestingly, when applying the Reynolds analogy $St=C^{\prime }_{f}/2$ for $\textit{Pr}_\infty =1$ , verified in figure 24, and using $C^{\prime }_{f}$ from the analytical prediction in figure 25(a), the estimated $St$ curve closely agrees with the DNS results in the turbulent regime. In conclusion, this behaviour suggests that even under the examined transcritical conditions with variable mixed Prandtl number, the Reynolds analogy may remain a robust tool for predicting the turbulent heat transfer in ZPG flat-plate boundary layers.