Direct Numerical Simulation of Turbulent Channel Flow on High-Performance GPU Computing System

The flow of a viscous fluid in a plane channel is simulated numerically following the DNS approach, and using a computational code for the numerical integration of the Navier-Stokes equations implemented on a hybrid CPU/GPU computing architecture (for the meaning of symbols and acronyms used, one can refer to the Nomenclature). Three turbulent-flow databases, each representing the turbulent statistically-steady state of the flow at three different values of the Reynolds number, are built up, and a number of statistical moments of the fluctuating velocity field are computed. For turbulent-flow-structure investigation, the vortex-detection technique of the imaginary part of the complex eigenvalue pair in the velocity-gradient tensor is applied to the fluctuating-velocity fields. As a result, and among other types, hairpin vortical structures are unveiled. The processes of evolution that characterize the hairpin vortices in the near-wall region of the turbulent channel are investigated, in particular at one of the three Reynolds numbers tested, with specific attention given to the relationship that exists between the dynamics of the vortical structures and the occurrence of ejection and sweep quadrant events. Interestingly, it is found that the latter events play a preminent role in the way in which the morphological evolution of a hairpin vortex develops over time, as related in particular to the establishment of symmetric and persistent hairpins. The present results have been obtained from a database that incorporates genuine DNS solutions of the Navier-Stokes equations, without superposition of any synthetic structures in the form of initial and/or boundary conditions for the simulations.


Introduction
The flow of a viscous fluid in a channel has been investigated numerically by several authors in the recent past, becoming a reference case for the study of wall turbulence with DNS.
In the aforementioned channel-flow simulations, interesting results have been obtained, as related in particular to the evolution of the fluctuating-velocity statistical moments with the Reynolds number (see also Alfonsi [24], Marusic et al. [25], Smits et al. [26], Kim [27], Jiménez [28]).
Though, there are several aspects of the channel-flow case that can be further investigated, in particular related to the processes of development of the wall-turbulence flow structures (see also Alfonsi and Primavera [29]).
In the present work, three DNS channel-flow database have been calculated, respectively, at friction-velocity Reynolds numbers Re τ " 200, 400 and 600 (one can see at the Nomenclature for the meaning of symbols and acronyms used).Statistical moments of the fluctuating-velocity field have been computed and compared with results obtained by other authors, obtaining a rather good agreement with the latter.The vortex-detection technique of the imaginary part of the complex eigenvalue pair of the velocity-gradient tensor (the λ ci or swirling-strength criterion) as introduced by Zhou et al. [30] has been applied to the computed fluctuating-velocity fields.As a result, and among other shapes, hairpin-like vortical structures are unveiled.The processes of evolution that characterize the hairpin vortices are investigated giving particular attention to the relationship that exists between vortex dynamics and the occurrence of Q2 and/or Q4 quadrant events.Interestingly, it is found that the physical condition for the development of a complete and stable hairpin is twofold, namely: (i) the development of ejections distributed on spheric-like isosurfaces behind an initial Ω-shaped vortex filament; (ii) the subsequent development of sweeps distributed on elongated isosurfaces adjacent to the external sides of hairpins' heads and necks.
The present work is organized as follows.Section 2 contains an outline of the numerical technique used for the solution of the Navier-Stokes equations in the plane-channel-flow computing domain, in Section 3 a concise presentation is given of the vortex-detection method used for flow-structure extraction, and in Section 4, the numerical simulations are described.In Section 5, the results are presented and compared with data obtained by other authors, mainly in terms of turbulence statistics, while in Section 6, the results of the simulations are presented in terms of vortical structures and quadrant events.Concluding remarks are given at the end.

Numerical Techniques
The three-dimensional time-dependent Navier-Stokes equations for incompressible fluids are considered in non-dimensional, conservative form (Einstein summation convention applies to repeated indices, i, j = 1, 2, 3): Variables and operators are nondimensionalized by the channel half-height (h) for lengths, the wall-shear velocity (u τ ) for velocities, the group (ρu 2 τ ) for pressure, and (h{u τ ) for time, being Re τ " u τ h{ν the friction-velocity Reynolds number (ρ is fluid density, ν is fluid kinematic viscosity.Note that, for simplicity, the symbols of both dependent and independent variables have not been altered in switching from the dimensional to the dimensionless formalism).The computing domain (Figure 1) is considered homogeneous along the x (streamwise) and z (spanwise) directions, so that Equations ( 1) and ( 2) are Fourier-transformed accordingly: where the superscript (ˆ) indicates variables in Fourier space, and k 2 " k 2 x `k2 z .The nonlinear terms in the momentum Equations (3a-c) are evaluated pseudospectrally by anti-transforming the velocities in physical space to perform the products (FFTs are used).Here, in order to avoid errors in transforming the results back to Fourier space, the discrete Fourier transforms are applied on "3n/2" points along each homogeneous direction.

  
where the superscript (^) indicates variables in Fourier space, and . The nonlinear terms in the momentum Equations (3a-c) are evaluated pseudospectrally by anti-transforming the velocities in physical space to perform the products (FFTs are used).Here, in order to avoid errors in transforming the results back to Fourier space, the discrete Fourier transforms are applied on "3n/2" points along each homogeneous direction.Due to the presence of the steepest variable-gradients near the walls, and in order to obtain a suitable spatial resolution in the calculations, a grid-stretching law of hyperbolic tangent type is used for the grid points along y (the direction orthogonal to the walls): where y indicates the uniform grid and PP, QQ are two parameters characterising the distribution.The partial derivatives along y are calculated according to grid-point distribution Equation (5) using second-order centered finite-difference expressions.For time advancement, the third-order Runge-Kutta procedure originally introduced by Le and Moin [31] is implemented.For each Fourier mode one has: where l = 1,2,3 denote the Runge-Kutta sub-steps, and where: Due to the presence of the steepest variable-gradients near the walls, and in order to obtain a suitable spatial resolution in the calculations, a grid-stretching law of hyperbolic tangent type is used for the grid points along y (the direction orthogonal to the walls): y str " PPy `p1 ´PPq ˆ1 ´tanh rQQ p1 ´yqs tanh rQQs where y indicates the uniform grid and PP, QQ are two parameters characterising the distribution.The partial derivatives along y are calculated according to grid-point distribution Equation (5) using second-order centered finite-difference expressions.For time advancement, the third-order Runge-Kutta procedure originally introduced by Le and Moin [31] is implemented.For each Fourier mode one has: where l = 1,2,3 denote the Runge-Kutta sub-steps, and where: are the advective-and diffusive terms, respectively.Both terms are treated explicitly and, in Equation ( 6), α l , β l , γ l , ζ l assume constant values, so that the time advancement results are third-order accurate in the convective part, and second-order accurate in the diffusive: ; (9a) Time marching is coupled with the fractional-step method.In Equation ( 6), the velocity and pressure fields are decoupled, so that two distinct expressions are generated.At each Runge-Kutta sub-step (l) and for each Fourier mode (ˆ), an intermediate velocity field is introduced (superscript *): where, by applying the divergence operator to Equation (11) (and so enforcing mass conservation) one obtains a Poisson equation for the pressure, to be solved at each sub-step (l): The final values of the velocity are obtained from Equation (11).No-slip boundary conditions at the walls, and cyclic conditions in the streamwise and spanwise directions, have been applied to the velocity (for further details on the numerical algorithm one can refer to Passoni et al. [32][33][34]).

Flow-Structure Extraction
Among the different techniques used for the extraction of the coherent structures of turbulence (see Wallace [35], Alfonsi [36], Alfonsi and Primavera [37,38], among others), the swirling-strength criterion as devised by Zhou et al. [30] has been used.The latter is concisely summarized here.By considering the system of the governing equations, an arbitrary point O can be chosen in the field, and a Taylor-series expansion of each velocity component can be performed in terms of space coordinates with the origin at O, so that the first-order pointwise linear approximation at that point becomes: (A ij " Bu i {Bx j is the velocity-gradient tensor).If O is located at a critical point, the zero-order terms in Equation ( 13) are zero.From the characteristic equation of A ij one has: where: are the scalar invariants of the velocity-gradient tensor (tr is trace, det is determinant).In the case of incompressible flow, P " 0, and: where the discriminant of the characteristic equation of A ij becomes: When Dsc ą 0, the velocity-gradient tensor has one real eigenvalue (λ 1 " λ r ), and a pair of complex-conjugate eigenvalues (λ 2 , λ 3 " λ cr ˘iλ ci ).Zhou et al. [30] adopted the criterion of identifying vortices by visualizing isosurfaces of prescribed values of the imaginary part of the complex-eigenvalue pair of the velocity-gradient tensor.The swirling strength (λ ci ) represents a measure of the local swirling rate inside a vortical structure, so that isosurfaces of the imaginary part of the complex eigenvalue pair of the velocity-gradient tensor can be used to visualize vortices (the strength of stretching or compression is given by λ r ).The method is frame independent.It automatically eliminates regions having no local spiralling motion (due to the fact that the eigenvalues are complex only in regions of local circular or spiralling streamlines), and has proven to give rather satisfactory results in several different cases (see [39], among others).As concerns the choice of the threshold value of the swirling strength chosen for structure representation ( λ ci | th ), one can refer to Alfonsi and Primavera [40].

Numerical Simulations
Direct numerical simulations have been executed in the plane-channel domain (Figure 1) with dimensions, grid points, resolutions and Reynolds numbers as reported in Table 2.As concerns the calculation of the Kolmogorov microscales, they have been evaluated by estimating the average rate of dissipation of turbulent kinetic energy per unit mass (ε).
This method was first introduced by Bakewell and Lumley [41], where in the case of the plane channel, one has: The calculations have been executed on a specially-assembled hybrid multicore/manycore computing architecture.The system includes: (i) 2 Intel Xeon 5660 exa-core CPU processors (12 cores) at 2.8 GHz, with 48 GB GDDR3 RAM; (ii) 3 Nvidia C-1060 (Tesla) 240-core GPU boards (720 computing cores) at 1.3 GHz, each with 4 GB GDDR3 RAM at 102 GB/s (12 GB available); (iii) 1 Nvidia GTS-450 (GeForce) 192-core GPU board at 1804 MHz, with 1 GB GDDR5 RAM at 57.7 GB/s (mainly used for visualization); (iv) storage system including 5 Hard Drives at 7200 rpm, for a total supply of 5 TB.
Each GPU Nvidia C-1060 board handles a multiprocessor unit, the latter organized in 30 processors.Each processor includes eight floating-point units, 16 kB shared memory, and 4 GB of GDDR3 memory at 102 GB/s bandwidth.The process of implementation of the numerical algorithm (Section 2) on the above computing architecture is described in detail in Alfonsi et al. [42].The possibility of running the computational code on different partitions of the hybrid computing machine has enabled the execution of the numerical simulations in remarkably limited runtimes (Table 3).According to a procedural viewpoint, the initial transient of the flow in the channel has been first simulated, the turbulent statistically-steady state has been reached, and thereafter (for each value of the Reynolds numbers tested) 500,000 time steps of the statistical steady state have been gathered, with temporal resolution ∆t " 1 ˆ10 ´4h{u τ .The flow fields have been saved every given ∆ t interval, finally giving a 500¨∆t `database for each of the Reynolds numbers tested.The adequacy of length and span of the computing domain has been tested by verifying that the velocity fluctuations at streamwise and spanwise separations on half the domain dimensions were uncorrelated.The adequacy of the grid resolution has been also tested, through the analysis of the one-dimensional energy spectra.It has been verified that the energy densities associated to the high wavenumbers are up to nine orders of magnitude lower than those corresponding to the low wavenumbers (for more details on these issues one can refer to Ciliberti [43]).

Turbulence Statistics
In Table 4 and Figures  one-dimensional energy spectra.It has been verified that the energy densities associated to the high wavenumbers are up to nine orders of magnitude lower than those corresponding to the low wavenumbers (for more details on these issues one can refer to Ciliberti [43]).

Turbulence Statistics
In Table 4 and Figures 2-6, results are presented in terms of turbulence statistics.Figure 2 reports the values of the Reynolds shear stress ( ) in wall coordinates, in a comparison with the data.
).In Figure 6, the values of the production terms ( K P ), the transport terms ( K T ), the diffusion terms ( K D ), and the          Figure 3 reports the rms values of the velocity fluctuations (u1 rms , v 1 rms , w 1 rms ) in a comparison with those of Moser et al. [6].In a similar manner, Figure 4 reports the skewness factors of the velocity fluctuations (S u 1 , S v 1 , S w 1 ), and Figure 5 the flatness factors (F u 1 , F v 1 , F w 1 ).In Figure 6, the values of the production terms (P K ), the transport terms (T K ), the diffusion terms (D K ), and the dissipation terms (ε K ) of the turbulent-kinetic-energy transport equation (K " u 1 i u 1 i {2), as calculated in the present work, are reported, again in a comparison with Moser et al. [6].
As concerns mean-flow quantities, in Table 4 the values of a number of quantities [u τ ,u b ,u c , (u b {u τ ), (u c {u τ ), (u c {u b ), C f b ] as calculated in the present work are reported in a comparison with the experimental data of Dean [44].As for mean-velocity distribution, the linear distribution (u `" y `) is satisfactorily followed in the viscous sublayer (y `ă 5), while at larger distances from the wall (y `ą 30), the logarithmic distribution (u `" 1 κ lny ``C) is also satisfactorily followed, with κ " 0.4 and C " 5.5 (not shown).
As concerns fluctuating velocity components (Table 4), the peak values of the Reynolds shear stress (´u 1 v 1 ˇˇp eak ), those of the rms-fluctuating streamwise velocities ( u 1 rms ˇˇpeak ), and the fluctuating streamwise velocities skewness factors ( S u 1 | peak ) [and their positions ( y `ˇ´u peak , y `ˇSu 1 peak )], exhibit a good agreement with the values of Moser et al. [6].Moreover, it has been found that the peak values of the rms-fluctuating streamwise velocities and their positions satisfactorily follow the expressions devised by Mochizuki and Nieuwstadt [45] as a function of Re τ , as deduced from the analysis of a large number of experimental works (see also Alfonsi [46]): y `ˇu 1 rms peak " 0.00020Re τ `14.6 The values of the skewness factors of the streamwise fluctuations (S u 1 , Figure 4) are rather close to zero at the position of rms-fluctuating-streamwise-velocities peak values ( y `ˇu of the wall-normal fluctuations (F v 1 , Figure 5) also significantly differ from the Gaussian ones, and assume remarkably-high values approaching the walls, so unveiling the highly-intermittent character of the normal-velocity fluctuations near the walls.In particular (Table 4), approaching the walls, the flatness factor F v 1 assumes values that result in an excellent agreement with the data of Moser et al. [6] (see also at Xu et al. [47] and Alfonsi [48]).By looking at the budgets of the mean turbulent kinetic-energy (K, Figure 6), it can be noted that, at y `ą 30, the homogeneous character of the flow is reasonably confirmed, while, moving toward the wall, the turbulent-transport rate becomes relevant.The turbulent-transport term is a consuming term approaching the wall, and a producing term near the wall.Close to the wall, the dissipation rate balances the viscous-diffusion and pressure-diffusion rates, and, at the wall, the dissipation rate is nonzero, while being almost equal to the viscous-diffusion rate.
Overall, the comparison of the results of present work with data obtained by other authors, both numerically and experimentally, is rather satisfactory.

Flow Structures
After the application of the λ ci vortex-detection technique (Section 3), the flow field in the channel appears to be highly populated by turbulent structures adjacent to both the upper and the lower wall of the computing domain, with a wide range of inclination angles.The majority of them has no definite shape.Side views of the phenomenon (at Re τ " 200 and Re τ " 600) are given in Figure 7 at a generic instant.It can be noticed how the turbulent structures are noticeably smaller in size and more numerous in the case of Re τ " 600 with respect to Re τ " 200, as expected.In particular, Figure 7 shows that several structures also exist outside the buffer layer, protruding toward the center of the channel.Figure 8 shows a general view of vortical structures at Re τ " 600 on both walls of the computing domain.Here, the external surfaces of the structures are colored with the values of the local streamwise velocity (reddish indicates high values, greenish indicates low values).It can be noted, as expected, how the streamwise velocities increase towards the center of the computing domain, with respect to the more peripheral zones.Figure 9 shows a general view of vortical structures in conjunction with isosurfaces of Q2 and Q4 quadrant events (ejection and sweeps) at Re τ " 400, on the lower wall of the computing domain, at a generic instant.It can be noticed how flow structures and quadrant events are more densely located in streamwise-elongated zones of the domain, the latter being separated from the adjacent by low-speed streaks (arrows, see also Kline et al. [49]).
Through Figures 10-13 the process of evolution in time is shown with a single hairpin-like vortical structure.At t `" 450 (Figure 10), the onset is represented by the process of formation of a single, isolated, two-leg, symmetric and stable hairpin.It can be noted how an ejecting surface is pushing the perspective hairpin upwards (actually a Ω-shaped vortex filament) while, mainly on its left side (the flow goes from left to right), a sweeping surface starts to develop.Moreover, Figure 10b shows that the head of the structure is subjected to stretching, while the neck and legs are subjected to compression.

Re
, on the lower wall of the computing domain, at a generic instant.It can be noticed how flow structures and quadrant events are more densely located in streamwise-elongated zones of the domain, the latter being separated from the adjacent by low-speed streaks (arrows, see also Kline et al. [49]).
Through Figures 10-13, the process of evolution in time is shown with a single hairpin-like vortical structure.At 450   t (Figure 10), the onset is represented by the process of formation of a single, isolated, two-leg, symmetric and stable hairpin.It can be noted how an ejecting surface is pushing the perspective hairpin upwards (actually a  -shaped vortex filament) while, mainly on its left side (the flow goes from left to right), a sweeping surface starts to develop.Moreover, Figure 10b shows that the head of the structure is subjected to stretching, while the neck and legs are subjected to compression.11 and 12) during which the ejecting surface keeps pushing the head of the structure upwards, and the sweeping surface, adjacent to the right side of the structure, further grows, so that both the right and the left portions of the structure neck become adjacent, externally to the sweeping isosurface, and internally to the ejecting isosurface.The head of the structure (Figures 11b and 12b) continues to be stretched under the action of the ejecting surface, while neck and legs are subjected to compression, due to the action of the sweeping surfaces.

At instant
453   t (Figure 13) the ejecting surface starts to extinguish, while the sweeping surface now exerts its characteristic stabilizing action onto the entire external structure of the now fully-formed hairpin vortex (head, neck and legs).From Figure 13b it can be also noted that the head of the vortical structure is subjected to a less intense stretching, due to the gradual process of extinction of the previously upward-pushing underlying ejecting isosurface, while the legs of the structure are subjected to compression when subjected to the action of the overlying sweeping surfaces.
Figures 14-17 show the evolution in time of a more complex aggregate of vortical structures in a different portion of the computing domain.11 and 12) during which the ejecting surface keeps pushing the head of the structure upwards, and the sweeping surface, adjacent to the right side of the structure, further grows, so that both the right and the left portions of the structure neck become adjacent, externally to the sweeping isosurface, and internally to the ejecting isosurface.The head of the structure (Figures 11b and 12b) continues to be stretched under the action of the ejecting surface, while neck and legs are subjected to compression, due to the action of the sweeping surfaces.

At instant
453   t (Figure 13) the ejecting surface starts to extinguish, while the sweeping surface now exerts its characteristic stabilizing action onto the entire external structure of the now fully-formed hairpin vortex (head, neck and legs).From Figure 13b it can be also noted that the head of the vortical structure is subjected to a less intense stretching, due to the gradual process of extinction of the previously upward-pushing underlying ejecting isosurface, while the legs of the structure are subjected to compression when subjected to the action of the overlying sweeping surfaces.
Figures 14-17 show the evolution in time of a more complex aggregate of vortical structures in a different portion of the computing domain.The process continues through instants t `" 451 and t `" 452 (Figures 11 and 12) during which the ejecting surface keeps pushing the head of the structure upwards, and the sweeping surface, adjacent to the right side of the structure, further grows, so that both the right and the left portions of the structure neck become adjacent, externally to the sweeping isosurface, and internally to the ejecting isosurface.The head of the structure (Figures 11b and 12b) continues to be stretched under the action of the ejecting surface, while neck and legs are subjected to compression, due to the action of the sweeping surfaces.
At instant t `" 453 (Figure 13) the ejecting surface starts to extinguish, while the sweeping surface now exerts its characteristic stabilizing action onto the entire external structure of the now fully-formed hairpin vortex (head, neck and legs).From Figure 13b it can be also noted that the head of the vortical structure is subjected to a less intense stretching, due to the gradual process of extinction of the previously upward-pushing underlying ejecting isosurface, while the legs of the structure are subjected to compression when subjected to the action of the overlying sweeping surfaces.
Figures 14-17 show the evolution in time of a more complex aggregate of vortical structures in a different portion of the computing domain.At 220   t (Figure 14) two main primary  -shaped vortex filaments are visible at the center of the field (arrows).Right below the head of each filament, the internal space of the structure is occupied by an ejecting isosurface, again showing that, in this initial phase, ejections represent the main mechanism according to which the heads of the structures are raised upwards.Correspondingly (Figure 14b), the process of progressive stretching of the vortex head (actually still a double head) also starts.
In Figure 15, the flow field at 221   t is shown.The heads of both structures 1 and 2 continue to raise (being subjected to stretching, Figure 15b) due to the upward-pushing action of the ejections.The sweeping surface mainly maintains its position adjacent to the right side of hairpin 1, causing the compression of the contiguous vortical-structure neck (Figure 15b).Note again that the ejecting isosurface below hairpin 1 is almost equally developed in the streamwise, spanwise and normal-to-the-wall directions (a spheric-like surface), so guaranteeing that the pushing action of the ejecting surface against the internal part of the vortical structure is exerted almost uniformly.It can be also noticed that the presence of sweeping isosurfaces adjacent to hairpin 2 is not so evident as in the case of hairpin 1, and this is the reason why vortex 2 will be destroyed much sooner with respect to vortex 1 (not shown).At 222   t (Figure 16), the vortical structures continue their development.The heads of hairpins 1 and 2 continue to raise under the residual influence of the ejection events (the underlying ejecting isosurfaces are now reduced) while a well-defined sweeping isosurface assumes its definite position, adjacent to the external side of head, neck and legs of hairpin 1.The action of compression (Figure 16b) exerted by the sweeping vortex-overlying surface, in particular onto the neck of hairpin 1, becomes more intense.At 223   t (Figure 17), the process of evolution of the vortical structures continues in a similar manner, with respect to the previous instants.The action of the ejection events is decreasing, while the sweeping isosurface acts towards the maintenance of the stability of hairpin 1.Also, in this case, the result is hairpin 1 becoming a two-leg, symmetric, and stable vortical structure.

Concluding Remarks
The Direct Numerical Simulation (DNS) of the turbulent flow of incompressible fluid in a plane channel has been executed at three values of the friction-velocity Reynolds number, using a hybrid CPU/GPU computing architecture, and an analysis has been performed of the characteristics of the vortical structures in the wall region of the turbulent channel flow.Turbulent flow structures have been extracted from the simulated flow fields using the ci  or swirling strength criterion, as devised by Zhou et al. [24].
The joint analysis of hairpin vortices and ejection/sweep quadrant events has led to the following conclusions: At t `" 220 (Figure 14) two main primary Ω-shaped vortex filaments are visible at the center of the field (arrows).Right below the head of each filament, the internal space of the structure is occupied by an ejecting isosurface, again showing that, in this initial phase, ejections represent the main mechanism according to which the heads of the structures are raised upwards.Correspondingly (Figure 14b), the process of progressive stretching of the vortex head (actually still a double head) also starts.
In Figure 15, the flow field at t `" 221 is shown.The heads of both structures 1 and 2 continue to raise (being subjected to stretching, Figure 15b) due to the upward-pushing action of the ejections.The sweeping surface mainly maintains its position adjacent to the right side of hairpin 1, causing the compression of the contiguous vortical-structure neck (Figure 15b).Note again that the ejecting isosurface below hairpin 1 is almost equally developed in the streamwise, spanwise and normal-to-the-wall directions (a spheric-like surface), so guaranteeing that the pushing action of the ejecting surface against the internal part of the vortical structure is exerted almost uniformly.It can be also noticed that the presence of sweeping isosurfaces adjacent to hairpin 2 is not so evident as in the case of hairpin 1, and this is the reason why vortex 2 will be destroyed much sooner with respect to vortex 1 (not shown).At t `" 222 (Figure 16), the vortical structures continue their development.The heads of hairpins 1 and 2 continue to raise under the residual influence of the ejection events (the underlying ejecting isosurfaces are now reduced) while a well-defined sweeping isosurface assumes its definite position, adjacent to the external side of head, neck and legs of hairpin 1.The action of compression (Figure 16b) exerted by the sweeping vortex-overlying surface, in particular onto the neck of hairpin 1, becomes more intense.At t `" 223 (Figure 17), the process of evolution of the vortical structures continues in a similar manner, with respect to the previous instants.The action of the ejection events is decreasing, while the sweeping isosurface acts towards the maintenance of the stability of hairpin 1.Also, in this case, the result is hairpin 1 becoming a two-leg, symmetric, and stable vortical structure.

Concluding Remarks
The Direct Numerical Simulation (DNS) of the turbulent flow of incompressible fluid in a plane channel has been executed at three values of the friction-velocity Reynolds number, using a hybrid CPU/GPU computing architecture, and an analysis has been performed of the characteristics of the vortical structures in the wall region of the turbulent channel flow.Turbulent flow structures have been extracted from the simulated flow fields using the λ ci or swirling strength criterion, as devised by Zhou et al. [24].
The joint analysis of hairpin vortices and ejection/sweep quadrant events has led to the following conclusions: (i) the physical condition for the development and subsequent morphological evolution of a stable hairpin-like vortical structure is the occurrence of ejections distributed onto an isosurface almost equally developed along the streamwise, spanwise and normal-to-the-wall directions (a spheric-like isosurface) behind an initially connected Ω-shaped vortex filament, lying near the wall.These ejections actually constitute the physical mechanism according to which the head of the hairpin is raised upward; (ii) the physical condition for the development of a complete and persistent hairpin is the subsequent occurrence of sweeps, as distributed on elongated isosurfaces adjacent to the external sides of the neck and legs of the hairpin.
The sweeps actually constitute the physical mechanism according to which: (ii/a) the legs of the hairpin are stably kept near the wall; (ii/b) the right portion (leg and neck) of the hairpin is characterized by local clockwise particle rotation, the left portion (leg and neck) by counter clockwise local particle rotation.
Author Contributions: All the authors contributed equally to this paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
Roman symbols (upper case) A ij " Bu i {Bx j velocity-gradient tensor C f b " 2τ w {ρu 2 b bulk-velocity friction coefficient Dsc discriminant of characteristic equation D K viscous-diffusion term of turbulent kinetic-energy transport equation F u 1 , F v 1 , F w 1 flatness factors of velocity fluctuations F v 1 | peak peak value of F v 1 K " u 1 i u 1 i {2 mean turbulent kinetic energy L x , L y , L z domain dimensions along x,y,z (h units) L x , L ỳ , L z domain dimensions along x,y,z (wall units) N x , N y , N z number of grid points along x,y,z N tot total number of grid points P K production term of turbulent kinetic-energy transport equation t DB ˇˇsaved actually saved database calculated time (wall units) u i pu, v, wq velocity components along x,y,z u 1 i `u1 , v 1 , w 1 ˘fluctuating-velocity components along x,y,z u 1 rms , v 1 rms , w 1 rms rms velocity fluctuations u 1 rms ˇˇpeak peak value of u Cartesian coordinates (wall units)
2-6 results are presented in terms of turbulence statistics.Figure2reports the values of the Reynolds shear stress (´u 1 v 1 ) in wall coordinates, in a comparison with the data.Computation 2016, 4, 13 8 of 21

Figure 7 .Figure 7 .
Figure 7. Side view of vortical structures in the computing domain: (a)

Figure 7 .Figure 8 .
Figure 7. Side view of vortical structures in the computing domain: (a)

Figure 8 . 21 Figure 9 .
Figure 8.General view of vortical structures on both walls of computing domain at Re τ " 600 (vortical structures are colored with values of local streamwise velocity, reddish indicates high values, greenish indicates low values).Computation 2016, 4, 13 14 of 21

Figure 11 .
Figure 11.λ ci ´isosurface representation of vortical structures in conjunction with quadrant events at t `" 451 (isosurfaces of ejections are shown in red, isosurfaces of sweeps are shown in yellow): (a) vortical structures are shown in cyan; (b) vortical structures are colored with the values of λ r (reddish indicates stretching, bluish indicates compression).Computation 2016, 4, 13 15 of 21

Figure 12 .
Figure 12. λ ci ´isosurface representation of vortical structures in conjunction with quadrant events at t `" 452 (isosurfaces of ejections are shown in red, isosurfaces of sweeps are shown in yellow): (a) vortical structures are shown in cyan; (b) vortical structures are colored with the values of λ r (reddish indicates stretching, bluish indicates compression).

Figure 13 .
Figure 13.λ ci ´isosurface representation of vortical structures in conjunction with quadrant events at t `" 453 (isosurfaces of ejections are shown in red, isosurfaces of sweeps are shown in yellow): (a) vortical structures are shown in cyan; (b) vortical structures are colored with the values of λ r (reddish indicates stretching, bluish indicates compression).

Figure 14 .
Figure 14.λ ci ´isosurface representation of vortical structures in conjunction with quadrant events at t `" 220 (isosurfaces of ejections are shown in red, isosurfaces of sweeps are shown in yellow): (a) vortical structures are shown in cyan; (b) vortical structures are colored with the values of λ r (reddish indicates stretching, bluish indicates compression).

Figure 15 . 21 Figure 14 .
Figure 15.λ ci ´isosurface representation of vortical structures in conjunction with quadrant events at t `" 221 (isosurfaces of ejections are shown in red, isosurfaces of sweeps are shown in yellow): (a) vortical structures are shown in cyan; (b) vortical structures are colored with the values of λ r (reddish indicates stretching, bluish indicates compression).

Figure 16 .
Figure 16.λ ci ´isosurface representation of vortical structures in conjunction with quadrant events at t `" 222 (isosurfaces of ejections are shown in red, isosurfaces of sweeps are shown in yellow): (a) vortical structures are shown in cyan; (b) vortical structures are colored with the values of λ r (reddish indicates stretching, bluish indicates compression).

Figure 17 .
Figure 17.λ ci ´isosurface representation of vortical structures in conjunction with quadrant events at t `" 223 (isosurfaces of ejections are shown in red, isosurfaces of sweeps are shown in yellow): (a) vortical structures are shown in cyan; (b) vortical structures are colored with the values of λ r (reddish indicates stretching, bluish indicates compression).
of u 1 rms ˇˇpeak (wall units) y `ˇ´u 1 v 1 peak y-position of ´u1 v 1 ˇˇp eak (wall units) y `ˇSu 1 peak y-position of S u 1 | peak (wall units) y `ˇFv 1 peak y-position of F v 1 | peak (wall units) Greek symbols (upper case) ∆t `time resolution of calculations (wall units) ∆x `, ∆z `space resolution of calculations along x,z (wall units) ∆y ẁall space resolution of calculations along y at channel wall (wall units) ∆y center space resolution of calculations along y at channel center (wall units) Greek symbols (lower case) ε average rate of dissipation of turbulent kinetic energy per unit mass ε K dissipation term of mean turbulent kinetic-energy transport equation η `Kolmogorov space microscale (wall units) λ eigenvalue λ r real eigenvalue λ cr real part of complex eigenvalue λ ci imaginary part of complex eigenvalue pλ ci q th threshold value of swirling strength ν

Table 2 .
Characteristic parameters of the numerical simulations.

Table 3 .
Runtime of the calculations with different computing-platform configurations (seconds per ∆ t).

Table 3 .
Runtime of the calculations with different computing-platform configurations (seconds per t  ).CPU/GPU Cores 200  τ Re 400  τ Re 600  τ Re

Table 4 .
Characteristic computed quantities of the numerical simulations.