{'fraction_non_alphanumeric': 0.022986663246986408, 'fraction_numerical': 0.005321877404462683, 'mean_word_length': 4.98599117251967, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical issues in Lagrangian tracking and topological evolution of fluid particles in wall-bounded turbulent flows

The determination of the local Lagrangian evolution of the flow topology in wall-bounded turbulence, and of the Lagrangian evolution associated with entrainment across the turbulent / non-turbulent interface into a turbulent boundary layer, require accurate tracking of a fluid particle and its local velocity gradients. This paper addresses the implementation of fluid-particle tracking in both a turbulent boundary layer direct numerical simulation and in a fully developed channel flow simulation. Determination of the sub-grid particle velocity is performed using both cubic B-spline, four-point Hermite spline and higher-order Hermite spline interpolation. Both wall-bounded flows show similar oscillations in the Lagrangian tracers of both velocity and velocity gradients, corresponding to the movement of particles across the boundaries of computational cells. While these oscillation in the particle velocity are relatively small and have negligible effect on the particle trajectories for time-steps of the order of CFL = 0.1, they appear to be the cause of significant oscillations in the evolution of the invariants of the velocity gradient tensor.


Introduction
The study of the temporal evolution of wall-bounded flows is complicated by the range of motion experienced by the near-wall and outer flow, and the dependence of many fluid properties on the relative motion of the observer. In order to overcome this challenge many studies have focused on Galilean invariant quantities such as the velocity gradient tensor A ij (VGT) which governs the formation of intense vortical filaments and whose invariants (Q A , R A ) can be used to classify the local topology at any point in the flow, within the framework of critical-point theory [1], as indicated in Figure 1. Following this methodology a number of researchers have calculated spiralling conditional mean trajectories (CMT) in Q A -R A space similar to that shown in Figure 2, with associated time-scales, which represent the non-local mean topological evolution [2]. Typically these calculations involve the solution of the Lagrangian evolution equations for the Q A , R A invariants using the right-hand side of Navier-Stokes equations [3] in an Eulerian frame of reference, evaluated from direct numerical simulations (DNS) [2,4,5]. Alternatively, experimental measurements of the Eulerian evolution of the invariants can be used, and the convective terms estimated [6]. This enables an investigation of the mean evolution of the invariants, conditioned by their location in the Q A -R A space, for a given fluid domain, but does not allow the investigation of individual fluid particle traces, or account for the movement of particles throughout the domain. Given that the evolution of the invariants and their associated time-scales have been shown to be strongly dependent on the wall-normal position and, in the outer region of a turbulent boundary layer at Re τ ≈ 1000, can have evolution times of the order of 2.5 eddy-turnover times and a convected distance of approximately 20 boundary layer thickness (δ) [5], it remains to be seen how representative these CMT are of the true mean Lagrangian evolution of the flow, and the extent to which individual fluid particles follow this mean.  The Lagrangian evolution of individual fluid particles is also of interest when it comes to the study of the entrainment of irrotational fluid into a turbulent flow, such as the interface between turbulent and non-turbulent (T/NT) regions of a turbulent boundary layer, which have been shown to extend far into the core of the boundary layer [8]. The dynamics of this interface play an important role in the entrainment of mass, momentum and passive scalars into a turbulent flow and in the growth of the boundary layer. Advances in experimental measurement techniques and access to the high-resolution data of direct numerical simulations has lead to an increasing interest in the geometry, role of coherent structures and mechanisms associated with T/NT interfaces in both free and wall-bounded shear flows [9][10][11][12][13]. To date, most of the research regarding the T/NT interface has focused on the investigation and the classification of the interface and the surrounding flow, with less attention given to the Lagrangian view of the entrainment process. For instance, what happens to fluid particles in an irrotational region? how do they interact with the rotational flow and gain enstrophy? what is the topological evolution and time-scale associated with this process?
In order to answer these questions about the true topological evolution experienced by fluid particles in various regions of a turbulent wall-bounded flow, or the evolution and dynamics of a fluid particle as it is entrained across the T/NT interface, it is necessary to implement fluidparticle tracking and thereby determine the Lagrangian evolution of the VGT as it follows a particle at each time step of a DNS of the wall-bounded flow. Fluid-particle tracking has received significant attention in DNS of homogeneous isotropic turbulence [14,15] but, with the exception of [16], presumably owing to the higher computational demands, it has received relatively little attention in DNS of wall-bounded flows. In this paper we investigate the implementation of particle tracking in both a turbulent boundary layer and a fully developed channel flow, and examine the influence of different interpolation and particle time-marching approaches on the calculated particle paths, and the Lagrangian evolution of the VGT.

Particle Tracking in Direct Numerical Simulations of Wall-bounded flows
Codes that perform direct numeric simulations of fluids flow can be modified to enable the calculation of the Lagrangian evolution of the fluid by treating element or packets of the fluid as if they were particles that exactly follow the fluid dynamics, so that where x is the particle's position and u is the fluid flow velocity [14]. Given the simplicity of the above equation there should be no need to modify the time-marching of the principle code, assuming that it resolves all the relevant time-scales of the fluid flow. However, if during the time-marching the code uses an implicit or predictor step, a decision must be made as to whether or not to increment particle position at each intermediate step, assuming continuity is enforced at each step, and then re-calculate the particle velocity, or to only update the particle position at the end of the time-step. This later option only involves the use of the higher-order-accurate velocity fields, but will likely require a small overall time-step and more computation in order to achieve an accurate particle trajectory. Ultimately the accuracy of the fluid particle dynamics and the calculated particle trajectory will depend on both the magnitude of the time-step and the accuracy with which the velocity u( x) at the particle location can be determined from the computational cells of the fluid simulation. Assuming that the velocity components at the faces of each cell are calculated to a sufficient accuracy to capture the dynamics of the fluid, the accuracy of the particle velocity will still depend on the numerical interpolation technique used to compute the sub-cell velocity components. Owing to the expense of direct three-dimensional interpolation, such interpolation is typically determined via successive one-dimensional interpolations of the required grid points along each axis, which can also be advantageous when it comes to the interpolation of particle velocity from a computational domain that is split across a series of processors nodes. Optimal implementation of particle tracking depends on the DNS code in which the tracking is to be performed. In the present case, fluid-particle tracking is to be implemented both in the high-resolution incompressible boundary layer code, and in the fully developed channel code developed by the School of Aeronautics at the Universidad Politécnica de Madrid. Simulation of the turbulent boundary layer is performed using a fractional-step method with a three-step Runge-Kutta scheme for time-marching that is second-order accurate for velocity, and enables the use of variable pressure gradients and time-dependent boundary conditions [17]. Convective and viscous terms in the streamwise (x) and wall-normal (y) directions are calculated using staggered three-point compact finite differences. The computational domain is periodic in the span-wise (z) direction with velocity and pressure components expanded in a Fourier series so that the span-wise gradient can be calculated using spectral methods with the 2/3 rule to prevent aliasing. A pressure correction is applied using the Poisson equation at each sub-step of the Runge-Kutta. Velocity gradients and divergence for the Poisson correction are calculated by centred second-order finite differences in the x, y directions and spectral in the z direction. The present version of the code uses a combination of multi-processing (OMP) and message passing (MPI) for hybrid multi-thread and multi-core parallel processing [18] with the decomposition of the computational domain alternating between cross-stream y, z slabs and stream-wise pencils. Tracking is carried out in the slab decomposition, so that communication is only required in the streamwise direction. The channel-flow simulations are performed using the established pseudospectral method of Kim et al [19] with particle tracking incorporated into the parallel code of Hoyas & Jiménez [20]. Homogeneity in both the stream-wise and span-wise directions allows the pseudospectral solution of the Navier-Stokes equations rewritten in terms of wall-normal vorticity ω y and the Laplacian of v, hence eliminating the need for pressure. While computationally efficient in terms of the fluid phase, this formulation requires the generation of a realisable velocity field u i (x, y, z) in order to interpolate for the velocity at each particle. As in the case of the boundary layer, particle positions and velocities are updated at each substage of a semi-implicit thirdorder Runge-Kutta time marching scheme, hence requiring the generation of the fluid velocity field at each sub-step. In this formulation, derivatives in both x and z directions are handled spectrally, while wall-normal derivatives are computed via a fourth-order seven-point compact finite difference.

Sub-grid Interpolation methods
For structured grids, a variety of techniques can be used for sub-cell interpolation. If the parent simulation code involves periodic dimensions in an incompressible flow field, a spectral interpolation can be performed with no interpolation error for all scales resolved by the grid. For non-periodic dimensions compact (spectral-like) interpolation schemes such as those found in [21] and [22] may be used. It is important to note that, while providing excellent interpolation accuracy, both the spectral and compact methods involve significant computational cost that must be repeated for every particle in the domain at each time-step or Runge-Kutta sub-step. Previous applications of fluid particle tracking in DNS of turbulent flows have typically used more computationally efficient and scalable explicit Basis-or Hermite-spline based interpolation schemes [14,16] with Choi et al. [16] demonstrating that a four-point Hermite interpolation (degree = 8) provides results similar to the far more expensive spectral interpolation, albeit with derivatives at each grid point evaluated spectrally. Recent work of Stegeman et al. [23] indicates that even higher-order Hermite interpolation is required if the interpolation error is to be limited to only a few percent at the smallest wavenumber of the computation. Only the use of Basis and Hermite splines will be consider in this paper.
While both Basis-and Hermite-spline interpolation involves fitting and evaluating a local piecewise polynomial, Hermite splines require the value of the function (in this case velocity) and its derivatives at the two grid points nearest to the target location of the interpolation. In contrast Basis or B-splines represent a recursion of only the target quantity at multiple grid points.
Consider a scalar function (e.g. a single component of the fluid velocity) f (x) sampled with N points over a uniform grid where the parametric positionx between two consecutive nodes is defined from the spatial position x asx The polynomial coefficients a p are determined from the boundary conditions of the scalar and its derivatives atx = 0 and 1: This ensures that the derivatives of the fitted spline correspond to the derivative of f (x) and are continuous across grid points up to an order K = M/2. The boundary derivatives of the sampled function are determined by a central differencing method with an order K such that the spline coefficients a p are a linear function of the set of sampled points Substituting the above expression for the spline coefficients into Equation (2) the spline can be presented as so that each spline of degree M has an array of (M − 1) Interpolation of a three-dimensional function can then be performed via successive one-dimensional interpolations of the required grid points along each axis. Given that the spline represents a piecewise polynomial, it can also be used to determine the gradient of a scalar at a given interpolated location, albeit at lower order of accuracy. Both this method and the interpolation of the velocity gradient from scalar fields of each component of ∂u i /∂x j are examined in section 5. From the above definition, a cubic Hermite spline consists of order M=4, and requires the calculation to second-order accuracy of derivatives up to first order at the nearest grid points. This can be obtained from a second-order central difference. Following the formulation of Equation (6) a cubic B-spline will have pre-determined coefficients An alternative form to this Hermite interpolation was considered in the turbulent channel flow, where rather than using higher-order derivative boundary conditions at the two nearest grid points, a higher-order interpolation was instead performed using the zero-and first-order derivatives at the four nearest grid points [24]. In one dimension, a four-point Hermite spline can be represented as where H and G are basis functions and f and all its first-order derivatives are evaluated at the grid points. The relative error of both Basis-and Hermite-spline interpolation is due to the truncation error of the polynomial approximation as well as the error in the finite-difference approximation of the boundary derivatives in the case of the Hermite spline. This error can be explored using a plane-wave error analysis which provides the transfer function of the spline interpolation operator. A continuous plane-wave with a biased phase can be expressed in terms of the parametric coordinate system as   Substituting Equation (8) into Equation (6), the interpolated value can be expressed as The difference between the true value of the plane wave and its interpolated value can be determined from a Monte Carlo simulation in the parameter space consisting of the parametric wavenumber k∆, the phase φ, and the location of the interpolated pointx relative to the nearest grid point. Figure 3 shows the magnitude of the interpolation error, defined as the magnitude of the difference between the interpolated value f I (x) and the true value of the scalar f (x) at the same location with respect to the grid points. Not surprisingly the error increases as the distance between the interpolation point and the grid points increases, for both Hermite and B-splines, particular at higher wavenumbers. This change is far more gradual in the case of the B-spline but, as shown in Figure 4, the total error in the Hermite spline is lower than that of the B-spine for the same wavenumber.  From Equation 9, a complex transfer function for this interpolation can be defined as which is only a function of the interpolation positionx and the relative wavenumber k∆. The modulus and argument of the transfer function defines the relative difference in magnitude and absolute difference in phase between the interpolated plane-wave and the original planewave, respectively. The contour plots in Figure 5 show the modulus and argument for the Hermite spline. The transfer function of the B-spline is similar, but shows less variation in interpolation position, as observed by the relative magnitude of the interpolation error. In both cases the argument of this transfer function suggests that both interpolation schemes may introduce oscillations in the higher wavenumbers. The absolute magnitude of these interpolation errors and the effect on the particle trajectories and the evolution the of VGT and its invariants will depend on the energy content of the scalar (the fluid-velocity field) at different wavenumbers, the grid spacing, and the chosen degree of the fitted spline. Figure 6 shows the effect of varying the degree of Hermite interpolation. Assuming negligible error in the scalar function at the grid points, a higher-order spline reduces the error at all wavenumbers. While the reduction in error from a linear (degree = 2) to a cubic (degree = 4) spline is significant, the incremental reduction in error decreases with higher order. If the effective cut-off wavenumber of the simulation is taken at kδ = 2/3, following the 2/3 rule for de-aliasing, Figure 6 indicates that the error at this wavenumber should be approximately 22% for a cubic Hermite spline, decreasing to 5% for a Hermite spline of degree 20.
Higher-order interpolation requires the use of higher-order gradients, either directly in the case of Hermite interpolation, or through a larger coefficient matrix in the case of B-splines. Both cases ultimately require a wider interpolation kernel of size N P = M , equal to the degree of the spline. This significantly increases the computation required for each particle. Given that the entire fluid phase must be computed for each time-step and that Lagrangian statistics require the calculation of numerous particles, it is therefore desirable to run the DNS with as many particles as possible. Naturally, this requires that the additional computation imposed by the particles be kept to a minimum. The benefits of higher-order interpolation are assessed in the context of the particle tracers in section 5.

Fluid-Particle Tracking Implementation
In order to test the benefits and possible need for high-order interpolation methods and to enable the processing of high-resolution boundary-layer simulations, the DNS code was modified to enable the use of Hermite spline interpolation of any specified order for both uniform and irregular grid spacing. Following the formulation for Hermite splines discussed in section 3, derivatives up to the required order are calculated at grid points nearest to the particle location using central or forward-or backward-biased finite-difference schemes, depending on the location of the particle with respect to the domain boundaries. Coefficients for these finite-difference schemes were computed for different-order schemes during runtime, following [25]. The velocity of each particle was determined at their initial location before the first Runge-Kutta time-step by converting one velocity component into physical space, and then performing three-dimensional interpolation of that velocity component at each particle location, before considering the next component. Interpolation was performed by first fitting one-dimensional Hermite splines of order M in the y, z cross-stream directions over a three-dimensional interpolation kernel of M ×M ×M points around each particle, with final interpolation performed in the x direction. Structuring the interpolation in this format has the advantage that when a particle's interpolation kernel extends across the boundaries between computational nodes, the cross-stream interpolation can be handled by separate nodes. This reduces the communication between nodes to at most M/2 quantities for each velocity component.
Particle positions where updated at the start of each Runge-Kutta sub-step using their velocity from the end of the previous step. Velocities were then re-interpolated at this updated particle location after the fluid velocity was updated by the right-hand side of the governing equation, and after the Poisson correction for the divergence. The positions of all master (centre of particle is contained in the node) and slave particles (part of the interpolation kernel of a particle from an adjacent node) are stored for each computational domain, as done by Uhlmann [26]. With each update of the particle position, a test is performed to see whether a particle or its interpolation kernel crosses the boundary of its parent node, forcing the transition of a slave to a master or the formation of a new slave, respectively.
The velocity gradient tensor (VGT) at each particle was evaluated by first determining the velocity gradients at the centre of each cell using second-order central finite differences, consistent with the methodology used for the Poisson correction. It was assumed that the error associated with this scheme would be insignificant in comparison to the error introduction by the interpolation (at least for lower-order schemes), although this assumption has yet to be tested. Each component of the VGT was then interpolated following the same procedure used for the  velocity field. Tests were also performed to see if the velocity gradients could be determined from differentiation of the spline fits as suggested by Yeung and Pope [14], which does not require the differentiation of the velocity field or additional interpolation, and is therefore far more efficient.
A similar process was performed in the channel flow using both cubic (degree = 4) B-splines and the four-point Hermite interpolation discussed in section 3. Particle positions and velocities were similarly computed at each Runge-Kutta sub-step. Pseudo-spectral methods were used to compute the velocity gradient at each grid point, from which each component of the velocity gradients tensor was calculated at each particle location via the same spline interpolation.

Results
To test the effect of the Hermite spline order and the size of the time-step on the particle trajectories and on the Lagrangian evolution of the invariants of the VGT, a turbulent boundary layer DNS was run over two moderate resolution grids whose dimensions and grid spacing are listed in Table 1. Initial testing was performed on the smaller grid in order to test the basic implementation of the code and the interpolation before extending the code to more realistic simulations.
In the absence of known particle trajectories, the accuracy of a given interpolation scheme was tested based on how much the particle position differed from that produced by a higherorder interpolation or a shorter time-step ∆t, given that position error is a function of both the interpolation error and the time-step. The time-step from this point on will be quoted in terms of Courant-Friedrichs-Lewy number CF L = U ∆t/∆x, where U is the free-stream velocity and ∆x is the length of the computational cell. In addition to this, for an incompressible flow, the divergence ∂u i /∂x i measured at the particle location can be used as an indication of the accuracy of the velocity gradients calculated for each particle. Figure 7 shows the difference in the divergence obtained when calculating the velocity gradients at each particle from the gradient of a spline fitted to the velocity field (as suggested by [14]) compared to fitting a separate spline to the gradients of the velocity field at each grid point. Separate interpolation of each component of the VGT provides a divergence that is typically five orders of magnitude smaller than that obtained by differentiating the spline. Interpolating each component of the VGT is considerable more expensive, but provides a far more accurate estimation of the sub-grid particle velocity gradients. The behaviour for the different interpolation and time-steps can be illustrated by considering the behaviour of an individual particle trajectory. Similar patterns were observed for most particles although, naturally, the degree to which a given particle trajectory diverges for different interpolation schemes will ultimately depend on the dynamics of the flow field at that time and position. Figure 8 shows the trajectory followed by a particle started at a wall-normal position of y + = 126.7 using a series of Hermite interpolations. Particle positions and velocities are updated for each Runge-Kutta time-step with the exception of the final trajectory, where the particle positions are instead updated using a first-order Euler time marching with CFL = 0.01. This final trajectory was generated in order to ensure that the use of the velocity fields inside the Runge-Kutta time-step did not have an adverse effect on the calculation of the particle's velocity. Comparison of the trajectories indicates that, for the same interpolation order, updating the particle position within the Runge-Kutta substeps enables time-steps at least an order of magnitude longer than for the Euler method, with negligible influence on the particle trajectory. Figure 9 shows the variation in particle velocity as a function of time for each of the interpolation methods discussed above, demonstrating that both velocity and the consequent position are highly sensitive to the order of interpolation scheme and to the size of the timestep used in the simulation. The use of linear (degree = 2) interpolation results in significant oscillations in each component of the particle velocity, which is ultimately responsible for the difference in the particle trajectory. The use of cubic (degree = 4) interpolation reduces the amplitude of these oscillations, but the benefit of moving to a degree = 8 is less clear. In all cases these oscillations maintain approximately the same period and are strongly correlated with the particle position relative to the Eulerian simulation grid points.
The Lagrangian evolution of the Q A , R A invariants of the VGT following the same particle are shown in Figure 10 where the trajectories originate in the top left corner of the domain, which corresponds to the particle being in a region of stable focus stretching (see Figure 1). The invariants shown have been normalised by the mean of the second invariant of the rate of rotation tensor Q W , which is directly proportional to the enstrophy [5]. While it is not clear what the  Figure 10. Lagrangian evolution of Q A , R A invariants of the velocity gradient tensor in the boundary layer simulation calculated for a single particle using Hermite splines of varying order over a series of different time-steps.
Re τ (L x , L y , L z )/δ ∆x + , ∆y + , ∆z + N x , N y , N z P articles 2000 2π × 2 × π 10 × 9 × 5 1536 × 633 × 1536 7.296 × 10 5 Lagrangian trajectory of the invariants of an individual particle should look like, results shows a similar pattern of oscillation to the particle velocity. This oscillation is again reduced by using a higher-order interpolation although, unlike in the case of the particle velocity, these oscillations continue to reduce up to an interpolation degree of at least 8. This decreased oscillation for higher-degree splines suggests that the oscillations are at least partly due to the interpolation method.
The channel flow has been tested most extensively in the configuration given in Table 2. Figure 11 illustrates oscillations in the time series of quantities interpolated using both cubic B-splines and four-point Hermite splines. As in the case of the boundary layer, these oscillations do not correspond to non-monotonic behaviour of the unsteady flow in the surrounding grid points. As shown by Figure 12, oscillations in the velocity gradients are more noticeable. It should be noted that the magnitude of the oscillations is typically small but, as shown in Figure  10, they can have a significant effect on the observed Lagrangian evolution of the invariants, and will have consequences when it comes to the calculation of Lagrangian accelerations. Despite significant differences between the resolution and the computation of velocity gradients in both the boundary layer and channel flow simulation, the similarity in the computed Lagrangian velocity and velocity gradient tracers suggests that the accuracy of the Lagrangian statistics is highly sensitive to the interpolation scheme. Given the shape of the transfer function of the interpolation method (see Figure 5) and the spacing of these oscillation with respect to the grid points, it is likely that these oscillations are purely the result of the sub-grid velocity interpolations. It may therefore be possible to low-pass filter these oscillations without effecting the Lagrangian statistics, assuming that the time-step of the simulation is small enough that these oscillations do not significantly influence the calculated particle trajectory, and given that there should be no appreciable energy in the fluid phase for parametric wavenumbers k∆ > 2π/3 under the 2/3 de-aliasing of the DNS. Technically, this filtering should be a function of positional variation in all three-directions and not a function of time. Given that the size of the computational cell varies with the wall-normal position, and that oscillations occur as the particle moves with respect to each face of the computational cell, the implementation of this filtering is non-trivial.

Conclusions
In order to enable the calculation of Lagrangian statistics and time-scale in wall-bounded flows, the Lagrangian tracking of fluid particles has been implemented in direct numerical simulation codes for both a boundary layer and for a fully developed channel flow. Sub-grid interpolation was implemented using both cubic B-spline, four-point Hermite and higher-order Hermite interpolation schemes. Despite significant differences in cell size, the computation of velocity gradients, and the implementation of the interpolation, the particle velocities and gradients in both the boundary layer and the channel flow showed similar oscillations in their Lagrangian tracers. These oscillations are highly correlated with the movement of particles with respect to the computation grid and, given the similar pattern in the phase error of the interpolation method, appear to be due to the sub-grid interpolation method. While these oscillations in the particle velocity are relatively small and have negligible effect on the particle trajectories for time-steps on the order of CFL = 0.1, they appear to be the cause of more significant oscillations in the evolution of the invariants of the velocity gradient tensor. The extent to which these interpolation errors may influence the Lagrangian statistics is not yet clear.