1 Introduction

Turbulence is everywhere, controlling the drag on cars, airplanes, and trains whilst affecting weather through its influence on large-scale atmospheric and oceanic flows. Even solar flares are a manifestation of turbulence as they are triggered by vigorous motions on the surface of the Sun. It is easy to be intrigued by a subject so pervasive in daily life.

Turbulence is one of the most difficult domains in nonlinear physics owing to its complex, daunting mathematical descriptions, profound difficulties of its inherent instabilities, and its chaotic processes. Until today, there has been no obvious breakthrough in the existence and uniqueness of turbulence solutions in mathematics. However, Chinese scholars have made some promising achievements in turbulence studies over the past 30 years. This article highlights some of the major developments that Chinese scholars have pioneered in the fields of turbulence theory, transition, and calculations, as well as in experiment and measurement techniques.

The paper is organized as follows. Section 2 reviews turbulence theory. Section 3 analyzes turbulence transition. Section 4 examines turbulence calculations. Section 5 evaluates experiment and measurement techniques. Finally, Section 6 concludes with future perspectives.

2 Turbulence theory 2.1 Zhou's turbulence theory

Zhou's (or Chou's) turbulence theory^[1-2] mainly consists of two parts, which can be applied≠wpage to turbulence models with arbitrary high-order statistics and to the statistical theory of turbulent vortex structures by first solving and then averaging. The first part of the theory was based on two studies conducted from 1940 to 1945, while the second one dealt with systematic studies during the 1950s and 1980s. It is worth noting that these two parts led to a number of advances both before and after Zhou's death in 1993, and constituted an important achievement in the field, namely, the basic theory of turbulence in China over the past 30 years.

Zhou's two 1940s studies proposed a step-by-step approximation method to solve the velocity correlation equations at various orders. Using this method and the quasi-similarity condition in the 1980s, one of Zhou's students, Chen, completed the calculations for the channel flow and the plane wake flow, thereby obtaining the result^[3-4] for the third-order correlation. This result was consistent with experimentation. In 1988, another Zhou's student, Meng, completed the calculation for a plane mixing layer and again obtained results that were consistent with experimentation^[5]. In 1991, Lin (also a student of Zhou) determined the relationship between the free shear turbulence and fourth-order correlation in a two-dimensional mixed layer, a turbulent wake, and a turbulent jet^[6]. Comparison of these results with those of previous studies showed that Zhou's theory was successful for these flows. Lin also obtained results for turbulence in a channel using the quasi-similarity theory whilst retaining the viscous term^[7].

Although credible results can be obtained, the odd-order truncated step-by-step approach is complicated and requires certain assumptions to be made at each level, which greatly restricts the theory's extension. Zhou had been thinking about this problem for some time, and in 1988 he proposed a step-by-step iterative method to simultaneously solve the equations for mean flow and turbulent velocity fluctuations^[8].

The following provides a basic understanding of this step-by-step iteration method. First, in the averaged equation of motion and in the second-order velocity correlation equation, the third-order velocity correlation is ignored. The results for the Reynolds stress and the mean velocity are obtained as zeroth-order approximations. These zeroth-order approximation results are then substituted into the turbulent fluctuation equation, which is solved for the fluctuating velocity. The Reynolds stress is then recalculated using this fluctuating velocity. Following this, the new Reynolds stress result is substituted into the averaged equation of motion, which is solved for the mean velocity. Finally, the new Reynolds stress and new mean velocity are taken as the first-order approximations. By repeating this procedure, it is possible to calculate the higher-order approximations of the mean velocity, as well as the velocity dependence. In particular, once the turbulent fluctuating velocity has been obtained, any higher-order velocity correlation can be calculated in principle. It is clear that the calculations in this step-by-step iteration method are simple compared to those in the previous step-by-step approximation method. Another attractive feature of this method is that the fast variables, such as the average velocity of the turbulent flow, the fluctuating kinetic energy, and the fluctuating velocity, are separated because of the assumption of the quasi-similarity theory and fast Fourier transform (FFT), which can be used to solve the fluctuating field in micro-scale normalized regions under a local reference system.

Using this method, Zhou, Huang, and their students conducted numerical calculations of various typical flow fields. Assuming quasi-similarity, Wang et al.^[9] used the FFT to solve the fluctuating field of homogeneous isotropic turbulence, and re-obtained the second- and third-order velocity correlation functions from late decay to a low Reynolds number. For a plane wake, Lin et al.^[10] used the calculated Reynolds stress equivalent as the zeroth-order approximation and applied the FFT to solve the fluctuating velocity field in a local reference system. Thus, they obtained satisfactory results for the fourth-order velocity correlation without any iteration. Furthermore, Lin et al.^[11] used this method to calculate high-order velocity correlations in a channel flow, and again obtained satisfactory results. Meng et al.^[12] obtained satisfactory high-order velocity correlation results for plane mixing layers in the doctoral dissertation. Fan et al.^[13] presented the results of calculations for a plane turbulent jet.

There have been recent developments in the statistical theory of turbulent vortex structures. Tsai and Lin^[14] proposed a new method, which combined Zhou's theory of vortex structures and Kraichnan's direct interaction theory. Considering that Kraichnan's direct interaction theory applies to waves rather than vortices, they first defined a single turbulent vortex as a moving fluid group with a self-similar structure. Following this, they solved the free vortex and used methods derived from quantum field theory to describe the vortex in terms of creation and annihilation operators. The "Schrödinger equation" for vortex interaction was constructed using the interaction term in the Navier-Stokes equation to describe the state change.

In the early 1990s, Huang and his collaborators observed that the spherical vortices of Beltrami flow, obtained as a constant solution of the Euler equation or from decay of the Navier-Stokes equation, can produce complex chaotic phenomena by superposition. In addition, they found an interesting double-ring vortex structure^[15]. They classified the exact vortex solutions described in the literature and found that a straight spiral vortex can be superimposed on an arbitrary coaxial Beltrami flow (again obtained as a constant solution of the Euler equation or from decay of the Navier-Stokes equation), citing Moffatt's well-known spherical vortex as a special case^[16]. Huang^[17] and Huang and Luo^[18] also took a non-axisymmetric vortex as a turbulence element and introduced a general mathematical formula for the direction average, thereby improving the work of Saffman and Pullin^[19]. As the longest-standing collaborator in Zhou's turbulence studies, Huang and Hu^[20] made a systematic summary of the statistical theory and recent developments in the turbulent vortex structure in 2002. This summary was based on Zhou's theory. The closure of the Reynolds equation obtained from the Navier-Stokes equation and Reynolds equation was also briefly discussed by Lee and Lian^[21].

In the 1990s, Zhou's goal of solving the averaged equation of motion and the equation for turbulent fluctuations was finally realized. The statistical theory of the turbulent vortex structure, based on solving following averaging, led to the foundation of a unique approach to turbulence research in China.

For important cases of high-order statistics that arise in engineering, the current popular mode theory cannot go beyond three orders of correlation. In this situation, Zhou's turbulence theory can play an important role. However, how this approach can be used in practice for more challenging turbulent flows (such as wall turbulence and separation flows) remains debatable. By mining direct numerical simulation (DNS) databases as deeply as possible, it is hoped that a combination of Zhou's theory and large eddy simulation (LES) may produce more favorable results.

2.2 Physical mechanisms of homogeneous isotropic turbulence

The study of the physical mechanisms of homogeneous isotropic turbulence is the most important and fundamental aspect of turbulence physics. Revealing and depicting these mechanisms are also the starting point of turbulence modeling.

Intermittency is an important characteristic of turbulent fields (in contrast to random fields) and has always been at the core of turbulence research. Studies that deal with intermittency in turbulence have concentrated on the singular scaling law of turbulence, with the degree to which the scaling index deviates from that of the Kolmogorov K41 theory used as a criterion for intermittency. In the history of theoretical turbulence research, several phenomenological models have been proposed to describe the singular scaling law of turbulence. However, no model to date provides complete information on turbulent structures. There are many reasons for this, one of which is a lack of effective means to identify such structures.

Chen and his collaborators improved the global zero-crossing and peak-valley recognition methods and proposed a local zero-crossing method, which used the number of crossing points as a marker for structural classification. Through analysis of incompressible turbulence, compressible turbulence, and transport of a passive scalar field, they showed that the probability of structural classes, classified by the local zero-crossing method, is scale-invariant in the inertial region. Moreover, the strongest intermittent structure is determined by the structural fluctuations of a single zero-crossing point. Thus, for the first time, structural information was introduced into turbulent batch models. From the analysis of the dynamics of vortex filaments, the above phenomena could be explained. Furthermore, the relations among the law of motion of a vortex filament, the structure of a single zero-crossing point, and the correlations of the strongest intermittent structure were revealed^[22].

In recent years, high-Mach-number turbulence has attracted much attention owing to the needs of the aerospace industry and high-energy-density physics research. The study of multi-scale and cascade processes has always been at the core of turbulence theory. However, there have been few multi-scale analyses from the point of view of multi-process decomposition, which explains the gaps in the understanding of the physical mechanisms of turbulence at a high Mach number. Research by Chen's team clearly showed that the multi-scale analysis must be based on the multi-process decomposition^[23-24]. Thus, a framework for understanding multi-scaling behavior and energy level series analysis in shear/bulging processes was proposed. This proved to be a feasible starting point for the theoretical study of high-Mach number turbulence. This study clearly revealed the physical essence of compressible turbulence at a high Mach number for the first time. It is clear that the modeling of high-Mach-number compressible turbulence should be based on the modeling of both the shear and bulging processes coupled together.

2.3 Spatiotemporal correlation model

Spatiotemporal correlation is the main method for studying the dynamic coupling of the spatial and temporal scales of motion in turbulence^[25-26]. Moreover, it determines the relationship between a transient turbulence disturbance in one position and an instantaneous disturbance in another position, thereby describing the behavior of turbulent disturbances on both space and time scales. Spatiotemporal correlation models of turbulence include the random scanning model and the local strain model of isotropic and uniform turbulence, the elliptical model of Taylor's freezing and turbulent shear flow, and the linear wave propagation and compressible turbulence models. He and Zhang^[27], He et al.^[28], Zhao and He^[29], and He et al.^[30] proposed an elliptical approximation model of Lagrange space-time correlations that included flow deformation factors. These factors modified Taylor's freezing hypothesis, which essentially represents degeneration of macro-scale vortices in turbulent shear flow^[27-30]. The elliptical approximation model was verified both numerically and experimentally, for example, via DNS of turbulent channel flow at a low Reynolds number^{[29, 31]} and experimental results on the Rayleigh-Bénard convection^[30]. The Taylor freezing hypothesis, as the first spatiotemporal correlation model, is widely used in all fields of turbulence research. For example, the spatial energy spectrum can be obtained from continuous-time measurements at a single point in space using this hypothesis. Indeed, many models of turbulent noise are dependent on Taylor's freezing hypothesis. The elliptic model is hugely significant to these related fields.

Spatiotemporal correlations can be used to provide the time scale and precise subgrid-scale turbulence models for LES^[32-33]. In contrast, the simultaneous-time precise subgrid-scale LES method can also be used to study spatiotemporal correlations^[34-35]. The turbulent noise problem^[36] and the particle-loaded turbulence correlation problem^[37] are based on the Eulerian spatiotemporal correlation of Lighthill's acoustic analogy and the Lagrangian time correlation. Hence, the accurate prediction of spatiotemporal correlations is the key to the study of these related problems.

Turbulent coherent structures play an important role in the study of turbulence as their dynamics depend on their specific characteristics. The various methods that identify turbulent structures based on the iso-surfaces of various physical quantities in the instantaneous flow field have clear defects. Thus, Yang et al.^[38] developed and applied a stable and topologically preserved particle inverse tracking method to perform multi-scale geometric analysis of isotropic turbulence. Later, a new definition^[39-40] was proposed wherein the vortex surface was defined by the equivalent surface of the vortex surface field. An associated vortex surface evolution equation was then obtained. The related phenomena of eddy dynamics, including vortex reconnection^[41-42], can thereby be explained. Yang and Pullin^[43-45] and Zhao et al.^[46-47] also studied and compared the Lagrangian and Eulerian structures in the channel flow using computational methods.

2.4 Qian's method

A number of researchers have sought to improve on the Kolmogorov^[48-49] (K41) analysis by matching solutions valid in the inertial subrange to those valid in the viscous dissipation range, i.e., considering finite Reynolds number effects on K41. Among them is Qian, who has published 27 single-authored international peer-reviewed papers dealing with aspects of the scaling properties of homogeneous, isotropic turbulence since 1983. An overview of "Qian's method" (see Refs. [50]-[52]) appears in McComb (2014, pp. 16-170)^[53]. Notably, Qian has investigated whether the second-order exponent corresponds to normal Kolmogorov scaling or anomalous scaling. Qian's method makes use of exact relationships incorporated with well-established data correlations available in the literature to extract as much physics as possible from the available experimental results.

2.5 Quantifying wall turbulence via a symmetry approach — a Lie-group theory

She et al.^[54], Chen et al.^[55], and Chen and She^[56] recently revealed the first unified analytic description of the non-uniform distribution of mean velocity and turbulence intensities across the entire domain of turbulent boundary layers for all Reynolds numbers. This had been a century-old goal dating back to Prandtl's pioneering work in 1904. This Lie group inspired the symmetry approach. Their method is particularly interesting in terms of methodology, as they proposed a universal dilation-symmetry-breaking form, which, as they argued, is necessarily present in non-equilibrium and complex systems (such as turbulence). In other words, they extended the classical symmetry analysis, which usually assumes symmetry preservation to deal with a non-equilibrium system with symmetry-breaking. Indeed, this method was successfully applied to describe a difficult problem involving rough pipe^[57], for which the roughness near the wall added a new dimension to the complexity and compressible flow with density and temperature variations^[58-59].

3 Turbulence transition

Luo et al.^[60] studied the effects of disturbances on both the basic flow profile and the transition process in the plane channel flow. They found that the growth of a disturbance can correct the basic flow, leading to a rapid increase in the range of wavenumbers of the unstable disturbance and a rapid increase in the growth rate. This causes the "breakdown" process. The results of this study revealed the importance of corrections to the basic flow and of the characteristics of nonlinear interaction between disturbances to the transition process. Zhou's research group then studied the systematic mechanisms. Tang et al.^[61] presented results for an incompressible boundary layer. Huang et al.^[62] performed a time model calculation and simulation for the mechanism of the laminar flow catastrophe in a supersonic boundary layer transition. Similar results were presented by Cao et al.^[63]. For a supersonic plate boundary layer, Dong and Zhou^[64] observed the same mechanism in the bypass transition, which was caused by large fluctuations of the incoming flow. These results confirmed that basic flow correction can promote the transition process. Zhang and Zhou^[65] studied the boundary layer receptivity of vortex perturbations in the free stream, while Gao et al.^[66] studied the receptivity of fast and slow sound waves in a supersonic boundary layer. The latter showed that a disturbance in the free stream can stimulate an unstable disturbance in the boundary layer through the outer edge. This provided a new way to understand the receptivity of the free stream disturbance in the boundary layer. Ren and Fu^[67] studied unsteady disturbance spectrum in a boundary layer with a curved surface at Mach 4.5. They found that a flow vortex/entropy disturbance could be transformed into a quasi-steady Görtler vortex, while a slow sound wave could be converted into an unsteady Görtler vortex. Jiang and Lee^[68] reviewed the receptivity of hypersonic boundary layers. In addition, Ren and Fu^[69] studied secondary instability of Görtler vortices and found the anti-symmetric mode to be more unstable than the symmetric mode at medium and high Mach numbers. Ren et al.^[70] used the parabolic stability equation (PSE) method to study the effects of streaks and Görtler vortices on the first and second modes, and found that the proper amplitude of the streak and the Görtler vortex can inhibit the growth of the first and second modes, thereby postponing the transition.

The e^N method is considered to be a transition prediction method with a relatively firm theoretical basis and has mostly been used in the aviation industry. However, this method has two weaknesses. The receptivity problem is not considered, and the geometry of the transition position is determined entirely by experience. A study by Su and Zhou^[71-72] showed that the traditional e^N method has problems in predicting transitions. Thus, they considered receptivity and proposed an improved e^N method. The results of this method are consistent with both numerical simulation results and those of wind tunnel experiments^[73-74].

Li et al.^[75] studied the natural transition in a blunt cone boundary layer at Mach 6 and at a small angle of attack using DNS. They found that low-frequency disturbances are more significant in the late stages of transition than in the second mode (see Fig. 1). Liang et al.^[76] used DNS to study the effect of wall temperature on stability in a blunt cone boundary layer at Mach 8. They found that reducing the wall temperature can significantly increase second-mode instability but can inhibit first-mode instability. In terms of the relationship between the characteristics of the flow disturbance and the transition position, Su and Zhou^[77] studied the changes in the transition position caused by changes in the disturbance amplitude over time. They found that the change in the transition position lags when the amplitude of the upstream disturbance changes over time. Both the turbulent and laminar flows appear alternately.

Liu and Luo^[78] studied the effects of different frequency disturbances on the transition position in a hypersonic blunt cone boundary layer at a small angle of attack. They found that the transition position depends on the distribution of the perturbation amplitude with frequency, which could explain the shift or forward movement of the transition position in a hypersonic blunt cone boundary layer.

These studies illustrate how the transition process depends strongly on the characteristics of both flow disturbances and wall disturbance, and reveals the scope for a wide-ranging investigation of the entire transition process.

The PSE approach is an efficient method used to calculate the evolution of perturbations in a boundary layer. Li and Luo^[79] studied these transitions in an incompressible boundary layer using this method. Zhang and Zhou^[80-81] applied the PSE method to study transitions in a compressible and supersonic boundary layer. All these investigations produced favorable results. In addition, rapid calculation methods combining the PSE and DNS were proposed to study both transitional and fully developed turbulence, which allowed for effective prediction of the evolution of perturbations and the transition to turbulence in supersonic and hypersonic boundary layers. Zhang and Su^[82] proposed an extension of the PSE method, and compared their results with those of the traditional PSE. They found little difference between them. Yu et al.^[83] proposed a stability analysis method called the EPSE. This method is also based on an eigenvalue problem; however, it can handle the stability analysis of nonparallel basic flows, thereby widening its range of application. Based on the PSE theory, Luo et al.^[84] developed a highly efficient numerical method to solve the nonlinear evolution of finite-amplitude Tollmien-Schlichting (T-S) waves, and analyzed the new modes produced by the nonlinear interaction of disturbance waves. Using the linear stability theory and multiple scales method, Tang et al.^[85] showed the influence of nonparallel flow on boundary layer stability. The stability problem of three-dimensional T-S waves in a non-parallel boundary layer was analyzed using the PSE method with the original-variable form of three-dimensional disturbance waves^[86]. Guo et al.^[87] analyzed and compared the multiple unstable modes of inviscid and viscous flows at a high Mach number through a numerical approach. Finally, Liu et al.^[88] studied the evolution of C-type instability in a nonparallel boundary layer using the PSE method. A systematic summary of the transition process, transition mechanism, and transition control in boundary layers was given in Ref. [89].

The research team headed by Lee at Peking University has also completed much experimental and theoretical work on hypersonic boundary layer transitions^{[68, 90]}. Zhang et al.^[91] and Zhang and Lee^[92] were the first in China to publish quiet wind tunnel experimental results. In particular, they used a CO₂ Rayleigh scattering technique to obtain clear images of the hypersonic laminar-turbulent transition, as shown in Fig. 2. Tang et al.^[85] studied the effect of two-dimensional roughness elements on the evolution of second-mode waves by means of the partially coherent beam (PCB) technique, particle image velocimetry (PIV), and other experimental methods. Zhang et al.^[93] and Zhu et al.^[94] successfully applied PIV to hypersonic boundary layer measurement and found that the second-mode waves first increased and then decayed, whilst the low-frequency disturbance amplitudes were large. Zhu et al.^[95-96] performed experiments and numerical simulations at Mach 6 in a cone boundary layer and studied aerodynamic heating, which was highly praised by Sun and Oran^[97]. They found that the bulging effect of second modes during the transition process can cause large rises in temperature, which can exceed those of turbulent-laminar transition effects. Chen et al.^[98] used the PSE, the two-time instability theory, and an energy analysis method to study the interaction between second and first modes in a conical boundary layer at Mach 6. They found that the second mode played the role of a catalyst, causing a rapid increase in the first mode.

4 Turbulence calculations

In recent years, the computational speed of computers has rapidly increased, and the development of supercomputers and parallel computing methods have provided better conditions for numerical simulations of turbulence. DNS provides direct solutions for the Navier-Stokes equation of control flow, meaning that the information on all turbulence scales can be obtained. DNS is one of the key methods currently used for the study of turbulence. The Navier-Stokes equations for both incompressible and compressible fluids were studied by Li et al.^[99-100] using DNS. Li et al.^[101] used DNS to study coherent structures in a Mach 6 flat-plate boundary layer. They found that quasi-directional vortices dominated at high Mach numbers, whereas entrapment vortices rarely appeared^[102]. In particular, OpenCFD, an open source compressible high-precision computational fluid dynamics software that Li developed, has been widely recognized and applied. Deng and Zhang^[103] developed a weighted compact nonlinear scheme (WCNS). Ma and Fu^[104] developed a group velocity control (GVC) scheme. Both schemes showed good performance. Dong et al.^[105], Dong and Luo^[106], Dong and Zhou^[107], Dong and Luo^[108], and Dong^[109] conducted a series of turbulent boundary layer studies under supersonic and hypersonic conditions using DNS and the Baldwin Lomax model. Finally, Dong and Zhou^[110] accomplished corrections to the eddy thermal conductivity near the peak of turbulent kinetic energy.

In addition, the gas kinetic method is a new approach based on the kinetic theory of gases. Li and Fu^[111], Li et al.^[112], Tan et al.^[113], and Fu and Li^[114] applied this method to DNS of compressible turbulence. Their results showed that the high-order kinetic scheme is superior to the second-order scheme with regard to capturing the fine structure of flow fields (such as in near-wall turbulence). However, the computational complexity is significant. Further research should be conducted to optimize the reconstruction quality and to improve the efficiency of this scheme.

Turbulence models are a key element in CFD and have a decisive influence on the accuracy of numerical simulations. A good turbulence model allows for significant savings on experimental costs. At a low Reynolds number, a nonlinear eddy viscosity model was established on the basis of the FRT model, which itself is based on the models of Fu and Guo^[115]. This Fu-Rung-Thiele (FRT) model does not contain wall parameters and is suitable for turbulence simulations and complex flow fields. In compressible turbulence, the fluid exhibits strong anisotropy. Based on a theoretical analysis of the dynamic pressure and fluctuating velocity fields in compressible turbulence and on DNS results, Zhang and Fu^[116] formulated a compressibility correction model for correlation of the pressure deformation rate in compressible and uniform shear turbulence.

LES has progressed significantly in recent years, and is likely to become the principal method for turbulence simulation in the future. Peking University has completed some pioneering work on the LES of turbulence. First, a constrained LES was proposed. Physical systems are subject to a number of fundamental conservation laws and symmetries. In turbulence, these are manifested as a set of physical constraints. The idea of constraints in this context can be traced back to Kraichnan and Chen's 1987 decimation theory^[117]. In addition, Meneveau^[118] enumerated the statistical relationships to be satisfied by the LES subgrid stress model. In the application of LES, these physical constraints are often destroyed. By introducing physical constraints, the Peking University team was able to solve a series of previously intractable problems in LES.

The dynamic subgrid scale stress model is currently considered to be a promising LES model. It takes advantage of the mathematical properties of the subgrid stress in the turbulent field. The model coefficient is determined by a variational principle, with no need for human intervention. After the kinetic model was proposed in the 1990s, it received great attention and has been widely used.

Chen's team extended the variational principle of the dynamical model to a constrained variational principle by introducing constraints, thereby solving the problem of imprecise simulation of the energy spectrum in a hybrid injection dynamics model^[119]. This work is considered to be an important development of the kinetic model^[120].

At present, the main challenge that LES faces is the simulation of near-wall turbulence at a high Reynolds number. To avoid high cost incurred in calculating near-wall turbulence in an LES simulation, the Reynolds-averaged Navier-Stokes (RANS)/LES method was proposed. This method combines the advantages of the RANS method with those of LES, while avoiding their respective defects. A hybrid RANS/LES method, called the detached eddy simulation (DES) method, has been successfully applied to aerodynamic simulation and is considered to be one of the most promising methods for turbulence simulation.

At present, the RANS/LES method is in its infancy, the most important remaining problem being the lack of small-scale fluctuations in the near-wall simulation, leading to the emergence of a nonphysical transition zone between the RANS simulation area and the LES simulation area and deleterious effects on simulation accuracy.

Chen et al.^[121] proposed a Reynolds stress constraint on the near-wall LES model. This method produces sufficient small-scale fluctuations in the near wall simulation, thereby dealing with the shortcomings of the traditional RANS/LES mixing method. They also proposed an advanced turbulence model that can be applied to aerodynamic simulations.

Constrained LES (CLES) has been applied to the simulation of a number of aerodynamic problems. The CFD team at Peking University used independent BXCFD software in 2012 to implement CLES for simulation of the flow around a large airliner and revealed a more complex flow structure than the DES. Not only does CLES provide higher precision in aerodynamic simulations than those achievable with DES, it also provides a fine vortex structure, which can be used as a source of aerodynamic noise to provide a more accurate simulation of the aircraft noise during take-off and landing. This study was a cover story in Science China^[122].

As previously mentioned, the CLES method solves the nonphysical interface problem that plagues RANS/LES, and, in addition, the simulation accuracy is clearly better than that of the DES method, which was developed by Boeing. At the 23rd World Congress of Mechanics that was held in August 2012, Chen et al.^[123] made an opening presentation entitled multi-scale fluid mechanics and modeling. The fluid mechanics and turbulence portions of this report describe the introduction of an energy dissipation constraint into the dynamic model and a Reynolds stress constraint into the high-Reynolds-number wall turbulence simulation along with the simulation results of the CLES method for a large aircraft. This work received significant attention.

In nature, the atmospheric boundary layer is mostly in a turbulent state. Li and Xie^[124] and Li^[125] studied canopy turbulent flow and developed a transient structure function (TSD) subgrid model. Cui et al.^[126] compared the effects of different numerical methods on predicting urban micro-environmental flow using the LES method. Their results showed that the Lagrangian dynamic model (LDM) is a cost-effective subgrid model for the simulation of complex turbulence. Thus, Cui et al.^[127] proposed a new eddy viscosity model for LES of anisotropic turbulence, which was applied to rotating turbulence and wall turbulence.

Many fluid flows of practical importance have unsteady characteristics on multiple scales. Among studies of complex flows, those around a cylinder are of particular importance. Hence, both Lu's group^[128-130] and Zhou's group^[131-140] investigated unsteady flow in the boundary layer and wake of a cylinder. Li's research group used PIV technology to study the two-dimensional flapping flag in a water tunnel, and obtained experimental phenomena such as the period and bifurcation. In addition, they obtained the flow field data around the flag surface^[141]. Zhou's group studied the flows around multiple cylinders and square columns, developing a complex flow theory^[142-144] of practical engineering significance. Xia et al.^[145] studied the dissipation function in turbulent plane Poiseuille flows and plane Couette flows subject to spanwise rotations. Research on multi-scale complex flows also includes the fluid dynamics of swimming and flying organisms^[146], unsteady vortex separation flows, vortex control, hemodynamics and physiological microcirculation, micro scale industrial flows, and heat transfer. Based on CFD analysis, Lu's team studied the propulsion performance of flapping wings and the law of vortex shedding related to biological movement^[147], and improved the understanding of the mechanism of fish swimming in terms of drag reduction and propulsion efficiency^[148].

5 Turbulence experiment and measurement techniques

Since the original Reynolds' experiment, significant contributions have been made to the fluid mechanics community. Because wall turbulence is easier to operate experimentally, it has become a standard geometric model for modern turbulence research. In recent decades, Lee and other researchers have conducted in-depth studies of incompressible boundary layer transitions^[149-153]. Furthermore, Lee improved upon the Hama^[154] experiment and observed the flow structure under a finer natural transition through the parallel setting of hydrogen bubble lines at different heights. For the first time, the soliton-like coherent structure (SCS) (formed at the initial stage of transition) was proposed^[155]. Kachanov observed the periodic "spike (peak)" structure from the quantitative test of hot wire, and confirmed that the "spike" refers to the velocity traces of the chain of ring vortices. Lee obtained the results of the first flow visualization of the chain of ring vortices, which was the result of interaction between the secondary closed vortex and the Λ-vortex (see Fig. 3). This interaction also produced the streamwise vortex. Further experiments showed that the secondary vortex is generated by the SCS. These experiments indicated that the SCS is ever-present and dominates the transition process^{[149, 156-157]}.

In addition, Lee proposed a new idea to aid understanding of the turbulent transition model and initiated a self-consistent dynamic model for wall turbulence generation. He devised a general framework for the mechanism of the onset of turbulence. The formation, development, and interaction of various coherent structures were explained in detail. Zhao et al.^[158] proved that a solitary-like wave is a common structure in the transitional boundary layer. In their study of wall turbulent flow structure, Lian^[159-160] proposed and explained the important characteristics of the flow structure of the long strip (long streak) and black spot (black spot) in wall turbulence. Guo et al.^[161] used the hydrogen bubble method to visualize the laminar flow turbulent transition in the plate boundary layer flow in the water tunnel. Thus, they discovered the inverse hairpin vortex. Finally, a new method of delayed transition^[162], called the spanwise discrete suction, was developed to produce a high amplitude stable streak.

Gong et al.^[163] conducted a series of work on the relationships between the turbulent energy series, kinetic process, and turbulence occurrence. Using the main dynamic process established in the boundary layer transition, the concept of turbulence cascading was further analyzed. Following this, a quantitative analysis of the measurement results of boundary layer transitions was conducted. From the energy spectrum, the relationship between the turbulent cascade and the dynamic process of the transition was clearly visible^[164]. Moreover, when the energy spectrum is used to analyze the energy level series, it must be clear that the frequency distribution caused by physical or non-physical factors is transmitted to either a low or high frequency. Otherwise, a true and false situation will occur and an error analysis result will emerge^[165]. Through detailed chaotic dynamic analysis of the transition boundary layer, the existence of the attractor in the transition process was confirmed and the relationship between chaos and transition in the boundary layer was established^[166-167].

The intensity of the vortex structure in a turbulent boundary layer is significant for the generation and development of hairpin vortices. Gao et al.^[168] found that there is a significant difference between the real eigenvalue direction and the vorticity vector direction of the identified vortex set in the turbulent boundary layer. Based on the PIV and DNS data, the intensity of three-dimensional swirling flow, which is based on a local velocity gradient tensor, was used to identify the vortex core location^[169-170].

Based on the PIV and hydrogen bubble flow display, Wang et al.^[171], Zhang et al.^[172], and He et al.^[173] used PIV and hydrogen bubble flow visualization to study the effects of roughness on the bypass transition induced by cylindrical wake, the evolution of vortex structure in the boundary layer, and the initial development of the boundary layer disturbance under the effects of the wake. Thus, a new transition control strategy was proposed based on these investigations. In free-stream-induced transition^[174], the two-dimensional cylindrical array-induced transition mode shows differences from the Klebanoff mode. Furthermore, detailed studies show that the evolution dynamics and maintenance mechanism of hairpins and hairpackets are consistent with those in turbulent boundary layers. Wake vortices play an important role in generating and destroying the secondary transverse vortices. Thereafter, the internal mechanism becomes dominant and leads to the establishment of the self-sustaining turbulent boundary layer. The aforementioned study found that the flow of the hairpin vortex determines the logarithm of the turbulent boundary layer^[175-177]. In addition, it was found that the downstream convection of the hairpin vortex determines the transport characteristics of the turbulent boundary layer inside and outside the logarithmic region. The proper orthogonal decomposition (POD) of PIV data indicates a low velocity streak associated with the hairpin pocket and the hairpin vortex in the transitional turbulent boundary layer, which is the main coherent structure near the wall^[178]. Finally, the evolution of the single low-speed streak in the laminar boundary layer^[179] and the gap between understanding the streak characteristics and the dynamic mechanism of the free flow induced bypass transition were also studied.

The reverse transition phenomenon in a flat plate boundary layer was studied by Zhang et al.^[180], while the changes in the turbulent energy generation term, dissipation term, and absorption term before and after heating were measured. In addition, the physical mechanism of the reversal transition was preliminarily explained.

In addition to wall turbulence, complex flow turbulence has also been the focus of significant research. Important results were obtained concerning circular turbulent jets based on flow visualization and PIV. It was discovered that the development of the flow vortex plays an important role in the wrapping and rolling process of the near field of round turbulent jets. Furthermore, through the study of dynamic excitation control (excitation frequency range, amplitude, and mode) and static control, the deformation degree of the turbulent quasi order vortices is considered to be controllable to a certain degree^[181]. In the study of particle turbulence, owing to the limitation of experimental technology, complete field research was mainly focused on simple cases. At present, the PIV experimental technique for simultaneous measurement of a two-phase flow is a possible breakthrough direction^[182]. Turbulent thermal convection is also a large domain of focus and mainly includes turbulent heat transfer, quasi sequence structure (plume and large scale circulation) in turbulent heat convection, statistics of the turbulent pulse momentum of an RB system, and the study of some non-traditional means (rotation, modulation, stratification, and multi-phase Rayleigh-Bénard convection)^[183]. In free surface geomantic boundary experiments, Wang et al.^[184] made a series of important advances. Moreover, Li's group made a breakthrough in the complex flow problem, which was caused by a falling disk. Furthermore, they fully explained the complex vortex structure during the falling process, both theoretically and experimentally^[185-186].

Due to increasing aircraft speeds, the transition of the hypersonic compressible boundary layer has become a bottleneck in the aerospace field.

For the study of hypersonic boundary layer transitions, Peking University established two hypersonic quiet wind tunnels. Their team developed temperature sensitive painting (TSP), near wall PIV, Rayleigh scattering imaging (RSI), and so on. The first and second modes of boundary layer transition and the correlation between temperatures were studied. The first flow visualization of the complete transition process of the hypersonic boundary layer was obtained (see Fig. 4), and the two mode wave structures in the hypersonic flow field were clearly captured by the PIV technique for the first time. Furthermore, it was found that the second mode wave generated violent aerodynamic heating at the early stage of transition and that this value may exceed that of the aerodynamic heating once in the developed turbulent state. This process is closely related to the high frequency compression-expansion process, which is caused by the second mode wave and constructs a new principle of aerodynamic heating^[90]. In the experimental domain, Jia et al.^[187-188] produced a successful near wall measurement for a moving wall by means of the PIV method in the air for the first time. Li and Zhang^[189] and Zhang and Li^[190] studied the effect of glow discharge as an artificial disturbance on the first mode and the second mode. Experiment and linear stability results showed that the second-mode instability waves are significantly stimulated by the artificial disturbances, and the boundary layer transition is effectively triggered. Finally, Zhu et al.^[191] developed a new PIV data processing method, which can be used to obtain high accuracy data in a near wall region, and found a strong dilatation process within the second-mode instability by this method^[192]. Testing techniques of nano-tracer planar laser scattering techniques were used to study compressible turbulence^[193]. Clear flow structures were obtained in typical flow fields such as wake flow, supersonic boundary layer transition, and turbulent shock wave boundary layer interaction.

Land desertification and dust storms have posed great threats to human environments. Thus, Zheng's team conducted a series of studies on wall turbulence with both air and sand at high Reynolds numbers. During wind-blown sand movement, sand grains are charged to form a wind-blown sand electric field. Moreover, Huang and Zheng^[194] and Zheng et al.^[195] studied the electrification of sand particles in the sand movement, both theoretically and experimentally. Their results showed that the charge of sand particles is mainly produced by collision and friction between grains of different sizes, and the electric field in the wind sand flow is mainly formed by moving charged sand particles. Zheng et al.^[196] proposed a new theoretical model based on the theory of electron/hole exchange to explain the mechanism of the effects of relative humidity on the contact charges of identical insulator particles. Some numerical simulation results showed that the greater the friction wind speed, the smaller the particle size, the greater the calculated concentration at the same height, and the greater the evolution time^[197]. Hence, wind speed and particle size have significant influence on desert formation^[198].

The definition of whirlpools in the flow field is not clear. However, the generation and evolution of vortices in whirlpools play an important role in the flow. It is possible to use a flow control method to clarify the dynamic properties of the whirlpool and to improve the properties of the flow field by inputting small amounts of energy. Hence, Chen et al.^[199] studied the response of disturbance to asymmetric vortices. This work has potential practical significance for the flow control of vortices.

6 Prospects

In the early twentieth century, China's turbulence research was outstanding in the world because of the great work of Peiyuan Zhou. Of course, the great contributions of C. C. Lin and other overseas Chinese were not covered in this article. However, the contributions of Chinese scientists from across the globe, particularly during the second half of the twentieth century, have been included.

In summary, the study of turbulence in China began with Peiyuan Zhou. After more than 70 years, several generations of effort have resulted in some remarkable achievements in the field. With the development of China's economy and the need for national strategic scientific projects, turbulence research in China has been greatly promoted. In 2017, a China national strategic turbulence research program, i.e., Physical Mechanisms of the Formation, Evolution, and Interaction of Turbulence Structures supported by the National Natural Science Foundation of China (No. 91752000), was proposed and initiated under the strong leadership of Professor Shiyi Chen. It can be expected that, with the continued deepening of turbulence research, Chinese scholars will make more and more significant contributions in the field. These contributions will certainly have a strong impact on the establishment of a China School of Turbulence.

Acknowledgements

The author would like to acknowledge the support of Cunbiao LEE for his useful help.