About the role of the Hanratty correction in the linear response of a turbulent flow bounded by a wavy wall

Abstract

Scallop patterns forming on erodible surfaces were studied historically using a linear analysis of the inner region of a turbulent boundary layer growing on a corrugated wall. Experimental observations show a phase shift between the shear stress at the wall and the wall oscillation that depends on the wavenumber. An ad hoc correction applied to the turbulent closure and due to Hanratty et al. (Thorsness et al., Chem. Engng Sci., vol. 33, issue 5, 1978, pp. 579–592; Abrams & Hanratty, J. Fluid Mech., vol. 151, issue 1, 1985, p. 443; Frederick & Hanratty, Exp. Fluids, vol. 6, issue 7, 1988, pp. 477–486) was systematically used to recover the reference experimental results. In this study, Reynolds-averaged Navier–Stokes (RANS) and direct numerical simulations (DNS) were performed and revealed the role of the Boussinesq assumption in the results obtained. We show that the Hanratty correction acts as a palliative to the misrepresentation of Reynolds stresses due to the use of the Boussinesq hypothesis. The RANS calculations based on a turbulence model using a second-order moment closure recovered the expected results obtained in the reference DNS calculations, in particular with respect to wall heat transfer. The analysis of these results highlights the critical importance of the anisotropy of the diagonal Reynolds stresses on the prediction of wall transfer under these conditions and their implication in the occurrence of scalloping.

1. Introduction

Scallop patterns are found in a large variety of situations, characterising the interaction of a fluid and an erodible surface. They are observed on meteorites, and called regmaglypts (Lin & Qun 1987; Claudin & Ernstson 2004), in pipes (Blumberg & Curl 1974; Villien, Zheng & Lister 2001, 2005), in karst or ice caves (Anderson et al. 1998; Sundqvist, Seibert & Holmgren 2007; Pflitsch et al. 2017) or with dunes (Best 2005; Vinent et al. 2019) and sand ripples (Bagnold 1941; Charru, Andreotti & Claudin 2013). Many examples of these scallop patterns are listed by Claudin, Duran & Andreotti (2017). Thomas (1979) gathered several experimental results and provides evidence of a unique scaling of the wavelength of the scallops with the boundary-layer viscous length. Similar patterns are also observed on atmospheric re-entry vehicles. During the re-entry phase in hypersonic conditions, the windward face of a vehicle is exposed to severe heat fluxes due to the post-shock environment. Carbon-based thermal protection systems are commonly used to guarantee the integrity of the payload. The carbon oxidation and sublimation processes lead to the ablation of the heat shield, and under some conditions, scallops may be observed on vehicle nosetips. Few in-flight experiments are published (Canning, Tauber & Wilkins 1968; Larson & Mateer 1968), the most important reference being the TATER test (Hochrein & Wright 1976) for which scallops approximately $1$ to $4$ mm long and a depth $10$ times smaller were observed on the ablated surface as shown in figure 1. Several on-ground tests (Laganelli & Nestler 1969; Nestler 1971; Williams 1971; Baker 1972; White & Grabow 1973; Shimizu, Ferrell & Powars 1974; Reineke & Guillot 1995; Mikhatulin & Polezhaev 1996; Powars 2011), involving lower heat fluxes and using surrogate ablative materials such as camphor or Teflon, have confirmed the formation of scallops. The ablation process depends on the material and may imply decomposition or fusion. To study the formation of scallops on re-entry vehicles, we therefore rely on existing approaches for which several fundamental unresolved issues related to turbulence models still remain. The occurrence of these patterns on the surface of erodible walls were studied for many years by performing linear analyses (Benjamin 1959; Thorsness, Morrisroe & Hanratty 1978; Abrams & Hanratty 1985; Fourrière, Claudin & Andreotti 2010; Charru et al. 2013; Claudin et al. 2017). Classically, the surface regression rate is assumed to be small enough so that the associated characteristic time scale is very large compared with the mean flow characteristic time. The problem is then first reduced to the investigation of an incompressible turbulent boundary layer developing over a sinusoidally perturbed static surface. The linear forced response for this flow was first studied by Benjamin (1959) and consists of solving the Orr–Sommerfeld equation for a laminar flow. This problem was explored again by Hanratty and co-workers (Zilker, Cook & Hanratty 1977; Thorsness et al. 1978; Abrams & Hanratty 1985; Frederick & Hanratty 1988) providing a new insight into the linear response while introducing a slight modification to the Orr–Sommerfeld equation and considering turbulent flows. Thorsness et al. (1978) introduced a metric function to transpose the equations into the ’boundary-layer coordinate system’ before the linearisation. However, the base flow was moved together with the coordinate system and displaced to the new origin. This crucial modification was carefully analysed and discussed by Luchini & Charru (2019). In the present work, we take up the work of Fourrière et al. (2010) and Charru et al. (2013) to derive the linear problem. This is equivalent to the approach of Hanratty et al. and gives exactly the same results. The equations set and notations are recalled in Appendix A. Since the flow is supposed turbulent, a closure relation is used to model the contribution of Reynolds stresses in the stress tensor $\tau _{ij}$. In all the studies cited, the Boussinesq hypothesis (A3) is used together with a Prandtl mixing length model.

Figure 1. Scallops observed on in-flight and on-ground experiments representative of hypersonic re-entry vehicles. From left to right, nosetip pictures of the TATER experiment (Hochrein & Wright 1976), on-ground tests using camphor (Larson & Mateer 1968) and Teflon (Powars 2011) as surrogate material.

Simultaneously to their initial linear analysis, experimental works were conducted by Hanratty et al. (Zilker et al. 1977; Frederick & Hanratty 1988) providing essential data to validate the results of the linear analysis. A series of measurements in a turbulent channel flow equipped with a wavy wall highlighted a modulation of the wall shear stress phase with respect to the wall deformation in a specific wavelength range. The existence of a phase shift between the wall shear stress and the wavy wall can be explained by the momentum budget (Charru & Hinch 2000). For laminar flows, or simply as long as the perturbations are in the viscous sublayer, the pressure gradient induced by the wall waviness is responsible for the phase shift. For turbulent flows, other contributions may come into play, notably the diffusion term related to the difference in stresses $\tau _{xx}-\tau _{zz}$. When comparing the experimental observations and the linear analysis, Hanratty et al. noticed the failure of the mixing length model. Interestingly, by introducing a dependence of the mixing length to a relaxed pressure gradient, denoted $\mathcal {C}$ hereinafter, Hanratty and co-workers (Thorsness et al. 1978; Abrams & Hanratty 1985) were able to reproduce the behaviour of the wall shear stress phase. This correction to the mixing length was further reformulated by Charru et al. (2013) and Claudin et al. (2017) and used successfully.

To further elucidate how scalloping forms on erodible surfaces, the wall profile is made time dependent and is related to a wall flux involved in the transport mechanism controlling the wall recession. For sand ripple formation, the particle flux is used and is shown to be lagged behind the wall shear stress. The lag of the particle flux has a stabilising effect that balances the inertial destabilising effect of the shear stress. A thorough discussion is given in the review by Charru et al. (2013). For dissolution or melting problems, Claudin et al. (2017) considered a passive scalar transport equation, representing, for example, the concentration of a chemical species or the temperature, and the wall profile evolution is controlled by the wall normal flux of the scalar transported. The ablation problem on the nosetip of a re-entry vehicle can be apprehended in the same way but several issues must be addressed first, among which one is of key importance.

The correction $\mathcal {C}$ proposed by Hanratty is a heuristic model, made to recover measurement (Zilker et al. 1977) data for the wall shear stress from a mixing length approach. However, in order to close the passive scalar transport equation in the approach followed by Claudin et al. (2017), the turbulent scalar flux is related to the eddy viscosity based on the mixing length and including the correction $\mathcal {C}$. Assuming that $\mathcal {C}$ is a valid and sufficient correction for the turbulent scalar flux closure is far from being trivial and there are no existing data enabling us to validate this model. The choice of the closure is yet a determining factor for the assessment of the wall normal flux that controls the surface regression rate. To shed light on this point, we follow the approach presented by Claudin et al. (2017) for the transport of a passive scalar and in § 2 we study the forced response of the energy equation for an incompressible fluid. At first, a fixed corrugated surface is considered and a dedicated mixing length is proposed to model the turbulent scalar flux. The choice of the base flow is also discussed in this section to remove doubts about the relevance of the validation cases performed. In § 3, direct numerical simulations (DNS) are carried out to establish some validation points to complete the experimental data of Hanratty, notably concerning heat flux. Additionally, Reynolds-averaged Navier–Stokes (RANS) computations with first- and second-order moment closures are performed to discuss the influence of the turbulent closures on the momentum and energy equations. In the last § 4, through the analysis of the different types of results, we will discuss the achievements and some limitations of the Hanratty correction. Finally, a simple wall regression model, assuming scale separation between the ablation mechanism and the flow response, is presented to try to establish a link with the Thomas correlation. In particular, we highlight the key role played by the closure relation for the turbulent heat flux.

2. Linear forced response

2.1. Turbulent closure for the linearised momentum equations

We take up the work by Charru et al. (2013) to solve the linearised momentum equations, considering a steady and incompressible fluid flow, the corrugated surface being fixed in time at this stage. The notation and the system of equations are recalled in Appendix A.1. The study is restricted to the linear response of the flow to the wall undulation, i.e. the amplitude $\zeta _0$ of the wall deformation is small enough compared with the wavelength ${2{\rm \pi} }/{\alpha }$, with $\alpha$ the wavenumber. The nonlinear limit is $\alpha \zeta _0 \approx 0.1$ (Charru et al. 2013) whereas flow separations are expected for $\alpha \zeta _0>0.3$ (Zilker & Hanratty 1979). A dedicated code based on a collocation method (Canuto et al. 2006) using Chebyshev polynomials was developed to solve the linearised system. The Reynolds stresses are modelled with the help of the Boussinesq hypothesis (A3) and the eddy viscosity $\nu _t$ is deduced from a mixing length approach (A2). Thorsness et al. (1978) first proved that a correction is required to recover the experimental results (Zilker et al. 1977) showing large phase shifts of the wall shear stress with respect to the wall undulation in a specific wavenumber range, as illustrated in figure 2. The idea is to introduce a dependence to a relaxed pressure gradient for the van Driest number $A$ inspired by the work of Loyd, Moffat & Kays (1970) or similarly by that of Cebeci & Smith (1974). Since the mixing length $l$ (A2) depends on the non-dimensional variables, the wall normal coordinate $\eta$, the Reynolds number $\mathcal {R}$ based on the wavenumber $\alpha$ and the van Driest number $A$, the disturbed part of the mixing length $\hat {l}$ obtained after linearisation contains three distinct contributions

(2.1)\begin{equation} \hat{l} =-\kappa\left[1-\exp\left(-\frac{\mathcal{R}\eta}{A^0}\right) \left(1-\frac{\mathcal{R}\eta}{A^0}+\frac{\mathcal{R}\eta^2}{A^0} \left(\frac{\hat{\tau}_{xz}}{2}-\beta\mathcal{C}\right)\right)\right]. \end{equation}

The first one due to $\eta$ is the linearised effect of the geometrical deformation. The second reveals the influence of the wall shear stress disturbance $\hat {\tau }_{xz}$. Finally, the dependence to $\mathcal {C}$ is brought by the van Driest constant $A$ with $\beta$ the relative variation of $A$ due to the relaxed pressure gradient $\beta =({1}/{A^0})({\partial {A}}/{\partial {\mathcal {C}}})$. The parameter $A^0=26$ is the standard van Driest constant and $\beta =35$ is found to be the value that best fits the measurements (Frederick & Hanratty 1988; Charru et al. 2013). The dimensionless correction $\mathcal {C}$ is given by a differential equation that reads

(2.2)\begin{equation} \gamma \frac{\partial{\mathcal{C}}}{\partial{x}}=\frac{1}{u_{\tau}^2}\frac{\partial{}}{\partial{x}} \left(\tau_{xx}-\frac{p}{\rho}\right)-\mathcal{C} ,\end{equation}

where $p$ is the pressure, $\rho $ the density and $u_{\tau}$ the friction velocity. The constant $\gamma$ determines the length over which the relaxation operates with respect to the streamwise gradient of $\tau _{xx}-{p}/{\rho }$. Originally (Thorsness et al. 1978; Frederick & Hanratty 1988), $\mathcal {C}$ was only related to the pressure gradient, with similar results. The dimensionless quantity $\mathcal {C}$ does not correspond to the whole correction introduced in $\hat {l}$, but it will be called Hanratty's correction thereafter for brevity. When only the geometrical dependence of $\hat {l}$ is kept, and so the dependences on $\hat {\tau }_{xz}$ and $\mathcal {C}$ are dropped in (2.1), the turbulence can be seen as ‘frozen’ regarding the perturbations. This will be referred to as the frozen turbulence assumption in the following. More details on (2.1) and (2.2) can be found in the supplemental material of the review of Charru et al. (2013).

Figure 2. Phase of the wall shear stress in the transitional regime. Filled black circles denote Hanratty's experimental results. Solid lines are results of the linear analyses with the Hanratty correction $\mathcal {C}$ (blue) and under the frozen turbulence hypothesis (orange). Rectangles are results of RANS computations with the $k\unicode{x2013}\omega$ model (orange) and the elliptic blending Reynolds stress model (blue). Forced responses in channel flow are plotted with dashed blue lines for $\alpha \delta =2{\rm \pi}$ and $\beta =40$: – – –; $\alpha \delta ={\rm \pi}$ and $\beta =45$: $-.-.-$; $\alpha \delta ={{\rm \pi} }/{2}$ and $\beta =50$: $-..-..-$. The dashed orange line corresponds to the linear analysis where the Hanratty correction is off but the dependence on $\hat {\tau }_{xz}$ is conserved.

Experimental results and those of the linear analyses of the wall shear stress phase $\psi _{\tau }=\arg (\hat {\tau }_{xz})$ are plotted in figure 2 for wavenumbers in the transitional regime. Indeed, three regimes can be distinguished with respect to $\mathcal {R}$ and the penetration depth of the perturbation $\delta _i$. The first regime corresponds to small values of $\mathcal {R}$ $(\mathcal {R}<100)$, and, according to Charru & Hinch (2000), $\delta _i \propto \delta _{\nu }\mathcal {R}^{{1}/{3}}$, where $\delta _{\nu }$ is the viscous length ${\nu }/{u_{\tau }}$. The perturbation is confined in the viscous sublayer so that the turbulent closure plays no role in this regime. The third regime corresponds to the long wave approximation ($\mathcal {R}>10\,000$) for which the flow disturbances extend far beyond the viscous region where the Reynolds stresses cannot be neglected anymore. As recalled by Charru et al. (2013), velocity measurements confirm the linear increase in mixing length with wall distance in the logarithmic region. Therefore, in this regime, the results are little affected by the choice of turbulent closure as long as the linearity of the eddy viscosity with respect to the wall distance is recovered in the logarithmic region of the inner layer. The intermediate regime, i.e. $\mathcal {R}\in [100,10\,000]$, often called the transitional regime, is far more complex and more challenging. The linear analysis with the standard mixing length model, i.e. without the inclusion of correction $\mathcal {C}$, does not recover the trend measured, but the use of the Hanratty correction improves the results remarkably. The evolution of $\psi _{\tau }$ with $\alpha ^+=\mathcal {R}^{-1}$ from the laminar regime to the fully turbulent regime is then faithfully reproduced.

2.2. On the importance of the choice of the base flow

Implicitly, all the linear analyses over the years by Thorsness et al. (1978), Abrams & Hanratty (1985), Fourrière et al. (2010), Charru et al. (2013) and Claudin et al. (2017) were derived from the base flow solution of the inner region of the boundary layer configuration. Actually, with the use of Prandtl's mixing length model (A2), the linear analyses were made on a semi-infinite domain covering the viscous sublayer, the buffer region and the logarithmic region. The obtained perturbation is therefore included in this domain, without any interaction with the outer region as long as the upper boundary condition imposes a zero perturbation field. Additionally, the problem is then independent of the friction Reynolds number and only depends on the dimensionless wavenumber $\alpha ^+=\mathcal {R}^{-1}$. However, the reference experiments of Hanratty et al. were obtained in a rectangular channel of height $2\delta$ with $\delta \alpha ={\rm \pi}$. Therefore, the friction Reynolds number $\delta ^+$ may then influence the flow response to the wall deformation, and the validation of the results obtained from Hanratty's experiments in a channel may be questioned. To elucidate this issue, we consider a modified version of our code with a mixing length model adapted to channel flow configuration and using the Nikuradse formula

(2.3)\begin{equation} l = \delta \left(0.14-0.08\left(1-\frac{z}{\delta}\right)^2-0.06 \left(1-\frac{z}{\delta}\right)^4\right) \left(1-\exp\left(-\frac{\sqrt{\tau_{xz}}z}{\nu A}\right)\right) \end{equation}

For $\alpha \delta ={\rm \pi}$, corresponding to Hanratty's experiments, similar results (figure 2) are obtained with both versions of the code when $\beta$ is increased to $45$ in the channel configuration. Considering the existing dispersion for the experiments, both results are satisfactory. When $\alpha \delta$ is lowered or increased by a factor of $2$, the magnitude $\beta$ of the Hanratty correction $\mathcal {C}$ must be modified accordingly to recover the experimental data. There is a real influence of the friction Reynolds number on the results but it can be compensated by adjusting $\beta$. It is nevertheless important to note that both versions of the code with the respective mixing length models (A2) and (2.3) provide close results for $\mathcal {R}<500$ ($\alpha ^+>0.002$) for a common reference value $\beta =35$, whatever the values of $\alpha \delta$. Therefore, the dependence to the friction Reynolds number $\delta$ in the transitional regime is small and the linear responses obtained by considering the inner region of a boundary layer can be legitimately compared with measurements or computations obtained in channel flow configurations. The results presented below have all been produced by the code based on the inner boundary-layer region to be consistent with previous studies.

2.3. The role of the vorticity

Another remarkable aspect in the evolution of the wall shear stress phase is the influence of the vorticity. The penetration depth $\delta _i$ depends on the Reynolds number $\mathcal {R}$ and its definition (Charru & Hinch 2000) is given by the vorticity disturbance $\varpi =\hat {u}_{,\eta }-{\rm i}\hat {w}$ at the wall (see Appendix A). The penetration depth must not be seen as the distance to the wall where the perturbation is not zero but a measure of the distance over which the vorticity acts. Actually, the perturbation fields for the velocity and the pressure are not zero above $\delta _i$ but the vorticity is. Figure 3 depicts the normalised vorticity profiles for $\mathcal {R} \in [10,1000]$. Vorticity peaks, almost independent of $\mathcal {R}$, are clearly visible around $z^+=7$ before the profiles tend to zero. The disturbance field can be divided into a vortical region, near the wall, and a non-vortical region far from the wall. In the non-vortical region, the phases of the perturbations are nearly constant and without offsets from the corrugated wall. Below, the induced vorticity impacts on the profiles and phase shifts appear. The vortical region has a determining influence on the evolution of $\psi _{\tau }$.

Figure 3. Vorticity disturbance profiles. Dark blue to light blue lines indicate increased Reynolds number $\mathcal {R}={10,100,200,500,700,1000}$.

2.4. Turbulent closure for the linearised energy equation

To tackle dissolution or melting problems, Claudin et al. (2017) introduced an additional transport equation for a passive scalar in the linear analysis. The model was intended to be applicable to a wide range of applications using a Robin boundary condition at the wall. In the present context, in order to compare results of the linear analysis with numerical Navier–Stokes simulations, the considered passive scalar is the total enthalpy associated with the linearised energy equation (A10). Again, for the sake of comparison with numerical simulations, the boundary condition at the wall is a Dirichlet type condition where the enthalpy is imposed. For large values of wall heat flux, the dissipation can be neglected and the energy equation (A10) reduces to an advection–diffusion equation identical to the dissolution equation considered by Claudin et al. (2017). The model (A10) is representative of ablative materials for which, in the context of re-entry vehicles, the surface regression may be directly related to the energy equation or to an oxidiser concentration transport equation (White & Grabow 1973).

The main difference with Claudin et al. (2017) lies in the closure relation for the turbulent scalar flux, which here is the turbulent heat flux (A11). Claudin et al. (2017) considered that the mixing length for the turbulent scalar flux, denoted $l_{\theta }$, can be simply taken equal to $l$. For this study, a more general form (Cebeci & Smith 1974) for $l_{\theta }$ is retained by separating the damping functions for the velocity and the enthalpy

(2.4)\begin{equation} l_{\theta} = \kappa z \left(1-\exp\left(-\frac{z \sqrt{\tau_{xz}}}{\nu A}\right)\right)^{{1}/{2}} \left(1-\exp\left(-\frac{z \sqrt{\tau_{xz}}}{\nu A_{\theta}}\right)\right)^{{1}/{2}} .\end{equation}

The mixing length disturbance $\hat {l}_{\theta }$ is given by

(2.5)\begin{align} \hat{l}_{\theta} &=-\kappa\left(1-\exp\left(-\frac{\mathcal{R}\eta}{A^0}\right)\right)^{{1}/{2}} \left(1-\exp\left(-\frac{\mathcal{R}\eta}{A^0_{\theta}}\right)\right)^{{1}/{2}} \nonumber\\ &\quad \times \left[1 +\frac{1}{2}\frac{\displaystyle\exp\left(-\frac{\mathcal{R}\eta}{A^0}\right)}{\displaystyle 1-\exp\left(-\frac{\mathcal{R}\eta}{A^0}\right)} \left(\frac{\mathcal{R}\eta}{A^0}-\frac{\mathcal{R}\eta^2}{A^0} \left(\frac{\hat{\tau}_{xz}}{2}-\beta\mathcal{C}\right)\right) \vphantom{\frac{\displaystyle\exp\left(-\frac{\mathcal{R}\eta}{A^0_{\theta}}\right)}{1-\exp \displaystyle\left(-\frac{\mathcal{R}\eta}{A^0_{\theta}}\right)}}\right. \nonumber\\ &\quad \left. +\frac{1}{2}\frac{\displaystyle\exp\left(-\frac{\mathcal{R}\eta}{A^0_{\theta}}\right)}{1-\exp \displaystyle\left(-\frac{\mathcal{R}\eta}{A^0_{\theta}}\right)} \left(\frac{\mathcal{R}\eta}{A^0_{\theta}}-\frac{\mathcal{R}\eta^2}{A^0_{\theta}} \left(\frac{\hat{\tau}_{xz}}{2}-\beta_{\theta}\mathcal{C} - \epsilon_{\theta} \frac{\hat{\tau}_{xz}}{2}\right)\right)\right] . \end{align}

The introduction of a second damping function in (2.4) makes it possible to introduce an additional correction to $\hat {l}_{\theta }$ in (2.5). From Cebeci & Smith (1974), we have $A_{\theta }^0=30$; $A_{\theta }$ is made dependent on $\hat {\tau }_{xz}$ with a coefficient $\epsilon _{\theta }=({2}/{A_{\theta }^0})({\partial {A_{\theta }}}/{\partial {\hat {\tau }_{xz}}})$. The dependence of $A_{\theta }$ on $\mathcal {C}$ is taken to be identical to that of $A$ in (2.1) and in the following we take $\beta _{\theta }=\beta =35$. The results obtained with the model retained by Claudin et al. (2017) are recovered when ${A_{\theta }^0}=A^0=26$ and $\epsilon _{\theta }=0$.

3. Navier–Stokes computations

3.1. The RANS computations

To enlighten the impact of the turbulent closure on the forced response, several RANS computations were performed. The numerical procedure is based on the second-order compressible finite volume code named CEDRE (Aupoix et al. 2011; Scherrer et al. 2011), developed at ONERA and designed for unstructured grids. The computational domain is a two-dimensional periodic channel where $\alpha \delta ={\rm \pi}$. In order to respect Hanratty's experimental conditions, the sinusoidal profile was only applied on the bottom wall. Constant and homogeneous source terms were added to reproduce the mean pressure gradient and to balance the energy budget. A constant temperature was imposed as a boundary condition at the walls so that the induced fluxes compensate for the energy source term. The source terms were designed to respect as much as possible the incompressibility assumption. The density fluctuations were found to be three to four orders of magnitude below the velocity and pressure fluctuations. Eight configurations with various values of the kinematic viscosity $\nu$ were explored, corresponding to $\mathcal {R}\approx \{100,150,200,300,400,500,700,1000\}$ covering the transitional regime.

Two turbulence models were used to analyse the impact of the order of the Reynolds stress closure. On the one hand, computations with the $k\unicode{x2013}\omega$ model (Menter 1994), $k$ being the turbulent kinetic energy and $\omega$ the specific dissipation, were performed to characterise the influence of the Boussinesq hypothesis (A3) while, on the other hand, the elliptic blending Reynolds stress model (EBRSM) turbulence model (Manceau & Hanjalić 2002; Manceau 2015) was retained to obtain representative results of the second moment closure. The Boussinesq hypothesis is expected to have a significant impact on the streamwise momentum balance (A1) through the term $\tau _{xx}-\tau _{zz}$ in the transitional regime. With the Boussinesq hypothesis $\overline {u'^2}-\overline {w'^2}$ is made proportional to ${\partial {u}}/{\partial {x}}$, which is not true with second-order models. In particular, the exact production term for $\overline {u'^2}-\overline {w'^2}$ is $P_{xx}-P_{zz}=-4\overline {u'^2}({\partial {u}}/{\partial {x}})-2\overline {u'w'}({\partial {u}}/{\partial {z}}-{\partial {w}}/{\partial {x}})$ and suggests a dependence on the shear stress $\overline {u'w'}$ for the growth of $\overline {u'^2}-\overline {w'^2}$. At the first order with respect to the wall oscillation, the production term $P_{xx}-P_{zz}$ is not only ruled by the pressure induced velocity gradient ${\partial {u}}/{\partial {x}}$ but also by the shear stress $-\overline {u'w'}$. The objective of the computations is to highlight the effects of these differences on the evolution of $\psi _{\tau }$ with respect to $\mathcal {R}$.

The closure relations for the turbulent heat fluxes $\overline {u'_ih'}$ completely differ between the $k\unicode{x2013}\omega$ model and the EBRSM. The standard approach associated with eddy viscosity models such as the $k\unicode{x2013}\omega$ model is to make use of a simple gradient diffusion hypothesis (SGDH) with a turbulent thermal diffusivity including a constant turbulent Prandtl number $Pr_t$, in a similar manner to (A11) for the mixing length model of the linearised problem. In all the following $k\unicode{x2013}\omega$ computations, $Pr_t$ is set to $0.9$. In the context of second-order models, several approaches can be contemplated but the most commonly employed model relies on the generalised gradient diffusion hypothesis (GGDH) with the relation taken from Daly & Harlow (1970), $-\overline {u_ih'}=c_{\theta }\xi _t\overline {u_i'u_j'}({\partial {h}}/{\partial {x_j}})$. The turbulent time $\xi _t$ derives from the turbulent kinetic energy and its dissipation. The EBRSM was run with the classical value $c_{\theta }=0.22$, close to that recommended by Dehoux, Benhamadouche & Manceau (2017). The choice for the closure relation of $\overline {u'_ih'}$ has a considerable influence on the enthalpy perturbation field and the wall heat flux $\phi _w$. A close look to the expressions for the streamwise component $\overline {u'h'}$ for both models SGDH and GGDH reveals the influence of the shear stress $\overline {u'w'}$. The GGDH closure relation for a non-parallel bi-dimensional flow gives $\overline {u'h'}^{GGDH}=-c_{\theta }\xi _t\overline {u'^2}({\partial {h}}/{\partial {x}})-c_{\theta }\xi _t\overline {u'w'}({\partial {h}}/{\partial {y}})\approx \overline {u'h'}^{SGDH}-c_{\theta }\xi _t\overline {u'w'}({\partial {h}}/{\partial {y}})$ since $\xi _t \overline {u'^2} \propto \nu _t$. The shear stress is known to be affected by the wall deformation which means that, at the first order, the turbulent heat flux will thus behave differently between the SGDH and the GGDH models. In the transitional regime, the wall heat flux $\phi _w$ depends on the contribution of the turbulent heat flux in the energy budget and ultimately its phase $\psi _{\phi }$ with respect to the corrugated wall will be influenced by the choice of the closure relation.

3.2. The DNS computations for validation

The experimental data of Hanratty et al. do not allow a comprehensive examination of all of the aspects regarding the perturbations due to the wall waviness. There are no available data on heat transfer at the wall. For applications, the analysis of the energy budget is determining since the wall regression is most often driven by transfers at the wall that can be represented without any loss of generality by heat transfer, as recalled in § 2.4. To access such data, DNS were conducted with the spectral difference Navier–Stokes solver named JAGUAR (Chapelier, Lodato & Jameson 2016) and developed at ONERA and CERFACS. The code is designed to handle triangle (Veilleux et al. 2022a) or tetrahedral elements (Veilleux et al. 2022b) but all the presented computations were performed with a fourth-order discretisation scheme using hexahedral elements. Time integration is made with a low-dissipation low-dispersion sixth-order Runge–Kutta scheme. The computational domain is $[0,3\lambda ]\times [0,6\delta ]\times [\zeta _0\cos (\alpha x),\delta ]$ with $\alpha \delta ={{\rm \pi} }/{2}$. The streamwise extent of the domain is $12\delta \approx 4{\rm \pi} \delta$ that fits the usual requirements for periodic channel flow simulations. A constant source term is added on the momentum equation that sets the friction velocity $u_{\tau }$. The wall temperature is kept constant and the level of the mean heat flux on the wall is determined by the balance with the viscous and turbulent dissipation. As a consequence the wall heat flux is $\phi _w=\rho u_{\tau }^2 U_b$, with $U_b$ the bulk velocity, providing rather low values of $\phi _w^*=U_b^+$. Two mesh resolutions are used depending on the targeted Reynolds number. The numbers of solution points are $240\times 240 \times 160 \approx 9M$ and $320\times 320 \times 240 \approx 24M$. With a fourth-order discretisation, the mean $y^+$ values in the wall cells are found to stay between $0.25$ and $0.5$. The friction Reynolds numbers $\delta ^+$ range from $150$ to $500$ with $\mathcal {R}\in \{100,150,200,300\}$. Here again, the velocity field is nearly divergence free and the density fluctuations are several orders of magnitude lower than the velocity and pressure perturbations. The amplitude of the wall ripple is chosen to give $\zeta _0^+ \in [2.9,6.6]$, ensuring linear behaviours with $\alpha \zeta _0$ always less than $0.03$.

These DNS configurations cannot be directly considered as a complementary material to the experimental results of Hanratty et al. since $\alpha \delta$ is twice as small in the computations as in the experiments. However, it was shown in § 2.2 that, for $\mathcal {R}<500$, the phase shift $\psi _{\tau }$ is hardly affected by this change in the product $\alpha \delta$. This choice for $\alpha \delta$ is a compromise between representativeness and cost. The main purpose of these simulations is to serve as a reference for RANS computations and the linear analyses, especially concerning the heat transfer. For this reason, RANS computations were also performed with strictly similar conditions. All computations used air as the fluid with perfect gas assumptions and, given the temperature levels encountered, the specific heat capacity $C_p$ can be reasonably considered constant. The computed temperature fields are directly comparable to the enthalpy fields. We note by $\theta$ the temperature difference with the wall and $\theta ^+={\theta }/{\theta _{\tau }}$ the associated dimensionless variable, where $\theta _{\tau }={-\phi _w}/{\rho C_p u_{\tau }}$ is the friction temperature. Mean velocity profiles $\langle u^+ \rangle$ (figure 4) compare favourably between the different computations for all Reynolds numbers, even though the $k\unicode{x2013}\omega$ model underestimates the profiles in the buffer layer. The reference data of Hoyas & Jiménez (2008) obtained in non-deformed channels are also depicted to prove the validity of the DNS computations presented here. Second moments also agree between the two DNS datasets. The DNS mean temperature profiles are well reproduced by the EBRSM while the $k\unicode{x2013}\omega$ model tends to underpredict the profiles above the linear region.

Figure 4. Mean velocity (blue) and temperature (orange) profiles. Empty symbols ($\circ$,$\scriptstyle \square$) are DNS results while solid lines (EBRSM) and dashed lines ($k\unicode{x2013}\omega$) present RANS computations. The full black symbols ($\diamond$) are DNS results from Hoyas & Jiménez (2008) at $Re_{\tau }=180$ and $Re_{\tau }=550$: (a) ${\mathcal {R}}_{\tau }=150$; (b) ${\mathcal {R}}_{\tau }=300$.

4. Analysis and discussion

4.1. Influence of the turbulent closures on RANS computations

The narrow differences in the mean quantities visible in figure 4 actually hide more vast discrepancies in the perturbation fields, which increase with the Reynolds number $\mathcal {R}$. Profiles of the velocity and temperature perturbation fields were extracted at several streamwise location ${x}/{\lambda }$ and plotted in figures 5 and 6. The amplitudes of the perturbation are divided by a factor $2$ when $\mathcal {R}$ is doubled, in accordance with the linear expansion (A4) stating that any quantity $q$ is such that ${(q^+-\langle q^+ \rangle )}/{\zeta _0^+} \propto \alpha ^+=\mathcal {R}^{-1}$. It is immediately apparent that the EBRSM compares better with the DNS results than the $k\unicode{x2013}\omega$ model. The agreement is better for velocity perturbations than for temperature perturbations, where a noticeable difference exists below $z^+=20$. Despite a good overall trend, the perturbation profiles presented by the $k\unicode{x2013}\omega$ model are lagged behind those of DNSs with smaller amplitudes. The higher the Reynolds number, the larger the lag. Another notable point that emerges from these figures is that the ordering between the profiles is modified from the centre of the channel to the wall. Figures 5 and 6 again illustrate the division between vortical and non-vortical regions. Around the centre of the channel, the phase of the perturbed field is not altered with respect to the wall and the ordering between profiles is aligned with the wall locations, i.e. in-phase or anti-phase, depending on the sign of the perturbation. Conversely, near the wall, the ordering is modified by the phase of the perturbed field. Moreover, DNS and RANS calculations have also revealed a perturbation peak in the velocity profiles around $z^+=10$, consistent with the vorticity peak revealed by the linear analysis (figure 3). A similar peak is also visible in the temperature profiles, but is less pronounced due to the high levels of perturbations observed in the non-vortical region. The wall shear stress disturbances ${(\tau _w^+-\langle \tau _w^+ \rangle )}/{\zeta _0^+}$ of figure 7 corroborate the previous observations with $k\unicode{x2013}\omega$ predictions delayed compared with those of DNS while the EBRSM provides better agreement. For the wall heat flux disturbances ${(\phi _w^+-\langle \phi _w^+ \rangle )}/{\zeta _0^+}$ presented in figure 7, the $k\unicode{x2013}\omega$ model underestimates the amplitudes and is not able to recover the phase shift. The EBRSM greatly improves the results but the phase shift on $\phi _w$ is a bit overpredicted. The RANS results for the wall shear stress phase $\psi _{\tau }$ are also reported in figure 2. The closure relations of the RANS computations are manifestly responsible for the prediction accuracy and the results evidence the failure of the Boussinesq hypothesis, as expected. Even though the wall deformation is very small, ensuring a linear behaviour of the perturbation, the flow field is heavily affected by the turbulent modelling. The error is even more pronounced on the perturbed temperature field and the wall heat flux. As explained above, the good behaviour of the EBRSM compared with the $k\unicode{x2013}\omega$ model is essentially due to the representation between the Reynolds stress difference $\overline {u'^2}-\overline {w'^2}$. Figure 8 shows the mean and disturbed profiles of $\overline {u'^2}^{~+}-\overline {w'^2}^{~+}$ and $-\overline {u'w'}^{~+}$ obtained with the EBRSM calculations compared with those from the DNS for $\mathcal {R}=300$. Although the forced response does not match that yielded by DNS, the profile of $\overline {u'^2}-\overline {w'^2}$ at leading order is in good agreement with DNS results, while for $k\unicode{x2013}\omega$ calculations, (not shown here) the normalised stress difference at the leading order is $\langle \overline {u'^2}^{~+}-\overline {w'^2}^{~+}\rangle =4\langle \nu _t ({\partial {u}}/{\partial {x}})\rangle =0$. Figure 8 also indicates that the perturbations due to the wall in the diagonal stress difference $\overline {u'^2}-\overline {w'^2}$ are four to five times larger than those induced in the shear stress $-\overline {u'w'}^{~+}$. The results show that the term ${\partial {(\tau _{xx}-\tau _{zz})}}/{\partial {x}}$ has a magnitude five times smaller than that of the term ${\partial {\tau _{xz}}}/{\partial {z}}$ in the streamwise momentum equation (A1). In the end, the ${\partial {(\tau _{xx}-\tau _{zz})}}/{\partial {x}}$ term contributes approximately 20 % in the budget of the momentum equation at the first order, showing its critical importance.

Figure 5. Profiles of velocity perturbations at stations $x/\lambda =0.0$ (blue), $x/\lambda =0.2$ (purple), $x/\lambda =0.4$ (green), $x/\lambda =0.6$ (orange) and $x/\lambda =0.8$ (red). Symbols are DNS results, solid lines present the RANS computations with the EBRSM while the dashed lines stand for the $k\unicode{x2013}\omega$ results: (a) ${\mathcal {R}}_{\tau }=150$; (b) ${\mathcal {R}}_{\tau }=300$.

Figure 6. Profiles of temperature perturbations. The legend is identical to that of figure 5: (a) ${\mathcal {R}}_{\tau }=150$; (b) ${\mathcal {R}}_{\tau }=300$.

Figure 7. Wall shear stress (left/blue colour) and wall heat flux disturbances (right/orange colour) at $\mathcal {R}=150$ (a) and $\mathcal {R}=300$ (b). Symbols are DNS results. Solid lines are the RANS computations with the EBRSM and the dashed lines represent the computations with the $k\unicode{x2013}\omega$ model.

Figure 8. (a) Mean profiles of the Reynolds stress difference $\overline {u'^2}^{~+}-\overline {w'^2}^{~+}$ (blue) and the shear stress $-\overline {u'w'}^{~+}$ (orange) for $\mathcal {R}=300$. Symbols are the DNS results and lines stand for the EBRSM computations. Corresponding forced responses profiles ${(\overline {u'^2}^{~+}-\overline {w'^2}^{~+}-\langle \overline {u'^2}^{~+}-\overline {w'^2}^{~+}\rangle )}/{\zeta _0^+}$ (b) and ${(-\overline {u'w'}^{~+}+\langle \overline {u'w'}^{~+}\rangle )}/{\zeta _0^+}$ (c) at several stations $x/\lambda$. Lines and symbols are those used in figure 5.

The RANS results are now used to further examine the results of the linear analysis and assess the effect of the Hanratty correction $\mathcal {C}$ on the prediction of the phase shift of the wall shear stress and the wall heat flux.

4.2. Achievements and limitations of the Hanratty correction

Previous work by Abrams & Hanratty (1985), Charru et al. (2013) and Claudin et al. (2017) proved the effectiveness of the correction $\mathcal {C}$ in recovering the wall shear stress phase evolution with respect to the wavenumber (solid blue line in figure 2). Although very efficient, this correction suffers from two main limitations. The first one is related to the application of the correction in the mixing length model. The RANS computations highlighted the failure of the Boussinesq hypothesis to predict the stress difference $\tau _{xx}-\tau _{zz}$, which is then of the order of the perturbation $O(\alpha \zeta _0)$, in the streamwise momentum equation (A1). However, the Hanratty correction acts on the shear stress $\tau _{xz}$ through the modification of the mixing length. In other words, the Hanratty correction does not correct the problematic term but balances the streamwise momentum equation, and, in that sense, it can be viewed as an ad hoc palliative to the failure of the Boussinesq hypothesis. The second limitation comes from the use of a relaxed pressure gradient to drive the correction $\mathcal {C}$. The RANS and DNS calculations have evidenced the role of the mean vorticity of the flow in creating the turbulent stresses that ultimately lead to the observed phase shift in the wall shear stress. However, the pressure gradient does not enter the vorticity equation and is not a relevant variable to control turbulence. Furthermore, the pressure gradient is not involved in the Reynolds stress transport equations, which does not prevent the EBRSM computations from correctly reproducing the phase shift of the wall shear stress. Despite these limitations, the Hanratty correction is very useful and effective for linear analyses.

A further demonstration of the positive impact of the correction $\mathcal {C}$ is shown in figure 9, where the amplitudes of the wall shear stress disturbance are presented. In the linear analysis the wall shear stress fluctuation ${(\tau _w^+-\langle \tau _w^+ \rangle )}/{\zeta _0^+}$ is given by $\alpha \hat {\tau }_{xz}(0)$ according to (A4). Calculations of the linear response with $\mathcal {C}$ and, to a lesser extent, the EBRSM results, follow the measurements remarkably well, while the $k\unicode{x2013}\omega$ model and results of the linear analysis without the Hanratty correction move further apart as $\alpha ^+$ is decreased. We now focus on the use of the Hanratty correction in the closure relation for turbulent heat flux of the linearised energy equation detailed in § (2.4). In the mixing length disturbance $\hat {l}_{\theta }$ (2.5), $\mathcal {C}$ is considered twice with respect to the two van Driest numbers $A$ and $A_{\theta }$. An additional dependence on $\hat {\tau }_{xz}$ was introduced for the van Driest number $A_{\theta }$. Best agreements were obtained with $\epsilon _{\theta }=4$. The results of the linear analysis for the evolution of the phase of the wall heat flux $\psi _{\phi }$ with respect to $\alpha ^+$ are shown in figure 10 and compared with RANS computations. Results corresponding to the original model proposed by Claudin et al. (2017) ($A_{\theta }=26$ and $\epsilon _{\theta }=0$) are also reported in figure 10. Values of $\psi _{\phi }$ are shifted from $180^{\circ }$ when the sign of $\phi _w^*$ is changed. When $\vert \phi _w^* \vert$ is large enough, practically when $\vert \phi _w^* \vert >100$, the dissipation term $\hat {u}$ of the equation for the mean enthalpy (A14) is almost negligible and the equation is symmetrical with respect to $\phi _w^*$. The Navier–Stokes computations with the $k\unicode{x2013}\omega$ model provide values of $\psi _{\phi }$ in good agreement with the linear analysis obtained with the frozen turbulence assumption, consistent with the observation made on $\psi _{\tau }$ in figure 2. Results equivalent to those of Claudin et al. (2017) provide overestimated phase values of approximately $40^{\circ }$ whereas, with $A_{\theta }^0=30$ and especially $\epsilon _{\theta }=4$, the linear forced responses match those of the EBRSM computations. This means that the Hanratty correction has a beneficial impact on $\hat {l}_{\theta }$ but it is not sufficient. An additional correction on $A_{\theta }$, with $\epsilon _{\theta }=4$, is required to recover the results obtained with the EBRSM computations.

Figure 9. Amplitude of the wall shear stress perturbation. The filled black circles are measurements of Hanratty et al.. Square symbols are the RANS computations with the EBRSM (blue) and the $k\unicode{x2013}\omega$ (orange) model. The solid lines are the results of the linear analysis with the Hanratty correction (blue) and when the frozen turbulence assumption is used (orange).

Figure 10. Wall heat flux disturbance phase $\psi _{\phi }$ as a function of the wavenumber $\alpha ^+$. Blue symbols are the phase computed with the EBRSM for $\phi _w^*=400$ ($\Diamond$) and $\phi _w^*=-400$ ($\square$). The orange square symbols are the results obtained with the $k\unicode{x2013}\omega$ model. The solid lines are the corresponding results of the linear analysis with (blue line) all corrections activated ($A=26$, $A_{\theta }=30$, $\beta =35$, $\epsilon _{\theta }=4$) and with (orange line) the frozen turbulence hypothesis ($A=26$, $A_{\theta }=30$, $\beta =0$, $\epsilon _{\theta }=0$). The black dashed line presents the results corresponding to the approach followed by Claudin et al. (2017) for $l_{\theta }$ ($A=26$, $A_{\theta }=26$, $\beta =35$, $\epsilon _{\theta }=0$). The thin horizontal dashed line corresponds to $\psi _{\phi }=-90^{\circ }$.

In figure 11 the comparison of the amplitude of the wall heat flux disturbance points out several divergences. The disturbances of the wall heat flux ${(\phi _w^+-\langle \phi _w^+ \rangle )}/{\zeta _0^+}$ obtained in the linear analysis, i.e. $-\alpha \hat {f}(0)$, are smaller than those of the RANS computations. The results produced by the $k\unicode{x2013}\omega$ model and the results of the linear analysis with the frozen turbulence assumption exhibit almost the same trends whereas the EBRSM and the linear analysis results diverge as $\alpha ^+$ decreases. This may be due to the closure relation used for the turbulent heat fluxes $-\overline {u_i'h'}$. The EBRSM uses the GGDH assumption while the linear analysis makes use of a SGDH hypothesis and is impacted by the Hanratty correction $\mathcal {C}$ in $\hat {l}_{\theta }$ (2.5). The results obtained with $A_{\theta }^0=26$ and $\epsilon _{\theta }=0$ are not better. For the energy equation, the Hanratty correction is not able to compensate the approximation made in the modelling of the turbulent heat fluxes. This is not surprising since $\mathcal {C}$ was implicitly designed to correct only for the misrepresentation of the Reynolds stresses. Although imperfect, the linear analysis using the model described in § 2.4 for the energy equation allows a good prediction of $\psi _{\phi }$. However, it was not possible with this type of closure (2.5) for $\hat {l}_{\theta }$ to also obtain a satisfying prediction of the wall heat flux amplitude.

Figure 11. Amplitude of the wall heat flux perturbation. Square symbols are the RANS computations with the EBRSM (blue) and the $k\unicode{x2013}\omega$ (orange) model. The solid and dashed lines are the results of the linear analysis. The blue lines correspond to results of the linear approach with $\phi _w^*=400$ (dashed) and $\phi _w^*=-400$ (solid). The orange line presents the analysis performed with the frozen turbulence assumption. The dashed black line shows the results obtained with $A_{\theta }=A_{\theta }^0=26$ and $\epsilon _{\theta }=0$.

4.3. Linear stability of an ablative surface

The surface elevation is now a function of time $\zeta (x,t)=\zeta _0 \exp ({ (\sigma _w t+ {\rm i}\omega _w t + {\rm i}\alpha x)})$ and is assumed to be ruled by the ablation process and controlled by the wall heat flux. For moving surfaces, the critical layer, below which the flow propagates more slowly than the surface, has a crucial importance on the flow dynamics (Belcher & Hunt 1998). For our re-entry applications (see Appendix B), the surface speed ${\omega _w}/{\alpha }$ is low compared with the friction velocity. In this slow waves regime (${\omega _w}/{\alpha u_{\tau }}\lesssim 15$) the critical layer is thin and plays no significant dynamical role. In other words, only the temporal growth rate $\sigma _w$ matters and controls the surface regression in direction $z$.

The model detailed in the Appendix A.2 can be applied to dissolution or melting problems since the energy equation produces similar results to the advection–diffusion equation used by Claudin et al. (2017) when $\vert \phi _w^* \vert$ is large. Any solid surface can be decomposed into a series of sinusoidal profiles and the linear response of the flow will be the combination of the responses for each wavenumber. The surface regression is assumed to be proportional to the wall flux (Claudin et al. 2017). Dropping the homogeneous part of the flux, the evolution of the elevation at the first order is ruled by ${\partial {\zeta }}/{\partial {t}}=-r u_{\tau }^3 \zeta _0 \alpha \vert \hat {f}(0) \vert \exp ({ ({\rm i}\alpha x+{\rm i}\psi _{\phi })})$, with $r$ a constant proportionality factor (s$^{2}$/m$^{2}$), controlling the regression rate. The temporal growth rate of the surface elevation is then governed by the real part of the dispersion relation, i.e. $\sigma _w=-r u_{\tau }^3 \alpha \vert \hat {f}(0) \vert \cos (\psi _{\phi })$. Function $\sigma _w(\alpha )$ changes sign when $\vert \psi _{\phi } \vert$ crosses the horizontal line $\psi _{\phi }=90^{\circ }$. In the case of a negative wall heat flux, the horizontal line $\psi _{\phi }=-90^{\circ }$ is plotted in figure 10. When the Boussinesq hypothesis is used without the Hanratty correction in the linear analysis and for the computations with the $k\unicode{x2013}\omega$ model, $\psi _{\phi }$ is always less than $-90^{\circ }$ and $\sigma _w$ remains negative for all wavenumbers $\alpha ^+$. For the EBRSM results or for the linear responses involving the correction $\mathcal {C}$, $\sigma _w$ becomes positive for $\alpha ^+\approx 0.006$. All wavenumbers below $\alpha ^+ \approx 0.006$ are unstable, in the range of wavenumbers covering the transitional regime. However, the growth rate $\sigma _w$ quickly decreases as $\alpha ^+$ decreases, mainly due to its proportionality with $\alpha$. In figure 12, growth rates $\sigma _w$ (normalised) obtained in the RANS computations and in the linear analysis are depicted with respect to $\alpha ^+$. In the $k\unicode{x2013}\omega$ computations and in the linear analysis with the frozen turbulence assumption, the growth rates are always negative. Both models predict stable modes regardless of $\alpha ^+$. But the EBRSM and the linear approach show an unstable region where $\sigma _w>0$ and the presence of a peak. The wavenumber associated with this peak indicates the most unstable mode for which the surface time growth rate is the highest. The error in the prediction of the amplitude of the wall heat flux with the linearised model using (2.5) for the closure of the turbulent heat flux leads to a shift in the position of the peak. The linear analysis indicates a peak at $\alpha ^+=2.4 \times 10^{-3}$ ($\mathcal {R}= 417$) whereas it is found at $\alpha ^+\approx 4\times 10^{-3}$ ($\mathcal {R}=250$) by the EBRSM. The location of the peak is almost independent of $\phi _w^*$ and is not modified by the sign of $\phi _w^*$ as long as $\vert \phi _w^* \vert > 100$. The Prandtl number $Pr$ has a limited influence on the peak position in the linear approach. The same tendency is expected in Navier–Stokes computations. The peak is moved to higher values of $\alpha ^+$ as $Pr$ is increased and, for example, when $Pr=100$, the peak is located at $\alpha ^+=4.2\times 10^{-3}$.

Figure 12. Normalised growth rate ${\sigma _w \phi _w}/{r\rho }$ with respect to $\alpha ^+$ in logarithmic and linear scales. Square symbols are the RANS computations with the EBRSM (blue) and $k\unicode{x2013}\omega$ (orange) model. The corresponding dashed lines are splines computed from the data. The solid lines are the results of the linear approach with (orange) and without the frozen turbulence assumption.

Thomas (1979) presented evidence in support of the hypothesis that the scalloping of soluble surfaces may be attributed to wall turbulence. By analysing bed morphologies where scallops occur, he showed that the longitudinal wavelength of the bedform is a multiple of the viscous length $\delta _{\nu }$ providing $\alpha ^+ \approx 6\times 10^{-3}$. The proportionality between these quantities was demonstrated over a range exceeding four decades of length and covering a wide variety of situations from the corrosive dissolution of steel (Schoch 1968; Schoch, Richter & Köhle 1969; Schoch, Richter & Effertz 1970a; Schoch et al. 1970b; Schuster 1971; Heimsch et al. 1978), brass (Sick 1972) and copper (Knutsson, Mattsson & Ramberg 1972), the plastic shear of bitumen (Brauer 1963) and aluminium (Brunton 1966) and the rippling of colloidal-particle deposits in a water main (Wiederhold 1949; Seiferth & Krüger 1950). In the context of atmospheric re-entry vehicles, the wavelength found in the TATER experiment (Hochrein & Wright 1976) aligns with the Thomas correlation. The orders of magnitude provided in Appendix B justify the use of the linear approach (A.1 and A.2) to study this type of flow, particularly with respect to compressibility effects. The location of the most unstable mode with the EBRSM computations or with the linear approach are closed to the value found in the Thomas correlation, confirming the role of turbulence in the occurrence of scallops. It is nevertheless premature to draw general conclusions from these results. Only the linear response was examined, with a high degree of hypothesis on the flow that restricts the scope of the approach. Further verification is needed to extend the approach to different types of erodible surfaces where scallops are observed. Nonlinear effects, notably related to flow separations, may also interfere in the scallop formation (Charru et al. 2013). This will certainly require further experimental or numerical data for validation. The results presented are a first step towards explaining the value of the slope of the Thomas correlation.

5. Conclusion and perspectives

The scallops observed on re-entry blunt bodies are similar to those encountered in many applications, the characteristic scale of which is given by the Thomas correlation of the viscous boundary-layer length. The study of these scallops was historically based on a linear analysis of the disturbances generated by a fixed wall corrugation on the inner region of a turbulent boundary. The success of this approach relies in particular on the use of the Hanratty correction, without understanding the underlying mechanisms requiring the intervention of this correction. Using RANS and DNS numerical simulations, an in-depth analysis of the perturbations generated by the corrugated wall has allowed us to clarify the implications of the different terms of the Navier–Stokes equations and to better understand the role of the Hanratty correction.

It is found that the disturbance profiles can be separated into two distinct regions. Away from the wall, the vorticity perturbation is zero and the velocity and temperature profiles are in phase with the wall undulation. In the vicinity of the wall, the vorticity disturbance is significant and a phase shift with respect to the wall is observed in the various perturbed quantities. The vorticity creation is directly related to the contribution $\overline {u'^2}-\overline {w'^2}$ in the streamwise momentum equation. The RANS computations using the $k\unicode{x2013}\omega$ model and the EBRSM, confronted with reference results from DNS, highlight the failure of the Boussinesq hypothesis in this context. The results for the velocity disturbances show that the $k\unicode{x2013}\omega$ calculations, which are based on the Boussinesq hypothesis, are not able to reproduce the DNS data correctly, unlike the EBRSM calculations, which are fairly accurate. The differences between the DNS results and the $k\unicode{x2013}\omega$ computations are even greater for the temperature profiles. The use of a SGDH closure for turbulent heat fluxes further increases the errors. In contrast, the EBRSM calculations, which use a GGDH closure, show very good agreement with the DNS calculations, notably for the parietal heat flux.

A comparative study of results from the linear analysis and RANS results highlights the role of the Hanratty correction. The latter serves in fact to compensate for the poor representation of the Reynolds stresses in the equations coming from the use of the Boussinesq hypothesis. The Hanratty correction was designed to act effectively on the momentum equation. Its indirect use in the energy equation does not make it possible to obtain the expected results for wall heat transfer. In particular, the phase shift and the amplitude of the wall heat flux fluctuation are poorly predicted by the linear approach, even with the Hanratty correction, unless a supplementary correction is also added to the mixing length governing the turbulent heat flux closure. Finally, the study of wall regression under the effect of an ablative flux is carried out. The surface elevation is supposed to be ruled by the wall heat flux and its growth rate, apart from the homogeneous contribution of the leading order, is governed by the phase shift and amplitude of the wall heat flux disturbance. When the Boussinesq hypothesis is used without compensation, the linearised problem is unconditionally stable. However, in the linear approach using the Hanratty correction and in the RANS EBRSM computations, the growth rate of the surface elevation is found to be positive for $\alpha ^+>0.006$ in the transitional regime. The most unstable mode is found for $\alpha ^+=2.4 \times 10^{-3}$ in the linear analysis and around $\alpha ^+=4 \times 10^{-3}$ in the EBRSM computations. The difference in location results from the errors made on the phase and amplitude in the linear analysis occurs because of the turbulent closure relations used. These values of the dimensionless wavenumber are close to that given by the Thomas correlation, providing a first indication on the mechanisms involved in the occurrence of the scallops in the linear phase.

Many questions are still open and studies are needed to evaluate the influence of compressibility, regression models including possible chemical reactions, real gas effects, roughness effects and finally nonlinear interactions. In parallel, as suggested in figure 1, a three-dimensional linear analysis taking into account surface curvature effects could provide additional information on the three-dimensional nature of the scallops.