Independent component analysis of streamwise velocity fluctuations in turbulent channel flows

Ting Wu , Guowei He , *

a The State Key Laboratory of Nonlinear Mechanics, Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China

b School of Engineering Sciences, University of Chinese Academy of Sciences, Beijing 10 0 049, China

Keywords: Independent component analysis Turbulent channel flow Proper orthogonal decomposition Third-order moment Localised modes

ABSTRACT Independent component analysis (ICA) is used to study the multiscale localised modes of streamwise ve- locity fluctuations in turbulent channel flows. ICA aims to decompose signals into independent modes, which may induce spatially localised objects. The height and size are defined to quantify the spatial po- sition and extension of these ICA modes, respectively. In contrast to spatially extended proper orthogonal decomposition (POD) modes, ICA modes are typically localised in space, and the energy of some modes is distributed across the near-wall region. The sizes of ICA modes are multiscale and are approximately pro- portional to their heights. ICA modes can also help to reconstruct the statistics of turbulence, particularly the third-order moment of velocity fluctuations, which is related to the strongest Reynolds shear-stress- producing events. The results reported in this paper indicate that the ICA method may connect statistical descriptions and structural descriptions of turbulence.

Turbulence is composed of flow structures at different scales. Spatially localised structures have important applications in turbu- lent cascades [1] , wall turbulence dynamics [2,3] and aerodynamic noise [4] . However, the structure of turbulence still lacks a com- monly accepted definition, and extracting spatially localised struc- tures at different scales from flow data remains challenging [5] .

One type of flow structure analysis method is vortex identi- fication criteria, such as theQcriterion [6] ,Δcriterion [7] ,λ2criterion [8] , etc. Most of these methods are based on local flow kinematics implied by the velocity gradient tensor and can iden- tify spatially localised structures, such as vortex worms in isotropic turbulence [9] and hairpin vortices in wall turbulence [10] .

Others type of flow structure analysis method includes Fourier analysis, proper orthogonal decomposition (POD) [11] , and wavelet analysis [12] , which decompose a flow into modes at different scales. Fourier analysis is used in homogeneous turbulence, and POD is a generalization of Fourier analysis in inhomogeneous tur- bulence [13] . Both Fourier analysis and POD consider the energy and Reynolds stress of a flow globally so that their modes are spa- tially extended. Wavelet analysis uses localized basis functions and can estimate the multiscale characteristics of the flow at different spatial positions [14] .

In addition, Schmid [15] proposed dynamic mode decompo- sition (DMD) to capture the dynamics of flow fields. The DMD method can extract modes with different temporal properties (i.e., frequencies and growth rates). Thus, the DMD method can effec- tively analyse flows that contain multiple instability mechanisms.

Recently, the empirical mode decomposition (EMD) method [16] has been applied in turbulence research. Agostini and Leschziner [ 17,18 ] used bidimensional EMD (BEMD) to analyse the effect of large-scale structures on near-wall turbulence. EMD par- titions the instantaneous flow into modes (IMFs, intrinsic mode functions) without relying on a priori basis functions; thus, EMD is an adaptive multiscale analysis method.

In addition to these methods, independent component analysis (ICA) has been developed to decompose mixed signals into inde- pendent components and is widely used in data analysis, such as predicting stock market prices [19] and analyzing RNA-sequencing experiments [20] . Carassale [21] used ICA in turbulent flows and described the difference between POD and ICA modes.

In this paper, we use ICA to analyse the localised multiscale modes in turbulent channel flows. We have two motivations for using ICA to analyse turbulence data. On one hand, independence may lead to the spatial locality of modes. If two structures are spa- tially localised (i.e., their intersection in space is relatively small), then they are more likely to be independent. In contrast, if two structures are spatially extended (i.e., their intersection in space is relatively large), then they are less likely to be independent. Jiménez [5] reported a similar explanation while considering the turbulent boundary layer over a wing. Therefore, spatial locality may be achieved through independence of modes. On the other hand, independence itself is also important to describe the statis- tical properties of turbulence. Mouri [22] described the kinetic en- ergy, Reynolds stress and two-point correlations in wall turbulence by assuming the randomness and independence of the attached eddies. Therefore, ICA may be able to connect statistical descrip- tions and structural descriptions of turbulence.

Considering an observed multidimensional variable o , ICA as- sumes that the variable o can be written as [23] :

where o = [o1,o2,···,on]Tis the observed signal, s = [s1,s2,···,sn]Tis a set of mutually independent source sig- nals, and A is a linear mixing matrix, indicating that the observed signalsoiare linear mixtures of independent source signals s . The task of ICA is to find the maximum independent set of source signals s and the mixing matrix A for the given observed signals o . Because both s and A are unknown, we can arbitrarily change the magnitudes and the order of source signals and suitably change the corresponding columns of A to generate the same observed signals. Thus, there are ambiguities in the magnitudes and the order of the independent components [23] . To eliminate these am- biguities, we assume that the source signals exhibit unit variances [24] and thus arrange the source signals by their heights.

In the ICA model, mutual information is used to measure the dependence between random variables. The mutual information betweennrandom variablessi(i= 1,2,···,n) is defined as:

whereH(·)is the differential entropy of the random variable:

wherep(si)is the probability density function (PDF) of random variablesi, andp(s)is the PDF of random vector s . Becausep(s)is ann-dimensional joint PDF,H(s)is difficult to estimate directly. Some pre-processing techniques, such as the whitening transfor- mation, can simplify the estimation of the differential entropyH(s). Mutual information is always nonnegative and equals zero if and only if the variablessiare mutually independent. Therefore, mutual information is a natural measure for independence, and the goal of ICA is to minimize the mutual information between the components (i.e., maximize the independence between the source signals).

There are various ICA algorithms, such as Infomax [25] , FastICA [26] , and MISEP [27] . In this study, we use the FastICA algorithm, which is described in detail by Hyv?rinen and Oja [24] and is im- plemented in Scikit-learn [28]. The Infomax algorithm is also used, and its results are similar to those of FastICA and thus not re- ported.

Considering then-dimensional zero-mean vector o , its covari- ance matrix is:

where the angular bracket 〈·〉 denotes the mathematical expecta- tion, and the superscript “T”indicates the transpose. According to POD, o can be represented by the modal expansion:

whereci(i= 1,2,···,n) are the modal coefficients and the vectors qi(i= 1,2,···,n) are the POD modes that are the eigenvectors of

the covariance matrix:

Among all available linear decompositions, POD is the most ef- ficient in the sense of containing kinetic energy for a given num- ber of modes if o refers to the velocity. Thus, POD is often used to derive reduced-order models. POD modes are orthonormal (i.e.=δij), and the POD coefficients are uncorrelated:

According to Eq. (1) , the ICA model can be rewritten as:

where giis theith column vector of the mixing matrix A . The source signalsiis the coefficient of the ICA mode gi. Due to the independence ofsi, we have:

Thus, the independent variables are uncorrelated. However, uncor- relatedness does not imply independence. Thus, the ICA modes are more advantageous in reconstructing higher-order statistics. Con- sidering the third-order moment as an example, we have= 0 (i/ =j) for the ICA coefficients, but there is no guarantee that= 0 (i/ =j) for the POD coefficients.

Considering the wall-normal profile of streamwise velocity fluc- tuationsu(y;x,z,t)in the turbulent channel flows,x,yandzde- note the coordinates in the streamwise, wall-normal and span- wise directions, respectively; andtdenotes the time. We only discuss the wall-normal modes in this paper and not the three- dimensional modes because the ICA model of three-dimensional modes is difficult to formulate into a well-defined and solvable form. The flow structures in turbulent channels can move ran- domly in the streamwise and spanwise directions together with the downstream convection. If the three-dimensional ICA modes characterize these randomly moving structures, the ICA model can be expressed as:

There arengrid points from the wall to the centre of the chan- nel in the wall-normal direction, and the coordinates of these grid points areyj, wherej= 1,2,···,n. We representu(y;x,z,t)as a mixture of ICA modes:

We first obtain the POD modes of the streamwise velocity fluc- tuations by:

wherehis the half-height of the channel,(y)denotes thejth POD mode,R(y,y′)= 〈u(y;x,z,t)u(y′ ;x,z,t)〉 is the two-point cor- relation in the wall-normal direction, and the angular bracket 〈·〉 denotes the ensemble average, which is performed in time and in the streamwise and spanwise directions due to homogeneity. Then, the velocity field can be represented by the POD modes:

whereojis the coefficient of thejth POD mode. When the dimen- sionnis large, we can neglect some POD modes that have little energy and use an appropriate number of modes to perform the following ICA decomposition, thereby reducing the dimension and making the ICA simpler.

Assuming thatojcan be expressed as an ICA model:

whereojdenotes the observed signal, which can be obtained by Eq. (15) ;si(i= 1,2,···,n) are independent source signals; andAjiis the element in the mixing matrix. Substituting Eq. (16) into Eq. (15) yields:

Thus, we can obtainnICA modes as:

This ICA procedure can be used for other variables, such as pressure and streamwise vorticity, and for all three velocity com- ponents, in which the POD modes of one velocity component must be replaced with the POD modes of all velocity components.

Some ICA modes described by Eq. (18) are not smooth and have small amplitudes, which can be considered as noise. To dis- card these noise modes, we define the following regularization cri- terion based on the local energy ratio:

ρiis the maximum ratio of the energy of theith ICA mode to the local total energy at different wall-normal positions. A smallρiimplies that theith ICA mode can be ignored at all wall-normal positions; then, this mode can be discarded. In this paper, we re- gard the modes ofρi<0.05 as noise and only retain the modes ofρi≥0.05 . We assume thatmmodes remain after the regular- ization, and the indices of these modes arei= 1,2,···,m. The ve- locity field can be exactly reconstructed bynlinearly uncorrelated modes, andmmodes after regularization can only approximately reconstruct the velocity field:

whereuICA(yj;x,z,t)is the approximate velocity field recon- structed bymICA modes ands′iis the modified coefficient. The velocity field of each mode in the channel is:

The reconstructed velocity field can be written as the sum of the contributions of themICA modes:

We letNji=uICAi(yj)and the matrix N =(Nji)n×m; then:

To make the reconstructed velocity fielduICA(yj;x,z,t)a good ap- proximation ofu(yj;x,z,t), we determine the modified coefficient as:

where(Qij)m×n= N+, and “+”in this study indicates the Moore- Penrose pseudoinverse, which provides a least-squares solution to a system of linear equations that lacks a unique solution [29] .

In the next section, the ICA modesfor the wall-normal profile of the streamwise velocity fluctuations will be obtained us- ing the above procedure.

In this study, we use datasets from the direct numerical sim- ulation (DNS) of turbulent channel flows at Reτ≡u(píng)τh/ν= 205 , whereuτis the friction velocity,his the half-height of the chan- nel andνis the kinematic viscosity. The pseudospectral method is used to solve the Navier-Stokes equation, and the 3/2 rule is used to remove aliasing errors. Periodic boundary conditions are used in the streamwise and spanwise directions, and no-slip bound- ary conditions are used at the bottom and top walls. The com- putational domain is 1.3 πh×2h×0.35 πh, and the grid numbers 64 ×129 ×32 are used in the streamwise (x), wall-normal (y), and spanwise (z) directions. The time step to advance the Navier-Stokes equations is considered to beΔt+= 0.012 , and “+”indicates nor- malization with viscous scales. The Navier-Stokes solver and the datasets used in this study have been validated in previous studies [30] .

One hundred snapshots of the instantaneous streamwise veloc- ity fluctuations are used. Because we only investigate the wall- normal modes from the wall to the centre of the channel, not the modes between two walls, the total number of samples for the ICA procedure is 409,600 (64 ×32 ×100 ×2). The velocity at the wall is zero, and the grid point at the wall is therefore not considered in the POD and ICA; there are also 64 grid points in the wall-normal direction. Therefore, 64 POD modes and 64 ICA modes are obtained according to the procedure. With regularization, 33 ICA modes re- main. To test whether the 33 ICA modes can reconstruct the ve- locity field well, we define the relative error of the reconstructed velocity field:

Results show that the relative errorε(y)is less than 0.2% at all positions and less than 0.1% at most positions; thus, the 33 ICA modes can be used to reconstruct the velocity field well.

Figure 1 a plots the velocity profile of four POD modes withi= 1,3,6,9 , and Fig. 1 b plots the velocity profile of four ICA modes withi= 5,11,17,24 . The POD modes are spatially extended and oscillatory throughout the region, similar to the trigonomet- ric functions in Fourier analysis. While the ICA modes are spa- tially localised, particularly for those modes whose peaks are near the wall, the amplitudes of these modes asymptotically tend to- ward zero at the centre of the channel. In addition, there are some modes whose peaks are near the centreline, and these modes con- tribute to the fluctuations at the centre of the channel. Figure 1 b implies that spatial locality may indeed be induced by indepen- dence, which can be described as follows. If two structures are localised in space (i.e., their intersection is relatively small), then they are more likely to be independent. In contrast, if two struc- tures are spatially extended (i.e., their intersection in space is rela- tively large), then they are less likely to be independent. Therefore, the ICA method provides a path to obtain spatially localised struc- tures.

Fig. 1. (a) Velocity profile of four POD modes with i = 1 , 3 , 6 , 9 . (b) Velocity profile of eight ICA modes with i = 5 , 11 , 17 , 24 (solid lines) and i = 30 , 31 , 32 , 33 (dashed lines).

Fig. 2. Schematic of the height and size of the ICA mode.

Fig. 3. (a) Sizes of the first 33 POD modes vary with height. The blue dashed line indicates that the size is h . (b) Sizes of the ICA modes vary with height. The blue dashed line indicates that the size is proportional to the height, and the black dotted circle indicates the ICA mode with i = 5 .

Fig. 4. Normalized velocity profile of the ICA modes with i = 10 , 11 , ···, 16 .

Fig. 5. Instantaneous isosurfaces of the streamwise velocity fluctuations of the ICA mode with i = 5 . The blue coloured objects are low-velocity streaks, = ?0 . 5 . The red coloured objects are high-velocity streaks, = 0 . 5 .

Fig. 6. (a) Velocity profile of the ICA mode with i = 5 . (b) Spanwise correlation of the ICA mode with i = 5 . The blue dashed line indicates the local maximum correlation, and the green dash-dotted line indicates the minimum correlation.

To quantify the spatial extension of a mode, the heightγand sizedare defined as two characteristic length scales. The height of a modeui(y)is defined as:

whereγirepresents the energy-weighted wall distance of theith mode. The size of a modeui(y)is defined as:

wheredirepresents the energy-weighted spatial extension of theith mode. If the energy of a mode is uniformly distributed in the region [0,h] , the size of the mode is calculated ashaccording to Eq. (27) . Figure 2 shows sketches of the height and size of a mode. Both the height and size of the mode shown in Fig. 2 b are larger than those shown in Fig. 2 a. Although the height of the mode in Fig. 2 c is larger than that in Fig. 2 a, the two modes have the same size. As mentioned before, the order of the independent compo- nents cannot be determined by the ICA procedure [23] ; thus, we sorted the ICA modes according to their heights from small to large.

Figure 3 a plots the sizes of the first 33 POD modes with the heights. The sizes of most POD modes are h, which means that the energy of most POD modes is distributed throughout the re- gion [0,h] ; thus, the POD modes approximate the velocity field in a global sense. Figure 3 b plots the sizes of ICA modes with the heights. The sizes of most ICA modes are smaller than those of POD modes, indicating the spatial locality of ICA modes. In addi- tion, the sizes of ICA modes exhibit a multiscale nature and are approximately proportional to the heights, particularly the modes withi= 10,11,···,16 . We define the normalized velocity profile of the ICA mode:

Figure 4 plots the normalized velocity profile of the ICA modes withi= 10,11,···,16 , where the wall distanceyis normalized by the height of each mode. A good collapse is shown in Fig. 4 , high- lighting the similarity of these ICA modes. This collapse implies that the height is the characteristic length scale of an ICA mode.

The mode withi= 5 is special and is marked by the black dot- ted circle in Fig. 3 b; its size is much larger than other modes at the height≈13 in the buffer layer. Figure 5 plots the instantaneous isosurfaces of the streamwise velocity fluctuations of the ICA mode withi= 5 . Streamwise elongated streaks are also shown. Figure 6 a plots the velocity profile of the ICA mode withi= 5 . The positions with large amplitudes of this mode are in the region ofy+<30 ; thus, this mode primarily exists in the buffer layer. We then calcu-

Fig. 7. Instantaneous isosurfaces of the streamwise velocity fluctuations of four ICA modes. The blue coloured isosurfaces,= ?0 . 5 . The red-coloured isosur- faces,= 0 . 5 . (a) i = 1 ; (b) i = 7 ; (c) i = 30 ; (d) i = 33 .

Fig. 8. (a) Sum of the second-order moments of all ICA modes compared with the DNS result of streamwise velocity fluctuations. (b) Sum of the second-order moments of ICA modes with i = 1 , 2 , ···, 16 and with other ICA modes.

late the spanwise correlation of this mode:

Figure 6 b plots the spanwise correlation of the ICA mode withi= 5 . The mean spanwise spacing between the adjacent low- and high-speed streaks isλ+z≈55 according to the minimum of the spanwise correlation, which is consistent withλ+z≈50 in Kim et al. [31] . The mean spanwise spacing between two adjacent low- speed or high-speed streaks isλ+z≈110 according to the local maximum of the spanwise correlation, which is consistent with≈100 in Smits et al. [32] .

Figure 7 plots the instantaneous isosurfaces of two near-wall ICA modes withi= 1,7 and two outer ICA modes withi= 30,33 . The isosurfaces ofi= 1 andi= 30 are disrupted and irregular, which are markedly different from the streamwise streaks of thei= 5 mode. The isosurfaces ofi= 7 andi= 33 are elongated in the streamwise direction, similar to the streamwise streaks, but their spanwise spacing is larger than that of thei= 5 mode.

If the ICA modes are strictly independent, we have the follow- ing equations for turbulence statistics:

Figure 8 a plots the sum of the second-order moments of all ICA modes compared with the DNS result of streamwise veloc- ity fluctuations. The two results are similar, which comes from the approximate independence between these modes. Figure 8 b plots the sum of the second-order moments of ICA modes withi= 1,2,···,16 and with other ICA modes. The energy of modes withi= 1,2,···,16 is distributed across the near-wall region. Figure 9 a plots the DNS result of the wall-normal correlation of stream- wise velocity fluctuations. Figure 9 b plots the sum of the wall- normal correlation of all ICA modes. These two results are sim- ilar. Figure 9 c plots the sum of the wall-normal correlation of ICA modes withi= 1,2,···,16 , and Fig. 9 d plots the result with other ICA modes. The wall-normal correlation of ICA modes withi= 1,2,···,16 is distributed in the near-wall region, which further shows the spatial locality of the ICA modes.

Fig. 9. (a) DNS result of the wall-normal correlation of streamwise velocity fluctuations. (b) Sum of the wall-normal correlation of all ICA modes. (c) Sum of the wall-normal correlation of ICA modes with i = 1 , 2 , ···, 16 . (d) Sum of the wall-normal correlation of ICA modes with i = 17 , 18 , ···, 33 .

Fig. 10. (a) Sum of the third-order moments of all ICA modes and that of all POD modes compared with the DNS result. (b) Sum of the third-order moments of ICA modes with i = 1 , 2 , ···, 5 and with other ICA modes.

Figure 10 a plots the sum of the third-order moments of all ICA modes and that of all POD modes compared with the DNS re- sult. The DNS third-order moment is positive in the near-wall re- gion and negative in the other region. The behaviour of the third- order moment is expected from the quadrant analysis and the most violent Reynolds shear-stress-producing events [31] ; thus, the strongest Reynolds shear-stress-producing events are the ejection events (u′<0) fory+>12 ; sweep events also occur (u′>0) fory+<12 . The sum of the third-order moments of all ICA modes has this feature, while that of all POD modes does not. The sum of the third-order moments of all POD modes is always negative through- out the region. Therefore, compared with POD, the ICA method is better for the higher-order statistics of turbulence. As mentioned before, the obtained ICA modes are not completely independent; thus, the sum of the third-order moments of all ICA modes is not exactly equal to the DNS result. Figure 10 b plots the sum of the third-order moments of ICA modes withi= 1,2,···,5 and with other modes. The modes withi= 1,2,···,5 contribute to the pos- itive third-order moment in the near-wall region, and the other modes contribute to the negative third-order moment.

In this paper, the ICA method is used to decompose the wall- normal profile of the streamwise velocity fluctuations in turbulent channel flows. ICA aims to decompose the signals into indepen- dent components, and the spatial locality may be induced indepen- dently. Therefore, spatially localised structures in turbulence may be obtained by the ICA method.

Using the DNS data of the turbulent channel flows atReτ= 205, we find that ICA modes are indeed spatially localised, and the second-order moments and the wall-normal correlation of some modes are localised near the wall and tend toward zero at the centre of the channel. In addition, the sizes of the ICA modes ex- hibit a multiscale nature and are approximately proportional to the heights. The size and height are well-defined length scales in ICA modes that neither extend to the entire region nor shrink to a point. Thus, ICA can be used as a data-driven method for lo- calised mode decomposition. ICA modes can also reconstruct the third-order moment of the velocity fluctuations well. The sum of the third-order moments of all ICA modes is positive in the near- wall region and negative in the other region, which agrees with the DNS result and the mechanism of the most violent Reynolds shear- stress-producing events. As a comparison, the sum of the third- order moments of all POD modes is always negative throughout the region.

As an application of the ICA method in turbulence, this study shows that the ICA may be able to connect statistical descriptions and structural descriptions of turbulence. Future work should anal- yse the exact relationship between ICA modes and the attached eddy model of wall turbulence in detail.

Declaration of Competing Interest

The authors declare that they have no known competing finan- cial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work is supported by NSFC Basic Science Center Program for “Multiscale Problems in Nonlinear Mechanics”(No. 11988102) and National Natural Science Foundation of China (Nos. 12002344 , 11232011 and 11572331). The authors would like to acknowl- edge the support from China Postdoctoral Science Foundation (No. 2020M670478), the Strategic Priority Research Program (No. XDB22040104) and the Key Research Program of Frontier Sciences of the Chinese Academy of Sciences (No. QYZDJ-SSW-SYS002).