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

Abstract

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.

1. 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.

Accurate DNS calculations of the turbulent channel flow have been carried out by Kim et al. [1], Lyons et al. [2], Antonia et al. [3], Kasagi et al. [4], Rutledge and Sleicher [5], Moser et al. [6], Abe et al. [7], Iwamoto et al. [8], Del Alamo and Jiménez [9], Del Alamo et al. [10], Tanahashi et al. [11], Iwamoto et al. [12], Hoyas and Jiménez [13], Hu et al. [14], Alfonsi and Primavera [15], Lozano-Durán et al. [16], Lozano-Durán and Jiménez [17], Vreman and Kuerten [18,19], Bernardini et al. [20], and Lee and Moser [21], at different values of the Reynolds number (see also at

Table 1. Outline of turbulent-channel-flow DNSs.

Author(s)	Year	Numerical Technique
Kim et al. [1]	1987	Spectral
Lyons et al. [2]	1991	Spectral
Antonia et al. [3]	1992	Spectral
Kasagi et al. [4]	1992	Spectral
Rutledge and Sleicher [5]	1993	Spectral
Moser et al. [6]	1999	Spectral
Abe et al. [7]	2001	Finite Difference
Iwamoto et al. [8]	2002	Spectral
Del Alamo and Jiménez [9]	2003	Spectral
Del Alamo et al. [10]	2004	Spectral
Tanahashi et al. [11]	2004	Spectral-Finite Difference
Iwamoto et al. [12]	2005	Spectral
Hoyas and Jiménez [13]	2006	Spectral-Finite Difference
Hu et al. [14]	2006	Spectral
Alfonsi and Primavera [15]	2007	Spectral-Finite Difference
Lozano-Durán et al. [16]	2012	Spectral
Lozano-Durán and Jiménez [17]	2014	Spectral
Vreman and Kuerten [18]	2014	Spectral
Vreman and Kuerten [19]	2014	Spectral
Bernardini et al. [20]	2014	Finite Difference
Lee and Moser [21]	2015	Spectral

).

The aim of these simulations is mainly that of calculating a given number of time steps of the statistically-steady turbulent flow in the channel, to build up DNS databases, and extracting from the latter useful information for a better comprehension of the wall-turbulence phenomena.

Overall, in the above-mentioned works, the Navier-Stokes equations system is mainly solved within a fractional-step-method framework, in conjunction with Runge-Kutta algorithms for time marching. In particular, in the milestone work of Kim et al. [1], and in Lyons et al. [2], Antonia et al. [3], Kasagi et al. [4], Rutledge and Sleicher [5], Moser et al. [6], Iwamoto et al. [8], Del Alamo and Jiménez [9], Del Alamo et al. [10], Iwamoto et al. [12], Hu et al. [14], Lozano-Durán et al. [16], Lozano-Durán and Jiménez [17], Vreman and Kuerten [18,19], and Lee and Moser [21], the unsteady three-dimensional Navier-Stokes equations are integrated in space by using either the fully spectral Fourier-Chebychev numerical technique originally introduced by Kim and Moin [22], minor variants of the latter, or fully spectral techniques introduced by other authors. In Abe et al. [7] and Bernardini et al. [20], the flow governing equations is integrated by means of a finite-difference algorithm in which a grid-stretching law is inserted orthogonally to the walls. In Tanahashi et al. [11] and Hoyas and Jiménez [13], mixed spectral-high-order finite difference numerical schemes are used. In Alfonsi and Primavera [15], the Navier-Stokes equations in conservative form are integrated by means of the mixed Fourier-finite difference method originally introduced by Alfonsi et al. [23], where a grid-stretching law of hyperbolic-tangent type is inserted along y (the direction orthogonal to the solid walls).

As for velocity boundary conditions, periodic conditions are generally imposed along the streamwise (x) and spanwise (z) directions, in conjunction with no-slip (and impermeability) conditions at the walls, while Neumann conditions are enforced for the pressure. In Table 1, an outline of the above-mentioned works is reported.

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 $R{e}_{\tau }=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 ${\lambda }_{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 $\Omega$-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.

2. 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):

$$\frac{\partial u_i}{\partial t}+\frac{\partial }{\partial x_j}\left(u_i{u}_j\right)=-\frac{\partial p}{\partial x_i}+\frac{1}{R{e}_{\tau }}\frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}$$

(1)

$$\frac{\partial u_i}{\partial x_i}=0$$

(2)

Variables and operators are nondimensionalized by the channel half-height ($h$) for lengths, the wall-shear velocity ($u_{\tau }$) for velocities, the group ($\rho u_{\tau }^2$) for pressure, and ($h/u_{\tau }$) for time, being $R{e}_{\tau }=u_{\tau }h/\nu$ the friction-velocity Reynolds number ($\rho$ is fluid density, $\nu$ 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:

$$\frac{\partial \widehat{u}}{\partial t}+i{k}_x\left(\widehat{u^2}\right)+\frac{\partial \left(\widehat{uv}\right)}{\partial y}+i{k}_z\left(\widehat{uw}\right)+i{k}_x\widehat{p}=\frac{1}{R{e}_{\tau }}\left(\frac{{\partial }^2\widehat{u}}{\partial y^2}-k^2\widehat{u}\right)$$

(3a)

$$\frac{\partial \widehat{v}}{\partial t}+i{k}_x\left(\widehat{vu}\right)+\frac{\partial \left(\widehat{v^2}\right)}{\partial y}+i{k}_z\left(\widehat{vw}\right)+\frac{\partial \widehat{p}}{\partial y}=\frac{1}{R{e}_{\tau }}\left(\frac{{\partial }^2\widehat{v}}{\partial y^2}-k^2\widehat{v}\right)$$

(3b)

$$\frac{\partial \widehat{w}}{\partial t}+i{k}_x\left(\widehat{wu}\right)+\frac{\partial \left(\widehat{wv}\right)}{\partial y}+i{k}_z\left(\widehat{w^2}\right)+i{k}_z\widehat{p}=\frac{1}{R{e}_{\tau }}\left(\frac{{\partial }^2\widehat{w}}{\partial y^2}-k^2\widehat{w}\right)$$

(3c)

$$i{k}_x\widehat{u}+\frac{\partial \widehat{v}}{\partial y}+i{k}_z\widehat{w}=0$$

(4)

where the superscript (^) indicates variables in Fourier space, and $k^2=k_x^2+k_z^2$. 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):

$$y_{str}=PPy+\left(1-PP\right)\left(1-\frac{tanh\left[QQ\left(1-y\right)\right]}{tanh\left[QQ\right]}\right)$$

(5)

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:

$$\frac{{\widehat{u}}_i^{\left(l\right)}-{\widehat{u}}_i^{\left(l-1\right)}}{\Delta t}={\alpha }_{l}D\left({\widehat{u}}_i^{\left(l-2\right)}\right)+{\beta }_{l}D\left({\widehat{u}}_i^{\left(l-1\right)}\right)-{\gamma }_{l}C\left({\widehat{u}}_i^{\left(l-1\right)}\right)-{\zeta }_{l}C\left({\widehat{u}}_i^{\left(l-2\right)}\right)-\left({\alpha }_l+{\beta }_l\right)\frac{\partial {\widehat{p}}^{\left(l\right)}}{\partial x_i}$$

(6)

where l = 1,2,3 denote the Runge-Kutta sub-steps, and where:

$$C\left({\widehat{u}}_i\right)=\frac{\partial }{\partial x_j}\left(u_i{u}_j\right)=i{k}_x\left(\widehat{u_{i}u}\right)+\frac{\partial \left(\widehat{u_{i}v}\right)}{\partial y}+i{k}_z\left(\widehat{u_{i}w}\right)$$

(7)

$$D\left({\widehat{u}}_i\right)=\frac{1}{R{e}_{\tau }}\frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}=\frac{1}{R{e}_{\tau }}\left(\frac{{\partial }^2{\widehat{u}}_i}{\partial y^2}-k^2{\widehat{u}}_i\right)$$

(8)

are the advective- and diffusive terms, respectively. Both terms are treated explicitly and, in Equation (6), ${\alpha }_l,{\beta }_l,{\gamma }_l,{\zeta }_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:

$$\begin{array}{c}\hfill {\alpha }_1=\frac{4}{15},{\alpha }_2=\frac{1}{15},{\alpha }_3=\frac{1}{6};{\beta }_1=\frac{4}{15},{\beta }_2=\frac{1}{15},{\beta }_3=\frac{1}{6};\hfill \\ \hfill {\gamma }_1=\frac{8}{15},{\gamma }_2=\frac{5}{12},{\gamma }_3=\frac{3}{4}\hfill \end{array}$$

(9a)

$${\zeta }_1=0,{\zeta }_2=-\frac{17}{60},{\zeta }_3=-\frac{5}{12};{\displaystyle \sum _{l=1}^3\left({\alpha }_l+{\beta }_l\right)}={\displaystyle \sum _{l=1}^3\left({\gamma }_l+{\zeta }_l\right)=1}$$

(9b)

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 *):

$$\frac{{\widehat{u}}_i^{*\left(l\right)}-{\widehat{u}}_i^{\left(l-1\right)}}{\Delta t}={\alpha }_{l}D\left({\widehat{u}}_i^{\left(l-2\right)}\right)+{\beta }_{l}D\left({\widehat{u}}_i^{\left(l-1\right)}\right)-{\gamma }_{l}C\left({\widehat{u}}_i^{\left(l-1\right)}\right)-{\zeta }_{l}C\left({\widehat{u}}_i^{\left(l-2\right)}\right)$$

(10)

$$\frac{{\widehat{u}}_i^{\left(l\right)}-{\widehat{u}}_i^{*\left(l\right)}}{\Delta t}=-\left({\alpha }_l+{\beta }_l\right)\frac{\partial {\widehat{p}}^{\left(l\right)}}{\partial x_i}$$

(11)

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):

$${\nabla }^2{\widehat{p}}^{\left(l\right)}=\frac{1}{\Delta t\left({\alpha }_l+{\beta }_l\right)}\left(\frac{\partial {\widehat{u}}^{*\left(l\right)}}{\partial x}+\frac{\partial {\widehat{v}}^{*\left(l\right)}}{\partial y}+\frac{\partial {\widehat{w}}^{*\left(l\right)}}{\partial z}\right)$$

(12)

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]).

3. 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:

$$u_i=A_i+A_{ij}{x}_j$$

(13)

($A_{ij}=\partial u_i/\partial x_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:

$${\lambda }^3+P{\lambda }^2+Q\lambda +R=0$$

(14)

where:

$$P=-tr\left(A_{ij}\right);Q=\frac{1}{2}\left\{{\left[tr\left(A_{ij}\right)\right]}^2-tr\left(A_{ij}^2\right)\right\};R=-\mathit{det}\left(A_{ij}\right)$$

(15)

are the scalar invariants of the velocity-gradient tensor (tr is trace, det is determinant). In the case of incompressible flow, $P=0$, and:

$${\lambda }^3+Q\lambda +R=0$$

(16)

$$Q=-\frac{1}{2}tr\left(A_{ij}^2\right)$$

(17)

where the discriminant of the characteristic equation of $A_{ij}$ becomes:

$$Dsc=\frac{R^2}{4}+\frac{Q^3}{27}$$

(18)

When $Dsc>0$, the velocity-gradient tensor has one real eigenvalue (${\lambda }_1={\lambda }_r$), and a pair of complex-conjugate eigenvalues (${\lambda }_2,{\lambda }_3={\lambda }_{cr}\pm i{\lambda }_{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 (${\lambda }_{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 ${\lambda }_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 (${{\lambda }_{ci}|}_{th}$), one can refer to Alfonsi and Primavera [40].

4. 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. Characteristic parameters of the numerical simulations.

Quantities	$R{e}_{\tau }=200$	$R{e}_{\tau }=400$	$R{e}_{\tau }=600$
$L_x$	$4\pi h$	$4\pi h$	$4\pi h$
$L_y$	$2h$	$2h$	$2h$
$L_z$	$2\pi h$	$2\pi h$	$2\pi h$
$L_x^{+}$	2513	5026	7540
$L_y^{+}$	400	800	1200
$L_z^{+}$	1256	2513	3770
$N_x$	256	343	512
$N_y$	181	321	451
$N_z$	256	343	512
$N_{tot}$	$11.9\times {10}^6$	$37.8\times {10}^6$	$118.2\times {10}^6$
$\Delta x^{+}$	9.82	14.65	14.73
$\Delta y_{wall}^{+}$	0.25	0.28	0.30
$\Delta y_{center}^{+}$	3.87	4.36	4.66
$\Delta z^{+}$	4.91	7.33	7.36
${\eta }^{+}$	1.89	2.19	2.42
$\Delta x^{+}/{\eta }^{+}$	5.20	6.69	6.09
$\Delta y_{wall}^{+}/{\eta }^{+}$	0.13	0.13	0.12
$\Delta y_{center}^{+}/{\eta }^{+}$	2.05	1.99	1.93
$\Delta z^{+}/{\eta }^{+}$	2.6	3.35	3.04
$\Delta t^{+}$	0.02	0.04	0.06
$\Delta t$	$1\times {10}^{-4}h/u_{\tau }$	$1\times {10}^{-4}h/u_{\tau }$	$1\times {10}^{-4}h/u_{\tau }$
$t_{DB}^{tot}$	50·$h/u_{\tau }$	50·$h/u_{\tau }$	50·$h/u_{\tau }$
${t_{DB}^{+}|}_{saved}$	500·$\Delta t^{+}$	500·$\Delta t^{+}$	500·$\Delta t^{+}$
${\tau }_{\eta }^{+}$	3.56	4.79	5.87
$\Delta t^{+}/{\tau }_{\eta }^{+}$	0.006	0.008	0.010

.

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 ($\epsilon$).

This method was first introduced by Bakewell and Lumley [41], where in the case of the plane channel, one has:

$$\epsilon \cong \frac{2L_x{L}_z{\tau }_w{u}_b}{2\rho h{L}_x{L}_z}=u_{\tau }^2{u}_b$$

(19)

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. Runtime of the calculations with different computing-platform configurations (seconds per $\Delta t$).

CPU/GPU Cores	$R{e}_{\tau }=200$	$R{e}_{\tau }=400$	$R{e}_{\tau }=600$
1 CPU/240 GPU Cores	0.37	1.71	-
3 CPU/720 GPU Cores	-	-	3.32

). 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 $\Delta t=1\times {10}^{-4}h/u_{\tau }$. The flow fields have been saved every given $\Delta t$ interval, finally giving a 500·$\Delta 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]).

5. Turbulence Statistics

In

Table 4. Characteristic computed quantities of the numerical simulations.

Quantities	$R{e}_{\tau }=200$	$R{e}_{\tau }=400$	$R{e}_{\tau }=600$
$u_{\tau }$ (nominal)	200	400	600
$u_{\tau }$ (present work)	200.23	399.94	600.55
$u_b$ (nominal, after Dean [44])	3390.71	7254.35	11,343.22
$u_b$ (present work)	3197.67	6966.95	11,106.17
$u_c$ (nominal, after Dean [44])	3918.71	8310.35	12,927.22
$u_c$ (present work)	3706.26	7978.80	12,701.63
$u_b/u_{\tau }$ (nominal, after Dean [44])	16.95	18.14	18.91
$u_b/u_{\tau }$ (present work)	16.97	17.42	18.49
$u_c/u_{\tau }$ (nominal, after Dean [44])	19.59	20.78	21.55
$u_c/u_{\tau }$ (present work)	18.51	19.95	21.15
$u_c/u_b$ (nominal, after Dean [44])	1.16	1.15	1.14
$u_c/u_b$ (present work)	1.16	1.14	1.14
$C_{fb}$ (nominal, after Dean [44])	$8.04\times {10}^{-3}$	$6.65\times {10}^{-3}$	$5.95\times {10}^{-3}$
$C_{fb}$ (present work)	$7.86\times {10}^{-3}$	$6.86\times {10}^{-3}$	$6.58\times {10}^{-3}$
$-{\overline{u^{\prime }{v}^{\prime }}|}_{peak}$ (present work)	0.739	0.828	0.866
$-{\overline{u^{\prime }{v}^{\prime }}|}_{peak}$ (from Moser et al. [6])	0.723	0.837	0.864
${y^{+}|}_{peak}^{-\overline{u^{\prime }{v}^{\prime }}}$ (present work)	30.238	40.170	43.938
${y^{+}|}_{peak}^{-\overline{u^{\prime }{v}^{\prime }}}$ (from Moser et al. [6])	30.019	41.882	44.698
${{u^{\prime }}_{rms}|}_{peak}$ (present work)	2.680	2.720	2.751
${{u^{\prime }}_{rms}|}_{peak}$ (from Moser et al. [6])	2.660	2.740	2.770
${y^{+}|}_{peak}^{{u^{\prime }}_{rms}}$ (present work)	14.909	14.199	13.444
${y^{+}|}_{peak}^{{u^{\prime }}_{rms}}$ (from Moser et al. [6])	15.281	14.209	13.268
${S_{u^{\prime }}|}_{peak}$ (present work)	1.003	1.096	1.141
${S_{u^{\prime }}|}_{peak}$ (from Moser et al. [6])	0.922	1.013	1.066
${y^{+}|}_{peak}^{S_{u^{\prime }}}$ (present work)	1.315	1.423	1.504
${y^{+}|}_{peak}^{S_{u^{\prime }}}$ (from Moser et al. [6])	1.339	1.446	1.591
${F_{v^{\prime }}|}_{peak}$ (present work)	26.679 at $y^{+}=0.249$	19.424 at $y^{+}=0.280$	20.882 at $y^{+}=0.298$
$F_{v^{\prime }}|$ (from Moser et al. [6])	26.712 at $y^{+}=0.249$	34.757 at $y^{+}=0.280$	37.653 at $y^{+}=0.298$

and Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6, results are presented in terms of turbulence statistics. Figure 2 reports the values of the Reynolds shear stress ($-\overline{u^{\prime }{v}^{\prime }}$) in wall coordinates, in a comparison with the data.

Figure 3 reports the rms values of the velocity fluctuations (${u^{\prime }}_{rms},{v^{\prime }}_{rms},{w^{\prime }}_{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^{\prime }},S_{v^{\prime }},S_{w^{\prime }}$), and Figure 5 the flatness factors ($F_{u^{\prime }},F_{v^{\prime }},F_{w^{\prime }}$). 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 (${\epsilon }_K$) of the turbulent-kinetic-energy transport equation ($K=\overline{{u^{\prime }}_i{u^{\prime }}_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_{\tau }$,$u_b$,$u_c$, ($u_b/u_{\tau }$), ($u_c/u_{\tau }$), ($u_c/u_b$), $C_{fb}$] 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^{+}=\frac{1}{\kappa }ln{y}^{+}+C$) is also satisfactorily followed, with $\kappa =0.4$ and $C=5.5$ (not shown).

As concerns fluctuating velocity components (Table 4), the peak values of the Reynolds shear stress ($-{\overline{u^{\prime }{v}^{\prime }}|}_{peak}$), those of the rms-fluctuating streamwise velocities (${{u^{\prime }}_{rms}|}_{peak}$), and the fluctuating streamwise velocities skewness factors (${S_{u^{\prime }}|}_{peak}$) [and their positions (${y^{+}|}_{peak}^{-\overline{u^{\prime }{v}^{\prime }}}$, ${y^{+}|}_{peak}^{{u^{\prime }}_{rms}}$, ${y^{+}|}_{peak}^{S_{u^{\prime }}}$)], 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 $R{e}_{\tau }$, as deduced from the analysis of a large number of experimental works (see also Alfonsi [46]):

$${{u^{\prime }}_{rms}|}_{peak}=-0.0000024R{e}_{\tau }+2.70$$

(20)

$${y^{+}|}_{peak}^{{u^{\prime }}_{rms}}=0.00020R{e}_{\tau }+14.6$$

(21)

The values of the skewness factors of the streamwise fluctuations ($S_{u^{\prime }}$, Figure 4) are rather close to zero at the position of rms-fluctuating-streamwise-velocities peak values (${y^{+}|}_{peak}^{{u^{\prime }}_{rms}}$) and, except for these latter positions, are significantly different from the Gaussian values. The values of the flatness factors of the wall-normal fluctuations ($F_{v^{\prime }}$, 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^{\prime }}$ 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.

6. Flow Structures

After the application of the ${\lambda }_{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 $R{e}_{\tau }=200$ and $R{e}_{\tau }=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 $R{e}_{\tau }=600$ with respect to $R{e}_{\tau }=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 $R{e}_{\tau }=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 $R{e}_{\tau }=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 Figure 10, Figure 11, Figure 12 and Figure 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.

The process continues through instants $t^{+}=451$ and $t^{+}=452$ (Figure 11 and Figure 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 (Figure 11b and Figure 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.

Figure 14, Figure 15, Figure 16 and Figure 17 show the evolution in time of a more complex aggregate of vortical structures in a different portion of the computing domain.

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.

7. 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 ${\lambda }_{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.