DNS of a turbulent boundary layer using inflow conditions derived from 4D-PTV data

Unsteady, 3D particle tracking velocimetry (PTV) data are applied as an inlet boundary condition in a direct numerical simulation (DNS). The considered flow case is a zero pressure gradient (ZPG) turbulent boundary layer (TBL) flow over a flat plate. The study investigates the agreement between the experimentally measured flow field and its simulated counterpart with a hybrid 3D inlet region. The DNS field inherits a diminishing contribution from the experimental field within the 3D inlet region, after which it is free to spatially evolve. Since the measurement does not necessarily provide a spectrally complete description of the turbulent field, the spectral recovery of the flow field is analyzed as the TBL evolves. The study summarizes the pre-processing methodology used to bring the experimental data into a form usable by the DNS as well as the numerical method used for simulation. Spectral and mean flow analysis of the DNS results show that turbulent structures with a characteristic length on the order of one average tracer particle nearest neighbor radius r¯NN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bar{r}}_{\text {NN}}$$\end{document} or greater are well reproduced and stay correlated to the experimental field downstream of the hybrid inlet. For turbulent scales smaller than r¯NN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bar{r}}_{\text {NN}}$$\end{document}, where experimental data are sparse, a relatively quick redevelopment of previously unresolved turbulent energy is seen. The results of the study indicate applicability of the approach to future DNS studies in which specific upstream or far field boundary conditions (BCs) are required and may provide the utility of decreasing high initialization costs associated with conventional inlet BCs.


Introduction
The relative advantages and disadvantages of computational fluid dynamics (CFD) and experimental fluid dynamics (EFD) as they pertain to the study of turbulent flows are well understood. CFD, in particular DNS, delivers complete information about a turbulent flow field and is highly reproducible, however is limited in Reynolds number Re due to computational cost and may be sensitive to the choice of boundary conditions. By contrast, EFD can reach much higher Re but is typically constrained in spatial and/or temporal resolution and contains unavoidable measurement uncertainties whose correction may be ambiguous.
While the pros and cons of DNS and EFD are still very present, the simultaneous increase in availability of computing resources and progress in the development of experimental techniques in recent years have led to more overlap in the datasets that DNS and EFD produce in terms of Re regime and spatiotemporal resolution. For example, advances in post-processing algorithms for PTV image sequences and advanced PTV measurement setups (Schanz et al. 2013;Schröder et al. 2015;Schanz et al. 2016;Schröder et al. 2018) have allowed for significantly increased tracer seeding densities. As a result, high-frequency, transient 3D datasets have been produced which are the first of their kind and quality. The tracer particle localization scheme termed 'Shake-The-Box' (STB) approximates individual particles' trajectories as continuous B-splines through 3D space and time, allowing for the extrapolation of particle locations to new frames thus greatly reducing identification effort and improving accuracy. In doing so, the occurrence of so-called 'ghost particles' is reduced, in turn reducing measurement error. Ghost particles are spurious peaks falsely identified as tracer particles during PTV image post-processing which come about largely due to optical ambiguities of a given multi-camera arrangement (see Maas et al. 1993;Elsinga et al. 2006). Because the STB technique uses spatiotemporal (4D) information in its efficient predictor-corrector formulation, it is also referred to as 4D-PTV. As will be explored, such 4D-PTV datasets pose a unique opportunity to be used as BCs for DNS studies, which often struggle with the economical definition of realistic and/or appropriate BCs. The BC can be defined in various ways, for example as an 'initial value' type BC as in the present study or as a global immersed BC. The dual advantage of such a technique is the integration of realistic BCs in a numerical study and in turn the potential extraction of higher resolution information about the experimental problem through analysis of the numerical results. Combined experimental/numerical approaches are broadly referred to as 'hybrid methods. ' Hybrid methods refer to techniques seeking to leverage the relative advantages of EFD and CFD for combinations of various purposes, including noise reduction through assimilation to governing equations, the derivation of unknown field quantities, resolution enhancement, and the application of highly specific boundary conditions (Suzuki and Yamamoto 2015). Such goals are by no means separate, and the pursuit of one of these objectives often implicitly involves others. By assimilating real-world data within numerical simulation, highly resolved and reproducible outputs having agreement with the governing equations of fluid physics may be attained that have the added credibility of being at least partially rooted in natural observation. For example, techniques for the derivation of pressure and/or body forces (Fujisawa et al. 2005;Van Oudheusden et al. 2006;Murai et al. 2007;Schneiders et al. 2016) or temperature (Ma et al. 2007) have been successfully implemented in multiple studies. Similarly, modern numerical weather forecasting models (Skamarock et al. 2005) incorporate real-time sparse measurement data to define an initial condition (IC).
Measurement data may be integrated into a simulation through various means, for example via projection of proper orthogonal decomposition (POD) modes (Ma et al. 2003;Sirisup et al. 2004) or the coupling of PIV data to a 2D DNS via a proportional feedback term in the governing equations (Suzuki et al. 2006Yamagata et al. 2008;Suzuki et al. 2010). More recently, optimal Kalman filtering of the input data has allowed for the extension of the hybrid approach to higher Re (Suzuki 2012). The hybrid studies referenced employ advanced coupling schemes that are at least partially active within a complete subdomain of interest.

Methodology and objectives
In the present study, unsteady PTV data from a state-ofthe-art PTV measurement are used as an inflow boundary condition in a DNS. To evaluate the applicability of such a hybrid experimental-numerical approach, the canonical case of a flat plate TBL at a momentum thickness Reynolds number of Re ≈4000 is examined, where Re =̄ū e ∕̄e . Here, is the incompressible momentum thickness, the index e denotes values evaluated at the edge of the boundary layer and overlined variables (◻) represent time-averaged values. The external data are coupled via a relatively straightforward proportional forcing scheme (see Sect. 5.2) in a 3D inlet subdomain, downstream of which the flow is allowed to evolve freely. The current work is unique in that it uses 3D, unsteady PTV data from a state-of-the-art PTV technique rather than a single 2D plane, an achievement made realizable thanks to the availability of such datasets in the first place. Hybrid studies typically focus on the enhancement of experimental data within the complete area of interest, such as determination of the pressure field or the reduction of experimental noise. The current study focuses on the implementation and evaluation of a realistic turbulent inflow condition in which the incoming flow is already at a developed turbulent state, quantified by Re , though not spectrally complete. For reference, conventional DNS studies of flat plate TBLs often require significant computational effort to arrive at an equivalent Re . By qualitatively evaluating the resulting flow field's agreement with the input of the original EFD contribution, a first estimation of the potential of the methodology is assessed. Additionally, through analysis of the spectral development of the TBL within the DNS, inferences can be made about the degree of spectral deficiency inherent to the measurement and pre-processing chain. To the knowledge of the authors, the present study is the first of its kind in terms of resolution and dimensionality of inflow.
The present paper addresses the primary question of whether a 3D DNS using coupled 4D-PTV data is possible and worthwhile. Since the utility of the coupling as it pertains to the definition of an upstream BC is of central focus, the success of the study is evaluated through comparison of the freely evolving DNS field with the experimental field. This evaluation will reveal the net effect of factors such as potential mismatches in the exact far field condition, for example the freestream pressure gradient and presence of periodic span-wise conditions, on the evolution of the TBL flow. The central question of how the limited resolution of the external data source, in the present case meaning the finite tracer density, affects the development of turbulent kinetic energy k and the arrival upon a state of constant k production will be addressed. In short, the limiting factors of the methodology are to be investigated, such as resolution criteria.
The paper is organized according to the following structure. Section 2 provides a compact summary of the experiment from which the data for the simulation were derived. Section 3 characterizes the experimental data directly. Section 4 presents the interpolation method used to assimilate the Lagrangian PTV data to the rectilinear DNS grid. Section 5 presents the numerical method used, including information about the flow solver NS3D, the applied boundary conditions, as well as the hybrid coupling scheme used to introduce velocity perturbations into the flow domain. The results of the simulation itself are shown in Sect. 6, focusing on the mean flow (6.1) and wavenumber domain analysis (6.2). The simulation and experiment are then directly compared in Sect. 7, including snapshots of the instantaneous flow fields in Sect. 7.1, a quick demonstration of the application of spatial filtering in Sect. 7.2 and a look at the correlation between the experimental and DNS in Sect. 7.3. Finally, Sect. 8 summarizes the results presented in Sects. 6 and 7 and discusses the conclusions of the study.

Experimental data source
The 4D-PTV data were generated as a product of experiments described in Schanz et al. (2019) using the state-of-the-art 4D-PTV technique. The experiment was run in the Atmospheric Windtunnel Munich (AWM) at the University of the Federal Armed Forces Munich (Universität der Bundeswehr München) in Neubiberg, Bavaria. Figure 1 (from Schanz et al. 2019) provides a schematic of the experimental setup. In the 1.8 × 1.8 m test section, a double ramp wind tunnel model from the DLR VicToria Deutsches 2020) project was installed, which has an initial, upward-sloped favorable pressure gradient (FPG) region, followed by a ZPG flat region that is 4.0 m in stream-wise length, followed by a favorable pressure gradient (FPG) region that immediately precedes the adverse pressure gradient (APG) region associated with the 18° downward-sloped wall. The test section was seeded with Helium filled soap bubbles (HFSBs), which were illuminated by 10 LED array panels.Images were captured by a total of 12 high-speed cameras, 4 of which were focused primarily on the aft-most half of the ZPG region, having a resolution of 4096 × 2160 pixels per camera and capable of producing 1000 exposures per second. Datasets for 19 runs were received, each containing around 1360 timesteps, representing 1.36s of measurement duration at 1 kHz. The flow was seeded with HFSBs from two locations, using a total of 250 nozzles to produce the bubbles. To maximize the seeding density and prevent a large portion of bubbles from bursting upon contact with the wind tunnel meshes, 200 of the nozzles were placed immediately upstream of the test section, which was noted as having reduced the span-wise uniformity of the averaged boundary layer as seen in Fig. 2. However, since neither the original experiment nor the present numerical study rely on comparisons to canonical flat plate TBL studies in literature, but rather focus on the evaluation of new methodologies, this feature is inconsequential.

Characterization of experimental data
To allow for precise comparisons between the DNS and experimental results, detailed characterization of the original experimental flow field is necessary. This involves quantifying the degree of span-wise variance present in the 4D-PTV data as well as the isolation of a ZPG region having minimal absolute stream-wise (x) pressure gradient. A ZPG region is defined where the mean stream-wise acceleration is virtually unaffected by the FPG or APG regions, derived through analysis of the mean acceleration with respect to the stream-wise (x) and stream-normal (y, z) magnitudes. The ZPG box pictured in Fig. 3 represents the boundaries of the region found to be free of significant mean stream-wise pressure gradients as well as corner flow effects. It has dimensions [Δx×Δy×Δz] = 1.5⋅̄9 9,multi ⋅[5×2×4] , where 1.5⋅̄9 9,multi represents an approximate average boundary layer thickness at the middle of the wind tunnel model for several experimental runs, multiplied with a factor of 1.5. All later comparisons between the experimental and DNS flow fields are contained within this region, as illustrated by the superposition of the DNS domain within the ZPG region in Fig. 3. The precise determination of the dimensionless flow parameters needed for the DNS flow definition are derived by bin-averaging the tracers within this rectangular region. The average wall-normal profiles generated from the binaveraged tracer data are presented in Sect. 3.1. The calculation of the average tracer nearest neighbor distance r NN is discussed in Sect. 3.2. This characteristic length plays an important role in the interpretation of the analyses presented in Sect. 6 and may be thought of as the effective spatial resolution of the experimental dataset.

Spatial bin-averaging
Bin-averaging is performed to derive averaged flow metrics in a stationary 3D reference frame. Binning is the process of Wall-normal profiles of nearest-neighbor distance r NN quantizing or aggregating non-uniformly distributed discrete data into ranges. Here, binning is necessary for deriving a spatiotemporal average of PTV data in the form of discrete 'point cloud' style tracer measurements with sufficient statistical confidence. Spatial bin boundaries are defined, and tracer measurements for all times residing in each bin are averaged. Figure 3 shows the near-ZPG region as well as the stream-wise and span-wise binning schemes. Additionally, a superposition of the eventual DNS domain boundaries is shown. In the wall-normal direction, a bin width of 200 µm was used. Figure 2a shows the wall-normal ū + and u � u � + profiles for the [n x ×n z ] = [5×1] binning scheme, where it is evident that the profiles consistently collapse when nondimensionalized in wall coordinates. Primes (◻ � ) denote fluctuations with respect to the Reynolds averaged mean value. Figure 2b shows the same wall-normal profiles for the [n x ×n z ] = [1×4] binning scheme, showing a degree of span-wise non-uniformity present within the experimental data due to the presence of the upstream bubble rakes mentioned in (Schanz et al. 2019). For reference, Fig. 2 contains ZPG profiles from (Eitel-Amor et al. 2014) for comparable Re , represented by gray lines.Additionally, the u + =y + trend for the viscous sublayer, u + =1∕ ⋅ log(y + ) + C (with =0.41 , C=5.2 ) trend in the logarithmic layer and an approximation for the transition region from (Spalding 1961) are plotted.

Nearest neighbor calculation
In an effort to characterize the datasets received, one metric of interest is the nearest neighbor (NN) distance r NN of every tracer at every timestep. A k-dimensional tree based method was used to establish the distance between every combination of roughly 550K tracers per timestep at each of 1359 timesteps. The minimum of these distances, the nearestneighbor distance (or radius), is then saved as an additional scalar which can be averaged in time and space. Additionally, the standard deviation of r NN for tracers at all times in a given spatial bin delivers information about the uniformity of the tracer density. Figure 4 shows that, on average, the mean tracer distance r NN in the 'well-seeded' part of the experimental flow is roughly 5.3 mm.The plot colors indicate the x-bin (see Fig. 3), demonstrating the consistency of the seeding density in the stream-wise direction in addition to the wall-normal y direction between 100 and 2000 dimensionless viscous wall distance units y + =yu ∕ w . Additionally, the standard deviation of r NN is roughly 1.8 mm.
Close to the wall ( y + <100 ), both the mean r NN and its variance increase. Since no new tracers are 'randomly walking' from underneath the wall, the tracer replenishment is only coming from above (elsewhere from above and below). When tracer bubbles hit the wall they burst, so for a given random walk there are more sinks than sources of tracers. This phenomenon is well illustrated in Fig. 4, where the mean NN radius r NN starts to change at r NN from the wall, which is where the effect of the wall acting as a tracer 'sink' would be expected to first become visible on average, being that for y<r NN a sink exists within the local sphere defined by r NN . Also, the standard deviation NN radius has a sharp increase also at r NN =5.3 [mm] , which is likely due to tracer bubbles bursting. When this happens, the sudden elimination of one tracer from the local population causes a jump in r NN between timesteps, hence increasing the variance of a time-collapsed population.

Spatial interpolation to a rectilinear grid
The experimental data are in the form of sampled tracer trajectories, meaning that the velocity of the turbulent flow field is only known at scattered point locations for a given timestep. Since the coupled hybrid region acting as an unsteady inlet will require information at a set of stationary points (rectilinear grid points) a spatial interpolation of the data is necessary. Interpolation of sampled, limited-resolution data will generally introduce some level of error with respect to the continuous (and unknown) 'ground-truth' data, in this case the turbulent flow in the laboratory. The choice of interpolation method is therefore important and in the present case was chosen with two priorities in mind, namely the minimization of interpolation error and the fulfillment of known conservation relationships such as mass conservation. Naturally, these priorities are aligned in that reduction of error implies fulfillment of conservation laws (and viceversa). For this reason, significant effort was invested into the comparison and validation of potential methods.
One well-known interpolation method often employed for the general problem of interpolating 'point cloud' style data is inverse distance weighting (IDW), whereby the value at a given output point is determined by the distance-weighted average of N nearest input points. While highly scalable and computationally efficient, IDW was tested and seen to be highly sensitive to both free parameters: the choice of weighting exponent and number of NN points used for the average. Furthermore, IDW does not directly take into account any physical conservation law in its mathematical formulation.
Another class of commonly used interpolation methods for tasks of this type is radial basis function (RBF) based interpolation. Practical details of RBFs and their application to interpolation are available in (Schaback 2007; De Marchi and Perracchione 2018) and a formal mathematical description can be found in Buhmann (2000Buhmann ( , 2003. Similar to polynomial and the special case of piecewise polynomial (spline) interpolation, RBF interpolation involves setting up and solving a system of linear equation in the classic form A =f , such that a coefficient vector =A −1 f is obtained and used to form a continuous interpolant through linear dot-product summation of the coefficients and a set of basis functions (e.g., [1, x, x 2 , x 3 , … ] in 1D polynomial interpolation). RBF interpolation differs from classical polynomial interpolation in that the interpolant is constructed from the set of distances between the known point coordinates rather than the spatial coordinates themselves. This has enormous advantages when extending the method to higher dimensions, where polynomial interpolation quickly becomes expensive, difficult, and is not guaranteed to produce an interpolant for any arbitrary set of points when using a given set of basis functions (Mairhuber 1956). Discussion on the differences between polynomial interpolation and RBF interpolation, especially concerning interpolation in multiple dimensions can be found in (Buhmann and Jäger 2019). Section 4.2 gives a summary of radial basis function interpolation for 3D scalar fields as well as the extension to a 'divergencefree' formulation for 3D vector fields.

Summary of radial basis function interpolation
In its conventional formulation, RBF interpolation seeks to derive a continuous approximation (or interpolant) function s of a scalar-valued function f: s( )≈f ( ) . Here, the scalarvalued function f is assumed to be a function of (and is interpolated in) 3D space ℝ 3 , though the method is equivalent for functions of any space ℝ n since it transforms the support point coordinates to a set of shifted scalar distances. The function f is therefore a multivariate function of three coordinate dimensions, or equivalently the coordinate vector =[x, y, z] . If f is known at N known points in space, the condition is set on s that the input f must be exactly recovered at these N points, i.e., s( i )≡f ( i ) for i=1, 2, … , N . The interpolation process has two main steps: the solution of the scaling vector and the formation and evaluation of the interpolant. In the first step, a linear system of equations (LSE) of the form A =f is set up: where the A matrix is the square ( N×N ) distance matrix with every Euclidean distance r=‖ i − j ‖ element evaluated by a scalar-valued, positive-definite RBF . Typical candidates for include =r , r 3 , r 2 log(r) or e − r , where is an adjustment parameter. The f vector contains the values of f ( i ) at the N known points. The vector is the unknown scaling coefficient vector, which can be solved for through inversion of A as =A −1 f . In the second step, the interpolant function s( ) is built and evaluated in the form s=B =BA −1 f , where the B matrix is an M × N matrix similar to A, but is now the RBF evaluation of the distance matrix between the N known and M 'new' points. The s vector is the resulting vector of interpolated values at the M new points. In the present case, the number of tracers at a given timestep within and in the vicinity of the inlet region fluctuates around the average value of N≈6000 . The tracer data are used to form the interpolant, which interpolates the PTV flow field to M=[33×424×1280]=17.9M grid point locations for use in the DNS.
In the present study, the interpolated scalar f would represent each individual component of the vector-valued velocity field =[u, v, w] T and the 'new' points would be the inlet grid points. Interpolating each velocity component independently would not leverage known physical relationships between these quantities, such as being approximately divergence-free ∇⋅ ≈0 due to mass continuity and the effectively constant mass density at the low experimental freestream Mach number of M ∞,exp =0.019 . Therefore, a vector field ( ) can be interpolated using a matrix-valued RBF kernel (Narcowich and Ward 1994). A significant advantage of this formulation is that the RBF kernel can be formulated in such a way that the resulting interpolated vector field inherits the quality of being divergence-free by virtue of how the RBF kernel is constructed. Formal mathematical foundations and analysis of this method are found in Lowitzsch (2002Lowitzsch ( , 2005 and further applied examples are presented in McNally (2011); Yang et al. (2014). This extension reduces interpolation error by eliminating ∇⋅ ≠0 error. Alternatively considered, the accuracy of the interpolant is increased because it now has more detailed constraints.
A vector-valued approximation function ( ) of a vector-valued field ( ) is now once again sought such that ( )≈ ( ) . Drawing on Helmholtz's theorem, it is known that any smooth, well enough behaved vector field can be decomposed as the sum of a curl-free (irrotational) and a divergence-free (solenoidal) field =∇p+∇× , where the gradient of the scalar potential ∇p is curl-free by identity ) and the curl of the vector potential ∇× is divergencefree by identity ( ∇⋅(∇× )=0 ). In the special case that the vector field is known a priori to be divergence-free ( ∇⋅ ≡0 ), then it can be shown that Helmholtz's theorem reduces to =∇× (Griffiths 2014), i.e., that is the curl of a vector potential field . The vector function can still not be uniquely defined as the curl of a vector field , i.e., ≡∇× because for any , where ∇× = , the offset of by the gradient of any arbitrary scalar potential field ∇p � would also yield a field where ∇×( +∇p � ) is also = (see (Jackson 2002)). This non-uniqueness is summarized by the concept 'gauge freedom' or 'gauge invariance,' i.e., that inconsequential additional degrees of freedom (realized by 'gauge transformations') are available which yield equivalent observed states of a system (Jackson and Okun 2001). To break the gauge invariance and hence the non-uniqueness problem, additional conditions are needed, i.e., the gauge must be 'fixed.' One common choice of gauge is the Coulomb gauge, equivalent to saying that the vector potential is itself divergence-free, i.e., ∇⋅ =0 , which can be implemented by setting =∇× (see (Jackson 2002;Stewart 2003)). Thus can be uniquely defined as ≡∇×(∇× )=∇(∇⋅ )−∇ 2 . This is the starting point for the definition of the divergence-free approximation function , which can be expanded in 3D Cartesian coordinates: The 1D derivative operator x i is shortened to x i above for compactness. The operator {−Δ +∇∇ T } maps the field to and is also used in the construction of the matrix-valued RBF with being the 3x3 identity matrix, Δ being the Laplace operator Δ=∇ 2 =∇⋅∇ and ∇ is the del operator being ∇=[ x , y , z ] T in 3D Cartesian coordinates. A scalar RBF with sufficiently high order can be modified by the differential operator {−Δ +∇∇ T } (see (Lowitzsch 2002)) to form the matrix-valued RBF ={−Δ +∇∇ T } . (3) In the present study, a Wendland polynomial (Wendland 1995) was chosen for , specifically 4,2 , which in 1D is where 4,2 (r) is defined as null for r>1 . The 3D Cartesian derivatives of 4,2 ( ) needed for the construction of are as follows: in which the components r x , r y and r z are 1D distance components of =[r x , r y , r z ] T and ̂ is the norm ̂ =‖ ‖ , i.e., the 3D Euclidean distance. The same linear system is set up as in the scalar case to find the coefficient vector , except now in block-structured form whereby A now has dimensions 3N×3N, and and f have dimensions 3N× 1. Each block of A, i.e., A ij is the Fig. 5 Schematic 1D illustration of interpolation scheme used to enforce no-slip condition at the wall corresponding Ψ ij evaluation of with dimensions N × N. The blocks of the vector f are the stacked vector field components, i.e., f =[u 1 … u N |v 1 … v N |w 1 … w N ] T in the present study.
The evaluation of the interpolant B =BA −1 f is analogous to the scalar case except now the block-structured B has dimensions 3M× 3N and each M × N block B ij is the corresponding Ψ ij evaluation of , where is the vector of 1D distances between the known and new points. When the interpolant B is evaluated, the resulting vector field (formed by reshaping f from 3M× 1 to M × 3) is divergencefree. Formal mathematical derivation of this feature is shown in Lowitzsch (2002) and additional discussion can be found in McNally (2011) and Yang et al. (2014). The divergencefree feature of the output vector field becomes conceptually clear when considering that [u, v, w] at every output point is realized as the linear dot-product sum of weighted basis functions. Since the basis function kernel is constructed as the curl of the curl of a single polynomial function, the divergence of their resulting collocated sum is automatically null by identity ∇⋅(∇×◻) = 0.
The divergence-free RBF interpolation method has been implemented and its increased accuracy with respect to the scalar formulation has been confirmed ( Appelbaum 2020). The interpolation of experimental velocity field to a grid block is the first step in the preparation of the unsteady inlet for the DNS.

Preparation of the unsteady velocity inlet
Following the spatial interpolation of the tracer data to the DNS inlet grid, a series of additional pre-processing steps are necessary to make the data suitable for numerical simulation. These include treatment at the wall boundary as well as periodic boundaries, the solution of the pressure and temperature fields, and finally the non-dimensionalization and subsequent Mach number scaling of the input data.
A major challenge is posed by the finite tracer density near the wall boundary ( y = 0 ) when the dimensionless wall-distance is considered. The dimensional tracer density is relatively constant, with a slight reduction near the wall due to the tracer bubbles bursting (see Fig. 4), however when considering the local tracer density relative to velocity gradient ∇ (or, similarly, spatial tracer density normalized by the local dimensionless wall distance) one sees that near the wall there relatively sparse tracer coverage. For reference, the average nearest neighbor distance r NN ≈5.3[mm] (see Sect. 3.2), is equivalent to a wall distance y + ≈100 , well into the logarithmic layer. This means that, if reduced to a 1D problem, there is on average one single tracer available to resolve the curvature of the velocity profile, including the no-slip condition, below y + ≈100 during the interpolation process, which must produce output for every grid point at every timestep.The no-slip condition is therefore not universally respected during the interpolation step (see Sect. 4.1), so measures must be taken to force y=0 ≡ to avoid possible discontinuities in the computational domain. For example, one can introduce 'extra' artificial tracers having null velocity at the wall grid points for every timestep during interpolation, however this can have the opposite effect in which the null wall velocity becomes over-represented spatially if no real tracers are momentarily in the vicinity to otherwise influence the result, meaning that the velocity profile approaches zero too sharply too far away from the wall. Instead, a wall blending scheme is applied in a post-processing step to force the velocity of the experimental domain to y=0 ≡ at the wall. Figure 5 schematically demonstrates this procedure and shows the true proportion between r NN and ̄9 9 in the experimental data. A third-order natural spline (blue) is lofted and blended via a window function to u(y) (red), replacing the lower portion of the velocity profile (red dotted). The black line represents a reference average velocity profile ū(y) at the experimental Re . This procedure is performed for all velocity components u,v and w. A blending function is also used to transition the unsteady field to the steady mean flow field at the span-wise boundaries of the inlet to enforce the periodic condition. Such boundary condition enforcement schemes introduce a certain amount of non-physicality to the incoming flow, the effects of which are discussed in Sect. 6.
In a final step, the flow field variables are non-dimensionalized according to the convention required by NS3D and split up into separate MPI domains for use in the simulation. The dimensionless problem is defined by the characteristic dimensionless quantities where Rē9 9 is the hydrodynamic boundary layer thickness Reynolds number, Pr is the Prandtl number, M ∞ is the freestream Mach number, k ∞ is the freestream thermal conductivity, ∞ is the freestream density, = c p c v is the specific heat ratio and c ∞ is the freestream speed of sound. The freestream velocity magnitude is signified by U ∞ and generally U indicates the velocity magnitude, i.e., U=‖ ‖=̂ =(u 2 +v 2 +w 2 ) 1 2 . The freestream dynamic viscosity ∞ is determined by Sutherland's Law (Sutherland 1893). The length, velocity and time scales are non-dimensionalized by (8) where ̄9 9 is the average inlet hydrodynamic boundary layer thickness. The 3D, unsteady dimensionless temperature, pressure and density fields are calculated as deviations with respect to 1D, dimensionless averaged wall-normal y profiles, which are calculated as follows. The averaged 1D dimensional temperature profile T (y) is calculated using the Crocco-Busemann relation (Crocco 1932;Busemann 1935;White and Corfield 2006) where T aw is the dimensional adiabatic wall temperature and r is the recovery factor r=Pr 1 3 (Stalder et al. 1950;Ackermann 1942;White and Corfield 2006) with Pr=0.71 for air. The low Mach number of the experiment ( M ∞,exp = 0.019 ) means that stagnation heating and compressibility effects in the measured flow field are negligible; therefore, the necessity of deriving the unsteady temperature and density fields is purely numerical and a product of the Mach scaling applied to the dimensionless simulation, as described in Sect. 5.3. With an assumed constant pressure p over the boundary layer, the density profile is then derived using the equation of state. The inputs to the equations above are from the full (interpolated) transient field which is averaged in the stream-wise (x) and span-wise (z) directions as well as in time so that only a 1D wall-normal profile remains. The 3D unsteady dimensionless temperature, pressure and density fields are then calculated by their unsteady dimensionless, mean-removed (◻ � ) components as where the 3D primitive variable field is at all points compared to the corresponding wall-normal (y) position of the averaged profile, indicated by the overlined variables (◻) . More details regarding the pre-processing toolchain can be found in Appelbaum (2020).

Flow solver: NS3D
The flow solver NS3D is used to simulate a ZPGTBL.

Computational domain and boundary conditions
The computational grid is composed of The bottom boundary of the DNS domain is designated as an adiabatic, no-slip wall such that ( T∕ y) y=0 ≡0 , y=0 ≡ and ( p∕ y) y=0 ≡0 . The downstream boundary is defined as a subsonic outflow, which requires a steady flow field ( ∕ t) N =( ∕ t) N−1 ≡0 and thus the absence of any streamwise velocity gradient ( ∕ x)≡0 , assured through the presence of sponge zones. Sponge zones numerically limit the flow field's deviation from a reference field. The time derivative of the vector of conservative numerical fluxes Q= [ , u, v, w, E] T is adapted by the relation (12) where G( ) is a gain parameter. Fluctuations are attenuated by means of an additional source term in the discretized Navier-Stokes (NS) equations which is proportional to the amplitude of the fluctuation of each element of Q . The gain parameter is a constant multiplicative factor on the attenuative body force. The reference field Q ref is typically the so-called 'baseflow' or steady 1D y-normal mean profile (projected in [x, z]) which serves as an initial condition and steady boundary condition where required. The steady baseflow profile is calculated from the averaged experimental flow field, as described in Sect. 4.3. In the present case, an additional unsteady disturbance field is present in the inlet region (the green box in Fig. 3), serving locally as Q ref . A unity gain parameter G=1 provides a relatively 'soft' damping of disturbances with respect to the reference field and is therefore used at the downstream, top and lateral boundaries. A null gain parameter G=0 means that no artificial damping of the simulated field is present and is therefore used in the bulk of the simulation domain. Within the inlet region, the reference field Q ref is the unsteady field derived from the 4D-PTV data, as described in Sect. 4.3. Here, a relatively high gain parameter 0<G<50 is used to enforce the externally derived flow field by strongly damping the solution vector's deviation from the reference field. Figure 6 (13) illustrates the distribution of the sponge gain G in the DNS domain.By x * ≈0.16 downstream of the inlet, the sponge gain is null and the DNS solution is thus no longer (directly) forced by the inlet field. The lateral ( −z,+z ) boundaries are defined as periodic. The top boundary is a characteristic freestream boundary, which minimizes the reflection of acoustic waves. Oblique waves are minimized through the transitional grading of G in the boundary-normal direction, as well as through grid stretching in which the growth rate of the grid is increased in a given dimension. A more detailed description of the boundary conditions used in the present case may be found in (Babucke 2009).
The domain is spatially discretized with a rectilinear grid, defined through the projection of three orthogonal 1D coordinate vector arrays. In the present case, the grid is uniformly spaced in the span-wise (z) dimension. In the wall-normal (y) direction, the grid is stretched in the direction of the top boundary ( +y ) from a minimal cell size Δy=y i+1 −y i adjacent to the bottom wall. The grid growth rate Δy i+1 ∕Δy i in the first block of 136 cells near the wall is limited to 1.016 (1.6%), after which it is decreased to 1.00167 (0.167%) until about y∕ 99 ≈1.5 , where it is allowed to grow at a rate of 1.008 (0.8%) up to the top boundary. In the stream-wise (x) direction, the grid has a uniform spacing until x * =x∕̄9 9 =3 , after which it is stretched in the direction of the downstream boundary, promoting numerical damping of turbulent structures as they approach the subsonic outflow boundary. Starred coordinate variables (◻ * ) indicate the spatial coordinate non-dimensionalized by the average hydrodynamic boundary layer thickness ̄9 9 of the experimental data at the inlet plane location and 'plus notation' (◻ + ) specifies that a length has been non-dimensionalized by the local viscous length scale w ∕u . The grid resolution has been confirmed to be sufficient through post-processing of the DNS results, which has revealed that the first cell height Δy 0 is sufficiently small such that Δy + 0 =(Δy 0 u )∕ w stays less than roughly 0.60 throughout the complete area of interest ( x * <3 ). In the stream-wise direction, the maximum Δx + is 12 and the maximum Δz + is 3.6. The sampling period of the experiment (0.001[s]) is roughly 413 times larger than the relatively small timestep size required for the explicit time integration scheme of NS3D. Accordingly, a temporal interpolation scheme is needed to provide data at every sub-timestep of the Runge-Kutta temporal integration scheme, for which a thirdorder Lagrange polynomial based interpolation is used.

Flow parameters
The dimensional and non-dimensional parameters of the experiment and DNS are listed in Table 1.The defined inlet hydrodynamic boundary layer thickness Reynolds number Rē9 9 ≈43e3 is kept constant, however the freestream Mach number M ∞ of the dimensionless problem is increased to M ∞ =0.85 to reduce the simulation cost.Compared to an equivalent simulation at M ∞ =0.3 , the computational effort at M ∞ =0.85 is reduced by a factor of roughly three due to the less strict stability constraints on the timestep Δt max related to the discretization of the energy equation, see Babucke (2009). Error introduced by compressibility effects at M ∞ =0.85 is predicted from similar studies to not be particularly significant (see Wenzel et al. 2018), and in the present case is likely to fall beneath the magnitude of error introduced during grid interpolation or boundary blending. The dimensionless simulation result is re-dimensionalized using the experimental M ∞ for dimensional comparisons. Quantification of the impact of M ∞ scaling remains to be evaluated in further studies, however the already high agreement between the experimental and DNS flow fields presented in Sect. 7 despite such potential error is encouraging for future studies with lower M ∞,DNS ∕M ∞,exp .

Performance and output
The NEC SX-Aurora TSUBASA vector platform at the High Performance Computing Center Stuttgart (HLRS) was used to run the simulation. Four vector hosts, each containing eight vector engines, were used to perform a total of 540K timesteps, corresponding to 96.1% of the total possible duration of the 1.36 s of experimental measurement duration. Roughly 175 h of wall time were needed for completion.

Simulation results and analysis
To provide a first impression of the resulting DNS flow field, Fig. 7 presents a typical snapshot of the DNS domain fed by the unsteady velocity inlet, illustrating the normalized stream-wise velocity u∕U ∞ (a) and the normalized, mean- normal plane. Note that only the bottom 15% of the BL thickness is pictured. The contour plots represent one fourth the resolution of the simulation. Large turbulent structures resolved by the experiment are introduced in the inlet region (the left-most gray lines), downstream of which the structures break up into finer characteristic scales. As a result, the characteristic turbulent structure size at the inlet tends to be larger than downstream. The 'breakup' of the near-wall structures introduced to the domain tends to fluctuate in stream-wise (x) location, however is typically initiated by x * ≈2.0 at the latest. Additionally, non-physicalities resulting from the interpolation scheme adopted to force a no-slip condition (Sect. 4.3) are (at least intermittently) present. Such structures in the near wall region do not necessarily fulfill the discretized transport equations and boundary conditions of the numerical method used in the DNS and therefore are quickly dissipated within the simulation, leading to a transient effect immediately downstream of the inlet. Between the first and second dotted gray line, the effective forcing of the unsteady inlet field onto the flow solution is reduced to zero, controlled by the sponge parameter (see Fig. 6). After the second gray line, the structures freely evolve with no artificial forcing.
In the following sections, the simulation results will be analyzed in more detail to better characterize and quantify the behavior of the flow when using the present hybrid seeding approach, however the general trend captured by Fig. 7 remains present. As will be shown, the large flow patterns resolved by the experiment are preserved in the simulation, and the unresolved high-frequency content absent in the experimental results due to limited tracer density is relatively quickly recovered, meeting the goals of the present initial study and showing promise for future studies. The flow's evolution will be characterized with the help of analysis of the mean flow in Sect. 6.1 and wavenumber domain analysis in Sect. 6.2. Afterward, the experimental and simulated flow fields are directly compared in Sect. 7. The DNS flow field in Sect. 7.1 is shown alongside a filtered version of the flow field which is discussed in Sect. 7.2, though it should be emphasized that the filtering is simply a postprocessing exercise and not a separate numerical simulation.

Mean flow development
The flow field at the DNS inlet is reconstructed from 4D-PTV data and lacks some unsteady high-wavenumber content of the laboratory flow. As such, the reconstructed inlet flow field can be analogously described as a spectrally low-pass filtered version of a turbulent flow, as will be demonstrated by filtering DNS results in a post-processing exercise in Sect. 7.2. As the coarsely resolved TBL flow at the domain inlet is introduced to the DNS domain, it undergoes an adjustment process before a fully developed turbulent kinetic energy k spectrum is once again reached (see Fig. 8). As previously mentioned, there are two main drivers for the adjustment effects seen immediately downstream of the inlet: the limited tracer density of the experiment, leading to a lack of high frequency content in the unsteady inflow, and the side effects of boundary blending at the wall.
Non-physical flow structures that don't fulfill the discretized transport equations are quickly numerically dissipated, as evidenced in Fig. 8b by the immediate reduction in the root-mean-square (RMS) dimensionless wall shear stress + w,rms , as well as in (a), where a relatively thick layer of low dimensionless turbulent kinetic energy k + is present near the wall. The dissipation of such structures and relaxation of shear stresses produces a relaminarizing tendency, as seen in reduction of the skin friction coefficient c f =2̄w∕( ∞ U 2 ∞ ) whose minimum occurs at x * ≈0.7 , downstream of the minimum in + w,rms at x * ≈0.3 . Due to its non-dimensional representation in wall units, the relaminarizing tendency causes k + to increased, having a peak 'bubble' at (x * ,y * )≈(1.0,0.02) . Subsequently, the off-wall fluctuations mix back into the near wall region to produce a peak in + w,rms around x * ≈1.2 . By 2.5≲x * ≲2.0 , the transient process begins to reach a more or less converged state. Interestingly, the normalized hydrodynamic boundary layer thickness 99 ∕ 99,exp grows relatively constantly throughout the domain (Fig. 8b) in spite of this transient process, indicating that the re-adjustment of non-physical conditions at the wall has little effect on the BL thickness.
Such transient effects are also visible in the progression of the various Reynolds numbers in Fig. 8c. The momentum thickness Reynolds number Re =(U ∞ )∕ ∞ initially increases as non-physical structures dissipate, then reduces slightly during the pseudo-relaminarization phase until x * ≈0.7 , consistent with the minimum in c f in (b). After x * ≳1.0 , where turbulent energy production begins to recover, Re once again increases as expected for a TBL. The progression skin-friction Reynolds number Re =(u 99 )∕ w closely resembles that of the skin-friction coefficient c f in (b), seeing as both the friction velocity u = √ w ∕ w and c f are proportional to the wall shear stress w = ( u∕ y) y=0 . The shape factor H 12 = * ∕ relates the displacement thickness * and the momentum thickness and indicates the state of turbulent development of a ZPGTBL. Laminar BLs have higher shape factors, meaning that more deficit in u(y) is present with respect to the freestream U ∞ over a larger range of the boundary layer thickness, and accordingly have less severe wall velocity gradients ( u∕ y) y=0 . The slight increase in H 12 in the range 0.3≲x * ≲0.7 is therefore consistent with the previously mentioned relaminarizing tendency.

Spectral analysis
The turbulent character of the DNS flow field can be efficiently summarized in a composite sense by inspecting the flow in the frequency domain. In the present case, spectral analysis will help to answer a primary question of interest, namely how quickly the recovery of unresolved turbulent scales occurs. Figure 9 presents a wavenumber power spectral density (PSD) of the u ′ u ′ component of the Re-stress tensor, where k ,99 is the normalized linear wavenumber k ,99 =k̄9 9 ∕(2 ) =̄9 9 ∕ . The 1D span-wise wavenumber spectral analysis is performed on thin volumes at seven  Figure 9a presents the PSD for every available wall-normal position vs normalized wavenumber. As the flow propagates through the DNS domain, the fine scales which are absent at the inlet begin to re-develop, as made visible by the extension of the high power region into the high wavenumber (right) part of the plot with increasing x. Eventually, between 2.0≲x * ≲2.5 the shape of the PSD in the y + vs k space converges to a final form, indicating that the flow has reached a developed turbulent state. Figure 9b provides a more quantitative representation of this convergent effect. Here, the spectrum at every selected x-position is plotted in the PSD vs k space for a number of y + locations, increasing in distance from the wall. At all chosen y + locations, the spectra show a convergent trend in the stream-wise direction. It's also clear that the net gap in energy density is not equal at all wall distances, quantified by the size of the curve envelope. Since the characteristic turbulent structure size is finer approaching the wall, a higher proportion of the spectral energy is concentrated in the high wavenumber regime closer to the wall. It therefore follows that spectral energy gap is wider in this region, as information about high wavenumber modes is limited relative to low wavenumber modes during the construction of the inlet. However, since the finer structures closer to the wall also develop over a shorter characteristic period, the larger spectral gap is closed more quickly, i.e., spectral Fig. 11 Correlation coefficient of u ′ between experiment and hybrid DNS result recovery occurs within a shorter transit distance. Closure of the tighter energy density gap in the outer log layer occurs over a longer distance, as less stress is present to drive the development process. Additionally, the flow speed is higher, so there is less residence time between the fixed stream-wise locations in which the development processes can occur. Nevertheless, by roughly 2.5⋅̄9 9 convergent arrival upon a power spectral density profile can be seen for all wall distan ces.

Instantaneous flow
The instantaneous flow fields of the experiment and DNS can be compared to help qualitatively visualize and support the findings of the spectral analysis in Sect. 6.2. Figure 10 depicts contour plots of the normalized stream-wise velocity u∕U ∞ on wall-normal slices and wall-parallel slices in the left column, as well as the difference between the DNS and experiment Δu∕U ∞ in the right column. Also present in Fig. 10 is a filtered version of the DNS result, which will be addressed in Sect. 7.2. Qualitative large scale structural similarity can be clearly seen when the interpolated experimental measurement and DNS domains are compared in Fig. 10, indicating that the goal of reproducing the main characteristics of the instantaneous turbulent field has been met within the considered domain. This result is consistent with the minimal changes seen in the PSD below a cutoff wavenumber of k ,99 =̄9 9 ∕(3⋅r NN ) in Fig. 9. The difference in normalized stream-wise velocity u∕U ∞ between the DNS result and the interpolated experimental result as seen in the right column in Fig. 10 indicates increasing disparity with increasing stream-wise position x and in the near-wall region. Difference within the inlet region itself is due to the temporal interpolation scheme. The large scale similarity and the ability to identify equivalent structures in both datasets is a highly encouraging result, indicating that at least the large scale, most energetic structures near the TBL edge may be accurately introduced and propagated through a numerical simulation domain of choice. In doing so, an already well developed, realistic boundary condition at the TBL edge is present, fostering the TBL's quick recovery to a statistically converged state.

Spatial low-pass filtering
A spatial low-pass convolution filter may be applied to the DNS results in an attempt to derive an analogy for the net attenuative characteristic of the experimental procedure and pre-processing toolchain, which are inherently limited in resolution by the tracer density.The aim of such an exercise is to show that the experimental field and its subsequent interpolation to the DNS inlet grid is broadly equivalent to a low-pass filtered DNS result. A simple Gaussian-type filter kernel was chosen, leaving the kernel width as a variable parameter. A snapshot of the lowpass filtered DNS result is presented in Fig. 10. The filter operation was run on an x-normal plane at x * =3.0 , where the DNS solution is appropriately turbulently developed. Figure 12 shows the effect of low-pass filtering on the span-wise wavenumber spectrum, namely that a simple Gaussian 'blur' filter with a kernel width ∕r NN = 0.4 applied to a statistically converged DNS result yields a flow that looks qualitatively similar to the interpolated experimental result both in the spatiotemporal (Fig. 10) and wavenumber domains.

Correlation between experimental and DNS flow fields
Wavenumber domain analysis gives composite information about how the energy density of turbulent fluctuations is distributed across wavenumber modes, however it does not necessarily deliver information about the agreement between the instantaneous flow fields. However, the correlation coefficient between the experimental and DNS flow fields quantifies their normalized covariance and is a metric for their degree of variational similarity in a time-dependent manner. Figure 11 shows the Pearson correlation coefficient r (Pearson 1901) on a span-wise normal slice where for every grid point in [x,y], r is calculated between two time series of velocity magnitude ̂ =(u 2 +v 2 +w 2 ) 1 2 . In Fig. 11a, ̂ (t) from the experimental measurement is correlated with that of the filtered DNS result (shown in Fig. 10 and 12) and 11b shows r(̂ ) between the experiment and the unfiltered DNS result. The time series of the mean-removed, normalized velocity magnitude ̂ � (t)∕U ∞ for two points is plotted in (c) and (d). One point is located near the wall where the DNS flow quickly decorrelates from the experimental flow, and the other point is farther from the wall where the flows remain highly correlated. The difference in correlation at the two locations is plainly visible when comparing (c) and (d), though the slightly better correlation of the filtered DNS field is fairly subtle.
In the near wall region, the tracer density is low compared to the local strain. As discussed in Sect. 4.3, this low relative information density necessitates a post-treatment of the interpolated experimental field to ensure that the no-slip condition is restored at the near wall region of the unsteady velocity inlet before simulation. In doing so, artifacts stemming from this treatment are introduced in the region y≲1⋅r NN or y + ≲100 , for example an unrealistically steady near-wall profile in which small-scale motions are not resolved. Additionally, the near wall region sees the highest production of turbulent kinetic energy k and the highest frequency of turbulent motion, serving to quickly amplify the deviation from any initially well-correlated state. For these two reasons, one expects that the DNS flow field will quickly become decorrelated with respect to the experimental field in the near-wall region, which is well illustrated in Fig. 11. The decorrelated near-wall flow is transported away from the wall through turbulent diffusive flux. Effectively, a new and statistically independent flow develops at the wall and mixes with the bulk flow. A similar but much weaker effect is seen at the upper boundary of the turbulent inlet, where a blending scheme similar to that used at the wall was used to guarantee a smooth transition to the steady baseflow profile. Accordingly, the artificially 'blended' flow at y * >1.2 no longer shows correlation to the experimental flow, and is also (albeit more slowly) mixed toward the middle of the BL. As seen in Fig. 10, the large scale turbulent structures in the mid and outer layer are largely preserved throughout the transit length of Δx * =3 due to their higher momentum and the lower turbulent stress present to drive turbulent mixing, which like the near-wall flow would amplify the degree of decorrelation.
To very roughly estimate the stream-wise distance at which the influence of the decorrelated wall and upper boundary flows permeate the entire simulated TBL, a least squares trend line was calculated using segments from the upper and lower arbitrarily chosen iso-contour of the correlation coefficient r=0.75 . The iso-contour is plotted in light gray in Fig. 11, the least-squares linear trend line is plotted as a dotted black line and the segment of the r=0.75 contour used to generate the trend line is indicated in dark gray. When extrapolated, the linearized contour segments intersect by roughly x * ≈13.1 , indicating that by this transit distance, one would expect the simulated TBL flow to be largely decorrelated from the experimental flow.

Summary and conclusions
The present study has demonstrated a proof-of-concept approach for defining an unsteady inlet boundary condition for the DNS solver NS3D and evaluated the results in the canonical case of a ZPGTBL. In doing so, the potential for such a methodology as well as possible difficulties are explored. In the current case, a 4D-PTV experiment was used to obtain a 3D transient description of a ZPGTBL, therefore attention to method-specific metrics such as tracer nearest neighbor distance and the inherent limitations of spatial interpolation at finite tracer seeding densities are investigated. Alternate methods of obtaining an input field such as LES, VLES or URANS could be explored in future studies and would face different limitations.
The cutting-edge 4D-PTV technique is effectively unparalleled in its ability to deliver highly resolved, 3D transient datasets of turbulent flows. Even so, 4D-PTV has a limited maximum tracer density due to hardware constraints which in turn restricts the minimal resolvable turbulent structure size, after which the measurement is spectrally incomplete. Since information about the flow field is localized at tracer locations, the information density relative to local velocity gradient becomes critical near the wall for wall-bounded viscous flows. Interpolation to the DNS grid near the wall is therefore problematic because of the high ratio of tracer nearest neighbor radius to local turbulent structure size or velocity gradient.
Steps must be taken to 'reinstate' a no-slip condition at the wall, introducing artifacts such as an unrealistically steady or otherwise unphysical near-wall velocity profile which quickly rebalance after entering the DNS domain. This rebalancing process is realized as an initial relaminarizing tendency near the wall, causing a 'tripping effect' and enhanced turbulent production slightly off the wall (Fig. 8). As a result, a statistically independent flow develops at the wall (Fig. 11), which mixes upward into the rest of the boundary layer with increasing stream-wise position x.
As the DNS flow develops, smaller scale turbulent structures previously unresolved in the 4D-PTV measurement appear. This effect is well summarized in the wavenumber spectrum ( Fig. 9), where high wavenumber energetic content is recovered within the DNS domain. The wavenumber spectrum reaches a converged state after roughly 2.0−2.5⋅̄9 9 . This result is highly promising, as similar DNS studies of ZPGTBLs require nearly 20−30⋅̄9 9 of transit distance to reach such a state when initialized from a randomized field. The high-momentum, large scale turbulent structures are adequately resolved by the current method, as is evident by the consistency of the low-wavenumber regime in Fig. 9 and qualitatively through the comparison of large scale structures in Fig. 10. The physically well-resolved flow of the outer log layer and wake region therefore serves as a physically consistent boundary condition under which the decorrelated wall flow can quickly recover and converge to a fully resolved turbulent DNS flow. This behavior suggests the aptitude of the present method for the fully resolved simulation of flows developing under unsteady conditions, or for the implementation of more efficient DNS initialization schemes utilizing externally derived turbulent fields.