Numerical tripping of high-speed turbulent boundary layers

The influence of turbulence inflow generation on direct numerical simulations (DNS) of high-speed turbulent boundary layers at Mach numbers of 2 and 5.84 is investigated. Two main classes of inflow conditions are considered, based on the recycling/rescaling (RR) and the digital filtering (DF) approach, along with suitably modified versions. A series of DNS using very long streamwise domains is first carried out to provide reliable data for the subsequent investigation. A set of diagnostic parameters is then selected to verify achievement of an equilibrium state, and correlation laws for those quantities are obtained based on benchmark cases. Simulations using shorter domains, with extent comparable with that used in the current literature, are then carried out and compared with the benchmark data. Significant deviations from equilibrium conditions are found, to a different extent for the various flow properties, and depending on the inflow turbulence seeding. We find that the RR method yields superior performance in the evaluation of the inner-scaled wall pressure fluctuations and the turbulent shear stress. DF methods instead yield quicker adjustment and better accuracy in the prediction of wall friction and of the streamwise Reynolds stress in supersonic cases. Unrealistically high values of the wall pressure variance are obtained by the baseline DF method, while the proposed DF alternatives recover a closer agreement with respect to the benchmark. The hypersonic test case highlights that similar distribution of wall friction and heat transfer are obtained by both RR and DF baseline methods.

on direct numerical simulation (DNS) and large-eddy simulation (LES) of the Navier-Stokes equations, are a suitable alternative to experiments for the study of such canonical flows.
Several numerical studies dealing with high-speed boundary layers have emerged in the last decades, probably starting with the work of Rai et al. [26] who first addressed the problem in a fully spatial setting. Pirozzoli and Bernardini [22,23] developed a DNS database for supersonic boundary layers under adiabatic conditions, which has been extensively used as a benchmark in later studies [43], which have also widened the Mach number range under scrutiny. Regarding hypersonic flow, pioneering DNS were carried out by Martín [18], with follow-up studies by Duan et al. [6,7]. DNS up to free-stream Mach number M 0 = 20 were carried out by Lagha et al. [15]. Zhang et al. [50,51] reported an extensive compilation of velocity, pressure and temperature statistics from M 0 = 2.5 to 14, for several wall cooling ratios, which is currently regarded as a reference in hypersonic flow studies.
A general issue encountered in all above-mentioned studies is the reliability and repeatability of the results. For instance, Poggie et al. [25] investigated effects related to grid resolution and domain width. A major source of discrepancy among different studies is associated with the imposition of the inflow conditions in spatially evolving simulations (e.g. Wu [47]). Sensitivity of boundary layer statistics to the upstream flow conditions is well known from studies dealing with incompressible flow [29,34], which also reflects experimental difficulties in achieving a fully developed state past tripping devices [8]. From a computational perspective, the most obvious way to overcome the problem has been using (more or less) long domains to allow attainment of an equilibrium state [33,34], and perhaps using several stages, whereby the solution at one stage is interpolated and fed as inflow condition to the next stage [51]. However, the use of very long domains may lead to unacceptable computational cost, which makes it certainly unsuitable for industrial problems, in which the development length shall be minimized. Several previous studies have therefore attempted to define an inflow length, namely the minimal distance from the inflow at which the flow properties become independent of the inflow forcing. This subject was dealt with regarding experiments [3], and regarding DNS data [28]. In the latter study, it was concluded that inflow effects tend to be concentrated in the wake region, but some effect is also observed in the near-wall region, with the conclusion that DNS data should be carefully scrutinized just like experiments. Schlatter and Örlü [29] investigated the effect of the influence of tripping on the development of a boundary layer with laminar inflow, concluding that tripping effects are mainly responsible for large spread of the numerical results. Wenzel et al. [43] attempted to extend the previous studies to compressible flows. By defining the inflow length based on fulfilment of mean momentum balance, those authors found that it increases monotonically with the free-stream Mach number, as also confirmed by Huang et al. [11]. The latter authors also noticed significant differences of numerical results obtained with the digital filtering (DF) and with the recycling/rescaling (RR) procedure, especially regarding the wake region. Adler et al. [1] and Dhamankar et al. [5] reported significant scatter across simulations with similar set-up and different inflow, concluding that the prediction of the wall properties is quite unreliable. Regardless of the turbulence seeding method, the inflow development length can be somewhat reduced by using sponge layers [16].
Given this background, it is clear that despite general acknowledgement, the problem of inflow effects on numerical simulations of high-speed wall-bounded flows is far from being quantitatively assessed. The goal of the present study is thus threefold. First, we aim at setting up a standardized procedure to evaluate achievement of equilibrium conditions in the numerical simulation of high-speed boundary layers. Second, and related to the previous item, we aim at developing a benchmark database to be used for the evaluation of the performance of inflow feeding techniques, which is devoid from inflow history effects, Third, we aim at developing improved versions of the most widely used inflow feeding techniques, namely the recyclingrescaling procedure [17,41,49] and the digital filtering procedure [14,39].
The paper is organized as follows: in Sect. 2, we provide generalities about turbulence seeding methods and suggest alternatives; in Sect. 3, we describe the numerical database, and in Sect. 4 we analyse the benchmark data set; the results of comparative tests with different inflow turbulence seeding are presented in Sect. 5, which serves as a basis for the estimation of the inflow length, as reported in Sect. 6.

Generating inflow turbulence
This study is focussed on the two most commonly used techniques for feeding inflow turbulence, namely the recycling/rescaling and the digital filtering procedure, and we discuss possible modifications. For later reference, all inflow properties (namely density ρ, velocity u i and pressure p) are split into mean and fluctuating parts, namely where 0 subscript is used to denote inflow quantities, overline and prime symbols denote ensemble-averaged quantities and fluctuations thereof, whereas tilde and double-prime symbols denote Favre-averaged properties and fluctuations thereof, e.g.φ = ρϕ/ρ, ϕ = ϕ −φ. Wall units, based on the friction velocity (u τ = √ τ w /ρ w , where τ w is the wall shear stress), and the viscous length scale (δ v = μ w /u τ /ρ w , with ρ w and μ w being the wall density and dynamic viscosity), will be denoted with the + superscript. The forthcoming paragraphs deal with specification of inflow fluctuations aiming to mimic realistic turbulence.

2.1
The recycling/rescaling procedure The procedure followed recycling/rescaling is sketched in Fig. 1. Following the original formulation of Lund et al. [17] and its compressible flow extension [49], the fluctuations of each flow variable are extracted at a suitable station (say x r ) and recycled to the inflow plane, after suitable rescaling.
The rescaling procedure is applied by dividing the boundary layer into two sub-layers: (i) the inner layer (superscript 'inn') where velocity is assumed to scale in wall coordinates (y + = y/δ v ) ; and (ii) the outer layer (superscript 'out') where flow properties scale in outer units (Y = y/δ, δ being the local boundary layer thickness). The fluctuation of a generic quantity (ϕ) is assumed to be a weighted combination of the inner-and outer-layer fluctuations where the weight function W (Y ) is defined as [17] with α = 4, b = 0.2. The inflow density and velocity fluctuations in each layer are rescaled from the recycling station according to  where γ is the rescaling parameter, γ = u τ 0 /u τ r , which we estimate as γ = (δ r /δ 0 ) κ , with κ = 0.13. The inflow conditions are softly enforced by carrying out characteristic decomposition at the computational boundary in which the incoming wave amplitudes are estimated from a target flow state [24], as defined in Eq. (4), with the additional assumption of uniform pressure.
To promote streamwise decorrelation, the inflow fluctuations at the spanwise location z are recycled from a staggered location (z + L z /2) at the recycling station [37], and random divergence-free disturbances with maximum amplitude 4% u 0 are added at the inflow to break any remaining symmetry. Following previous studies [22], the recycling station is placed at x r = 53δ 0 . Preliminary simulations carried out by moving the recycling station to x r = 30δ 0 have shown negligible differences in the computed flow statistics. This baseline form of the recycling/rescaling procedure is hereafter referred to as RR1.
Possible simplifications of the RR procedure are also considered here. Instead of taking a linear combination of inner-and outer-rescaled fluctuations, one can more simply consider a remapping of the wall-normal coordinate with a stretching factor defined to be δ ν,r /δ ν,0 close to the wall (to achieve correct rescaling in the inner layer) and δ r /δ 0 away from it (to achieve correct rescaling in the outer layer). Here, we consider a simple tanh blending of the typeŷ where Y c = 0.08 and = 1.1, the former defining the crossover between the inner and outer transformations, while the latter controlling the smoothness of the blending. The fluctuations at the inflow station are then obtained as follows: This variant is hereafter referred to as RR2. Finally, we also consider a simple variant which includes remapping of the wall-normal coordinate based only on the ratio of the boundary layer thicknesses δ r /δ 0 (hence,ŷ/y = δ 0 /δ r ) and which is hereafter referred to as RR3. It is noteworthy that recycling/rescaling procedures require specification of the mean flow distributions in (1). Different from the original formulation of Lund et al. [17], in which both mean flow and fluctuations are recycled, in the present RR implementation the inflow mean distributions are specified by the user and left unchanged during the simulation. In our experience, this approach prevents long-time numerical drift, thus accelerating statistical convergence. As a consequence, the inflow boundary layer thickness does not change in time, with the additional advantage that no continuous control of δ r /δ 0 is needed.

The digital filtering procedure
The digital filtering implementation relies on an extension of the technique originally introduced by Klein et al. [14], which makes use of the strong Reynolds analogy (SRA) [39], as introduced by Kempf et al. [13]. Specifically, synthetic velocity fluctuations are first generated from a white noise sample, which is then filtered based on a sequence of one-dimensional convolutions, thus obtaining a correlated signal in space and time, with arbitrarily prescribable integral length and time scales. Different wall-normal distributions of the spanwise integral length scale ( z ) are assigned for each velocity component. Following Xie and Castro [48], two length scales are used for the inner ( z,inn ) and for the outer ( z,out ) wall layer, which are suitably blended to give The use of two distinct length scales can potentially cause issues within the blending region, where eddies generated from digital filtering can be truncated. This potential problem is, however, not observed here as the Reynolds numbers are low enough that the outer length scale is selected throughout. Values of the integral length scales are reported for each velocity component in Table 1. As reported by Xie and Castro [48], the digital filtering technique has rather small sensitivity to variations in the integral length scales. The authors found that, for an incompressible channel flow, 50% change of y and z yields less than 10% change in the normal stresses, and less than 13% in the shear stress at a distance of 10δ 0 from the inflow. The wall-normal integral length scales are then selected as y (y) = 0.67 z (y). Finally, constant values of the longitudinal integral length scales x are assigned for each velocity component as in Table 1. The streamwise length scales are then converted to time scales using Taylor hypothesis, assuming that eddies are convected at the free-stream velocity. The resulting velocity fluctuations are then rescaled to match a desired wall-normal distribution of Reynolds stresses. Temperature fluctuations T are obtained by applying SRA (namely, T = −ũu /c p , where the tilde denotes Favre averages, and c p is the specific heat at constant pressure), and converted to density fluctuations assuming zero pressure fluctuations. This baseline implementation of the DF algorithm is hereafter referred to as DF1.
Modifications of the baseline DF algorithm are also considered in the present work. In preliminary attempts of improving the performance of DF, we found that removing the streamwise velocity fluctuations at the inflow yields a reduction in spurious pressure disturbances generated by the DF1 implementation. This variant is referred to as DF2 in the following. A further improvement includes introducing a suitable stream function for the cross-stream velocity fluctuations to make the inflow turbulence solenoidal. The stream function is defined such as and we assume (y, z, t) = C(y) r (y, z, t), where r (y, z, t) is obtained from application of a low-pass filter to a white noise generator, and C(y) is a suitable scaling function. The cross-stream velocity fluctuations are then obtained from Eq. (8). The shape of the scaling function is chosen so as to achieve a prescribed distribution of the wall-normal Reynolds stress, It is noteworthy that the variance term at the right-hand side of (10) can be computed beforehand, depending only on the spanwise length scale z in this formulation. As the present method hard enforces the distribution of the wall-normal velocity variance, the spanwise variance cannot be controlled independently to fit the target distribution. This variant is referred to as DF3 in the following. It should be pointed out that various choices can be made to introduce the stream function . Whereas defining the stream function based on ρu i0 would allow us to find a flow field that has ∂ρ/∂t = 0, defining the stream function based on u i0 helps us find a flow field that has zero dilatation. As the goal of the DF3 method is to minimize spurious acoustics, the second option is the most suitable one owing to the strong connection between dilatation and acoustics (while acoustic waves have nonzero ∂ρ/∂t). In fact, regardless of how we define the stream function , a fully developed field has significant ∂u/∂ x, but small dilatation ∂u j /∂ x j , and thus is far from zero. This means that the DF3 inflow condition is not a good approximation of the corresponding fully developed flow field at the inflow location. However, combined with the choice of u = 0, this approach produces a random incoming velocity field that has zero dilatation and thus small spurious acoustic waves.

The numerical database
Two reference flow cases have been selected, one representative of supersonic adiabatic boundary layers [22] and the other of hypersonic cooled boundary layers [51]. The former work refers to free-stream Mach number  The suffix RR denotes cases run with the recycling-rescaling procedure, and DF denotes cases run with the digital filtering procedure. The suffix L denotes benchmark simulations, run in long domains. M 0 is the free-stream Mach number, Re δ 0 = ρ 0 u 0 δ 0 /μ 0 is the Reynolds number based on the inflow boundary layer thickness, L x , L y , L z are the domain streamwise, wallnormal and spanwise sizes, x + , z + are the grid spacings in the wall-parallel direction, y + w is the minimum wall-normal grid spacing, and N x , N y , N z are the number of grid points, and δ f is the outflow boundary layer thickness M 0 = 2 and nominally adiabatic wall conditions (T w = T r = 1 + r (γ − 1)/2M 2 0 , where r = Pr 1/3 is the recovery factor and Pr = 0.72 the Prandtl number), and the latter has M 0 = 5.84, T w /T r = 0.25.
All simulations are carried out using a GPU-accelerated solver [2], which combines the energy-preserving properties of a sixth-order skew-symmetric central difference scheme [20] with the shock-capturing properties of a fifth-order weighted essentially non-oscillatory (WENO) scheme, through a modified Ducros shock sensor [21]. The sensor is disabled for the supersonic flow cases, whereas we have found that, although the flow does not include shocks, minimal numerical dissipation provided by WENO is needed for stability in the hypersonic flow cases. In all simulations, the diffusive fluxes are discretized with sixth-order central formulas, and time integration is carried out using a third-order Runge-Kutta method [46].
For both flow cases under scrutiny, two series of DNS have been carried out, one on relatively short domains, which serve to quantify effects of inflow seeding (RR-or DF-type), as compared to benchmark simulations, carried out in very long domains, which are verified to be yield to a healthy state of developed turbulence. The full list of DNS and the key computational parameters are reported in Table 2.
The supersonic data set includes six DNS in short domains and three DNS in long domains. All cases with RR seeding share the same mesh, which is identical to that used by Pirozzoli and Bernardini [22], whereas cases with DF feeding use a 50% longer mesh to reach comparable Reynolds number. In the short-domain simulations, the mean flow properties (as well as the turbulent stresses needed in DF) are taken from Pirozzoli and Bernardini [22]. The mean profiles for the long-domain DNS have instead been determined by applying the Van Driest transformation [36] to boundary layer profiles of the Musker family [19]. The hypersonic data set includes four DNS in short domains, and three in long domains. The choice of the inflow mean profiles is crucial in this case, as the Van Driest transformation is known to be inaccurate in the case of cold walls [40,51], and to perform poorly in the wake part of the wall layer, even under adiabatic conditions [44]. Hence, the inflow mean profiles for the M5.84-L1 RR simulation have been extracted from a coarser precursor DNS, and the simulation is used to provide mean inflow profiles for the other DNS. We note that the DF2 inflow condition is not considered for hypersonic cases as it fails to yield sustained turbulence. RR3 feeding is also disregarded as found to provide similar performance to RR1 in supersonic cases.
The statistical properties of the boundary layers are hereafter reported in terms of either standard or densityweighted (Favre) averages. Time averages are denoted asf , whereas Favre averages are defined asf = ρ f /ρ, with double primes denoting fluctuations thereof, f = f −f . For the short-domain simulations, time averages are accumulated over at least 800 and 1750 convective time units (δ 0 /u 0 ), for the M 0 = 2 and M 0 = 5.84 flow cases, respectively, taking advantage of spanwise averaging.
For the long-domain simulations, at least 1450 and 1750 convective time units have been used, respectively. Time averaging is initiated after a statistically homogeneous condition is reached, as inferred by monitoring the evolution of the spanwise-averaged wall properties. For all computations, spatially developing boundary layers, based on superposing van Driest-transformed velocity profiles with organized eddies resulting from the digital filtering procedure, are used as initial conditions.
For later guidance in the interpretation of the results, the distributions of the friction and momentum thickness Reynolds numbers for DNS in short domains are reported in Fig. 2, which makes it clear that

Results of benchmark simulations
A mandatory step for the evaluation of inflow turbulence feeding methods is the establishment of high-fidelity scaling laws for the main quantities of interest, as turbulence statistics, wall friction and heat transfer. This is of course important in its own sake for use in engineering analyses. In this respect, the first step is the identification of one or several diagnostic parameters, which allow to ascertain that an equilibrium turbulence state is achieved. The first diagnostic implies fulfilment of mean momentum balance [36], where σ xy = μ(∂u/∂ y + ∂v/∂x), and in which we have retained the normal Reynolds stresses u u and v v , the latter possibly being relevant at high Mach number [10].
Wall-normal integration of Eq. (11), combined with the continuity equation and the assumption of zero external pressure gradient, leads to the compressible equivalent of the von Kármán equation, originally derived for incompressible flow [27], where is the momentum boundary layer thickness and C f = 2τ w /(ρ 0 u 2 0 ) is the friction coefficient. Equation (12) also includes a term depending on the streamwise variation of the turbulent stresses. It can be assumed that an equilibrium state is reached if the residual of Eq. (12) is smaller than a certain prescribed tolerance [11,30,43]. We thus define a relative error in this metric as Two additional metrics refer to the development of the peaks of the turbulent shear stress, τ pk xy = max y (−ρ u v /τ w ), and of the streamwise velocity variance, τ pk x x = max y (ρ u u /τ w ). Correlations for the peak turbulent shear stress were suggested by Chen et al. [4] with B 1 = 13.7 for incompressible flow. As regards the streamwise velocity variance, Pirozzoli and Bernardini [23] found it grows logarithmically with Re τ , with A 2 = 3.35 and B 2 = 0.725. Two obvious additional diagnostic parameters include the friction coefficient (C f ) and the heat transfer coefficient, where q w = k∂ T /∂ y| w is the wall heat flux, with k = μ/(Pr C p ) the thermal conductivity. Power-law relationships versus Re θ are generally assumed for these parameters, namely where A 3 = 0.024 and B 3 = 0.25 for incompressible flow [35]. Finally, we consider the wall pressure variance, which is relevant in aero-vibroacoustic analysis. As shown by Jiménez et al. [12] and Pirozzoli and Bernardini [22], this quantity also includes contributions from distant eddies residing in the outer layer, hence incomplete boundary layer development may reflect into inaccurate prediction of this indicator. Farabee and Casarella [9] found a convenient representation in the form, with A 4 = −4.30, B 4 = 1.86 for incompressible flow. Figure 3 shows the distribution of the relative error in mean momentum balance, as defined in (13), for the three benchmark simulations in long domains. Large errors are found close to the inflow, similar to what reported by Schlatter et al. [30] for a incompressible boundary layers. In fact, the metric (13) provides a measure for the degree of development of a boundary layer in terms of its adherence to the thin-shear-layer equations [36]. Therefore, each neglected term in Eq. (12) can in principle be responsible of the observed imbalance. Momentum balance is satisfied to within 1% error only at a distance of about 100δ 0 . For further safety margin, we only consider the last third part of the domain to develop reference correlation laws for the various diagnostic quantities. The peak shear stress is considered next in Fig. 4a. The DNS data support validity of Eq. (14), with B 1 = 13.62, which is in remarkable agreement with what reported by Chen et al. [4]. Hence, it appears that density scaling is quite effective in mapping this particular metric from adiabatic boundary layers into their incompressible counterparts. The streamwise evolution of the peak streamwise velocity variance shown in Fig. 4 very well conforms with Eq. (15), with A 2 = 3.35, B 2 = 0.725 [23]. Figure 4 also reports the results obtained by Wenzel et al. [45], with corresponding best fits. Slight underprediction of both the pressure variance and the turbulent shear stress peak is generally found, as compared to the present data. Nonetheless, their results fall within a ± 2% uncertainty band from our inferred trends. On account of the different numerical approaches involved, the agreement is quite good.

Supersonic benchmark flow cases
Regarding the friction relation, Fig. 5a shows that the benchmark simulations follow Eq. (17)  23. Almost perfect agreement is observed between these correlations, and trends found by Wenzel et al. [45], throughout the Reynolds number range. Fig. 6, which shows longer inflow length as compared to the supersonic cases. In fact, about 35 − 70δ 0 are needed in order to fall within the 4% accuracy band. This finding is consistent with previous results of Huang et al. [11] and Wenzel et al. [43], who reported a monotonic increase in the inflow length with the free-stream Mach number. We additionally find that DNS with higher inflow Reynolds number require longer distance to adjust.

Fulfilment of the von Kármán integral equation for the benchmark hypersonic flow cases is verified in
The peak shear stress and the wall pressure variance are shown in Fig. 7a. As for the supersonic case, the shear stress peak still varies approximately as given in by Eq. (14), with B 1 = 29.8. The wall pressure variance still exhibits logarithmic increase with Re τ as predicted by Eq. (18), with A 4 = 3.37 and B 4 = 1.32, hence the growth is slower than in the supersonic case. The pressure variance is larger than in the supersonic case, as the result of finite heat transfer at the wall, as previously noticed by Huang et al. [11] and Zhang et al. [50].
The two most important flow properties in the study of hypersonic flow are obviously the friction coefficient and the wall heat transfer coefficient, which are reported in Fig. 8. A power-law dependence on Re θ is assumed as in equation (17), which yields the fitting constants A 3 = 0.0131, B 3 = 0.268, A 5 = 0.00774, B 5 = 0.272. Hence, the decay of C f is a bit steeper than in the supersonic case. All the fitting coefficients in equations (14)- (18) are listed for convenience in Table 3.

Effect of inflow turbulence seeding
Data from DNS in short domains are shown here, with the goal of evaluating the influence of the inflow seeding technique on the establishment of an equilibrium turbulence state, using the correlations established in the previous paragraphs as a benchmark.  Past this point, some scatter across curves is found, with RR generally performing better than DF, in that it yields smaller error for given distance from the inflow. No systematic differences are observed among different implementations of the two techniques. Figure 10 shows the distributions of the peak turbulent shear stress, which highlights some differences among the various approaches. Whereas the shear stress peak never exceeds unity in all DNS with RR, the DF1 and the DF2 cases exceed this threshold in a large region near the inflow, which makes the flow statistics quite unphysical. On the other hand, the DF3 implementation yields results more similar to RR, and quicker adjustment to equilibrium. The peak shear stress appears to monotonically approach unity at Re τ ≈ 400, but careful inspection of the figure suggests that DF tends to overestimate the trends returned from the benchmark DNS (dashed line), and to occasionally exceed unity.

Supersonic flow cases
The peak streamwise velocity variance is shown in Fig. 11. Notably, the DF has a more benign behaviour regarding this variable, as all DF implementations attain the correct behaviour at Re τ 300, despite large excursions from the reference trends in the initial transient. On the other hand, values in excess of Re τ = 400 are needed for RR statistics to fall within the ± 2% error band. As a result, RR requires longer fetch (x/δ 0 ≈ 80) for proper development of the streamwise turbulent stress.
The DF technique also appears to perform better than RR as regards the prediction of the friction coefficient, as shown in Fig. 12. In fact, despite much smoother initial transient, RR tends to consistently overpredict C f , with values falling within ± 2% error only at Re θ 3000, corresponding to the end of the computational domain employed for these DNS. The results of all RR-based DNS are quite similar, meaning that implementation details are unimportant in the observed behaviour. Despite substantial differences in the inflow recovery region, all DF implementations yield a C f within a ± 2% error band at Re θ 2000, i.e. at x/δ 0 50, and basically landing on the benchmark distribution at Re θ ≈ 3000. The DF3 implementation seems to be a bit slower than the other in this respect. It is noteworthy that the strong dips next to the inflow section seen in the DF simulations could in principle be mitigated as suggested by Larsson [16].
The von Kármán equation (12) can be used to connect observations made regarding the friction coefficient and the peak streamwise turbulent stress. In fact, Fig. 12 highlights irregular distribution of C f next to the inflow in the DF simulations, with sharp increase, followed by abrupt drop. The von Kármán equation clarifies that those regions are connected with u fluctuations, as the initial increase in the momentum flux deficit (namely dθ/dx) and decrease in wall shear stress compensate the growth of the streamwise stress. This mechanism is reversed in the ensuing region, featuring strong reduction in the streamwise stress an growth of the wall shear stress.
The distributions of the wall pressure variance are shown in Fig. 13. All RR implementations exhibit roughly similar behaviour, as they achieve a monotonically increasing trend compatible with theory [9,22] and with the benchmark distribution at Re τ 450 (namely x/δ 0 60). The behaviour of the DF-based DNS is richer as different trends are found depending on the implementation details. In fact, the baseline DF1 implementation largely overpredicts pressure fluctuations, throughout the computational domain. This large disagreement is partially cured through suppression of the streamwise velocity fluctuations in the DF2 implementation, and even more by making the inflow cross-stream velocity divergence-free, as in the DF3 implementation. Both DF2 and DF3 indeed exhibit reasonably correct trends at Re τ 300 and Re τ 400, respectively, which is in line if not better than for the RR cases. Still, some consistent overprediction of pressure disturbances in the developed region is observable, which we interpret as the result of acoustic disturbances arising at the inlet owing to the unrealistic structure of the inflow velocity fluctuations.  Figure 14 shows the relative error for the von Kármán equilibrium condition for the RR (left panel) and DF (right panel) techniques. The two RR implementations herein considered have quite similar behaviour, with ± 4% error band attained at x 50δ 0 , and the stricter threshold ± 2% is reached at x 60δ 0 . Regarding DF cases, the baseline DF1 implementation is slower to adjust towards zero error, whereas the DF3 implementation (suppression of u and stream function formulation) readjusts more quickly, falling within the ± 2% error band quite early (x 50δ 0 ). The peak turbulent shear stress is shown in Fig. 15, as a function of the friction Reynolds number. Unlike the supersonic cases, the peaks for both RR and DF cases remain strictly below unity, throughout the computational domain. Near equilibrium conditions are reached for Re τ 500 (corresponding to x ≈ 55δ 0 for RR1, RR2, DF1 and x ≈ 75δ 0 for DF3). At higher Re, the RR cases consistently follow the benchmark correlation, as the DF1 case also does. On the other hand the DF3 method seems to exhibit a plateau, followed by decline at Re τ 600 (corresponding to x ≈ 100δ 0 ). Figures 16 and 17 show the distributions of the friction coefficient and the heat transfer coefficient as a function of the momentum thickness Reynolds number. Not surprisingly, the two quantities exhibit a similar behaviour. We find that RR inflow feeding yields results which, past an initial dip, fall quite rapidly within the ± 2% error band. The dip is found to be much larger when DF is used. However, as also found for the shear stress, the baseline DF1 implementation adjusts quite quickly at Re θ ≈ 2000 similar to the RR cases, whereas the DF3 implementation requires longer fetch to reach equilibrium, at Re θ 3000 (corresponding to x 95δ 0 ). In Fig. 18, we finally consider the wall pressure variance. Regarding this parameter, the RR method is capable of attaining a monotonically increasing trend with the Reynolds number at Re τ 600 (corresponding to x 80δ 0 ); hence, the inflow adjustment length is much more than for supersonic cases. Past that location the data are in fair agreement with the The solid lines denote DNS data using RR inflow (a) and DF inflow (b); the dashed lines denote reference distributions obtained from the benchmark DNS, and a ± 2% error band benchmark correlation. The DF seeding exhibits results similar to the supersonic case. Whereas the baseline DF1 implementation starts from very large values and it tends to adjust to the expected behaviour towards the end of the computational domain, the use of solenoidal inflow fluctuations in the DF3 implementation yields a realistic behaviour starting much closer to the inflow (Re τ 350, corresponding to x 30δ 0 ). In any case, pressure fluctuations in the developed region appear to be a bit larger than expected, thus corroborating the notion that the DF technique is inherently more noisy than RR.

Discussion and conclusions
We next attempt to draw a quantitative comparison of the predictive capabilities of the various inflow turbulence seeding techniques in terms of achieving accurate representation of a state of fully developed turbulence. For that purpose, in Figs. 19 and 20 we show, for each diagnostic parameter herein identified, the inflow distance needed to target the benchmark correlations previously determined to a given accuracy. For the sake of the present analysis, the flow statistics shown in the previous sections have been smoothed by using a Savitzky-Golay causal filter.
Inspecting the results of supersonic DNS in Fig. 19, one may infer that DF turbulence seeding is capable of achieving faster adaptation of the friction coefficient than RR. In fact, all DF implementations achieve ± 2% confidence within 80δ 0 from the inflow, with DF2 substantially faster, with inflow length of about 30δ 0 . RR feeding does not seem to reach the same level of accuracy within the selected computational domain, the error on friction being still 2% error at 100δ 0 . The same error is obtained by all RR implementations, which also exhibit Solid lines denote DNS data using RR inflow (a) and DF inflow (b); the dashed lines denote reference distributions obtained from the benchmark DNS, and a ± 4% error band very similar trends. Regarding the peak streamwise velocity variance stress, one can likewise conclude that the DF1 and DF3 implementations are marginally better than RR. The DF2 case is obviously disadvantaged in this respect, as it has the lowest inflow turbulence kinetic energy. This criterion is, however, not very restrictive, as all RR and DF simulations require at most 40δ 0 to achieve ± 2% accuracy. The wall pressure variance exhibits quite a different behaviour. In fact, whereas RR shows consistent convergence towards the reference equilibrium value, DF in its variants consistently overpredicts this property. Improvements over baseline DF1 are obtained from formulations DF2 and DF3, which nevertheless require an inflow length of about 100δ 0 to achieve 5% accuracy. Proposed variants of the RR technique (mainly the RR2 implementation) yield small but observable improvement in this metric as may be argued upon close inspection of Fig. 13a. Reducing the intensity of streamwise turbulent fluctuation intensity of DF thus yields beneficial effects on the wall pressure variance, but on the other hand it yields poorer prediction of the peak turbulent shear stress. As a result, at least 60δ 0 are required to achieve 3% error in this quantity, whereas DF1 and all RR implementations require 30δ 0 for 2% tolerance. Baseline DF1, however, outperforms alternative implementations regarding the prediction of the peak shear stress, offering similar performance as RR in this respect. These conclusions need some adaptation for the hypersonic flow cases (see Fig. 20). In particular, we find that this time the baseline DF inflow feeding, while retaining minor advantage over baseline RR in the prediction of the friction and heat transfer coefficients, also yields similar predictions of the wall pressure variance and of the peak turbulent shear stress. Regarding the proposed modifications, we find that the RR2 implementation yields some advantage in making RR less noisy, while worsening the prediction of C h a bit. The beneficial effects of removing the streamwise velocity fluctuations in the DF3 implementation are instead lost, and instead we find significant deterioration in the prediction of friction and wall pressure variance. This is likely a result of very slow transition of the flow to a fully turbulent state, as suggested in Fig. 16b. Additional insight into the tendency of the flow statistics towards equilibrium can be gained by expressing the streamwise fetch from the inflow in terms of eddy turnover times, by exploiting Taylor's hypothesis [34] to define This is done for the error in the mean momentum balance in Fig. 21a, b, which highlights that all simulations collapse after four ETT, regardless of the inflow seeding. All supersonic cases fall within ± 4% error band after about one ETT, while the hypersonic cases after 1.5-2 ETT, in agreement with the literature results [11].
Regarding the friction coefficient, Fig. 21c, d shows that, whereas different inflow seeding yields different degree of accuracy in the initial transient, convergence to the benchmark distribution is basically complete after seven eddy turnover times for supersonic cases, and after about four turnover times in hypersonic cases, regardless of the seeding and of the inflow Reynolds number. As noted previously, the DF2 seeding yields faster adjustment of this particular parameter in the supersonic flow case, whereas baseline RR and DF seem to be most effective in the hypersonic flow case.
Recalling the goals of this study as originally stated, we believe that the main conclusions can be summarized as follows. First, we have established a standardized procedure to evaluate achievement of equilibrium conditions in numerical simulation of high-speed boundary layers. One of the conclusions in this respect is that no single criterion can be used to define the inflow length for arbitrary flow conditions, but rather each metric is subject to a different inflow length, which can also change as a result of the flow conditions. Perhaps surprisingly, we have found that the friction coefficient is particularly sensitive to inflow seeding, and it can bear memory of inflow seeding quite far from the inflow. An even more sensitive parameter is the wall pressure variance, which is typically overestimated from any inflow seeding technique, and which deserves special attention. Having each parameter to within a, say 1% error from true, does require extremely inflow lengths, certainly exceeding 100 inflow boundary layer thicknesses.
Surveying the current literature, it seems that very few studies have used such long domains, which raises questions about the reliability of reference data. For instance, we estimate that one of the commonly used databases [22], reports estimates of the friction coefficient which are about 1 − 2% too high. On the other hand, other studies as by Wenzel et al. [43] may be regarded as essentially devoid of spurious inflow effects. In this respect, DNS in very long domains as reported here can provide a more robust benchmark. In particular, we have developed simple correlations for several flow diagnostics which we believe can be confidently used as a benchmark in code development, as well as for the development of alternative inflow seeding techniques.
Regarding the latter subject, we have attempted to develop modifications of the baseline RR and DF techniques. The results in this respect are only partially conclusive, but perhaps can provide inspiration for follow-up studies. Regarding the RR method, we have studied simplifications of the baseline algorithm which facilitate practical implementation, and which seem to offer predictions very close to the baseline approach. Regarding the DF method, which may be more relevant for application to practical engineering problems, we have found that the baseline version suffers from large overestimation of the wall pressure variance, owing to large noise generated at the inflow. Whereas this issue can be partially cured by suppressing the streamwise velocity fluctuations in the supersonic case, the same approach is not equally effective in hypersonic flow cases, as a result of delay in numerical transition. A possible strategy to mitigate the drawbacks of suppressing u can be matching the turbulence kinetic energy content of the baseline DF implementation.
We expect that the present results can provide useful guidelines for the selection of inflow feeding techniques suitable for LES/DES of flow of industrial relevance, or in hybrid RANS/LES approaches [31,32,38]. Turbulence seeding is then used either at the inflow, or at the interface between scale-resolving and non-scaleresolving regions. In both cases a short inflow length is mandatory, not to interfere with physical phenomena of interest, e.g. shock wave impingement, flow control devices, wall curvature. The DF method is certainly a good candidate for the purpose as they can be adapted to complex geometries, even in the absence of homogeneous flow directions [1,42]. More importantly, the DF method can be easily incorporated into existing software. However, some weaknesses of DF should be overcome, including spurious pressure fluctuations from the inflow which can cause severe errors in aeroacoustic applications, and which can be minimized as shown here in the DF2-3 implementations. Finally, the computational overhead required by each inflow method must be taken into consideration as an important performance index. In this respect, we have found minimal difference in the wall clock time per iteration between DF and RR methods.
Overall, we find that the most important message from this paper is reiterating a caveat in the analysis of spatially developing flows, which despite careful numerical treatment require some inevitable adjustment length. In this respect, the error maps provided in Figs. 19 and 20 can provide useful guidance and feeling about the error bars to be expected for a given domain size, and for a given flow variable. We believe that all this information should be accounted for in the design of future large-scale DNS of high-speed wall-bounded flows, and in the development and assessment of novel inflow seeding techniques.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Conflict of interest
The authors report no conflict of interest.