An Improved Synthetic Eddy Method for Generating Inlet Turbulent Boundary Layers

Abstract

An improved synthetic eddy method (SEM) is proposed in this paper for generating the boundary layer at the inlet of a computational domain via direct numerical simulation. The improved SEM modified the definition of the radius and the velocities of the eddies according to the distance of the eddies from the wall in the synthetic region. The regeneration location of the eddies is also redefined. The simulation results show that the improved SEM generates turbulent fluctuations that closely match the DNS results of the experiments. The skin friction coefficient of the improved SEM recovers much faster and has lower dimensionless velocity at the outer of the boundary layer than that of the traditional SEM.

1. Introduction

The rapid growth of computation ability has had a great impact on the field of computational fluid dynamics (CFD). The great increase in computational power has brought the possibility of turbulence simulations (LES or DNS) of many complex flow problems. A fully developed turbulent boundary layer is required before the region of interest in the simulation domain or the authenticity of the whole simulation will be affected. A good turbulent inflow needs to be able to faithfully reconstruct all known turbulence information of a specific problem without introducing pseudowaves, independent of the computational grid, geometric configuration and numerical format, and does not affect large-scale parallel computing. Therefore, establishing turbulent boundary conditions is a very challenging task.

The generation methods of turbulent inflow conditions mainly include library-based methods, recycling-rescaling-based methods and synthetic turbulence methods [1,2,3,4]. Library-based methods adopt the turbulence library from experiments or DNS results. The phase information of the fluctuations can be well restored by library-based methods due to the reality of the turbulence library from experiments or DNS results [5]. However, these methods rely on periodically reading the library during the iterations. Periodic fluctuations may be generated in the streamwise direction if it is not treated well, and it is usually not easy to obtain the libraries corresponding to the calculation conditions. Therefore, library-based methods are rarely used now.

Recycling-rescaling-based methods generate turbulent inflow conditions by readjusting the velocity field downstream and reintroducing the velocity field at the inflow. Different similarity laws are used to rescale the mean velocity and fluctuations in the inner and outer regions of the boundary layer [6]. This method has a good application in the calculation of straight flow fields [7,8,9], but it is difficult to calculate the flow field with complex configurations. In addition, this method involves the correlation between a plane inside the flow field and the inlet plane, which will carry out additional communication in parallel calculations.

The simplest idea of the synthetic turbulence methods is to superimpose the desired average velocity distribution and random white noise on the inlet. However, due to the lack of any spatial or temporal coherence, the fluctuations generated by this simple method usually disappear within a short distance downstream [10]. Therefore, spatial and temporal coherence is added in practical synthetic turbulence generators to restrain fluctuations. The most common applications are the digital filter method [11] and synthetic eddy method [10]. The digital filter method establishes a digital filter of the Gaussian two-point correlation function. The filtered random signal has the required spatial length scale and can be further transformed to meet the specified Reynolds stress tensor. After continuous improvement, this method has been applied to many academic problems. Its obvious disadvantage is that the filtering operation is mainly suitable for finite difference separation and discretization using uniformly spaced Cartesian structured grids, which will be difficult in the application of complex configurations such as circular inlets.

In the synthetic vortex method, the superposition of coherent structures is simulated by introducing artificial eddies (the word “eddy” refers specifically to artificial eddies in the paper) into the inlet plane. Each eddy is a vorticity source with a three-dimensional structure represented by a specific shape function describing its spatial and temporal characteristics [12]. The velocity fluctuations at a given grid point on the inlet plane are calculated by superimposing the effects of all eddies. The synthetic eddy method has the advantages of the synthetic turbulence method; that is, it does not depend on the database, and vortex points do not depend on a specific grid, so it is easy to use in complex configurations and has no additional cost for parallel computing. Therefore, the synthetic eddy method is currently the mainstream generation method of turbulent inflow conditions.

Jarrin et al. [12,13] proposed the original SEM based on the principle of the vortex method. The core idea of the SEM is to generate turbulence through a series of eddies that satisfy the nondivergence condition. The spatial and temporal correlations of the turbulent structure are restored by the superposition of random pulsations and the flow directions of different eddies. SEM has been widely utilized and continuously improved since it was proposed. Pamiès et al. [14] used information in the available literature about the coherent structure in the turbulent boundary layer to design a shape function to simulate the vortex strip in the near-wall region and a hairpin vortex in the logarithmic layer and wake region. Poletto et al. [15] developed a nondivergent SEM by modifying the shape function for the action range of eddies. Tabor et al. [16] compared the recycling/rescaling method and SEM, analyzed their advantages and disadvantages, and reported that SEM is more flexible. Mankdadi [17] implemented the SEM and numerical filtering method based on the LES technique and discovered that the SEM can better reproduce the experimentally obtained turbulent stress. Lee et al. [18] used the turbulent kinetic energy at the corresponding position of each eddy to replace the Reynolds stress; accordingly, the turbulent kinetic energy and velocity time-averaged profile can be obtained by the Reynolds-averaged Navier–Stokes (RANS) equations, thereby reducing the difficulty in obtaining the Reynolds stress.

Several practical problems are found in our application of SEM:

• The traditional SEM seems to overestimate the eddy radii in the near wall region. The radii reach the maximum $0.41\delta$ ($\delta$ is the boundary layer thickness) at the approximate position of $0.2\delta$. The grids near the wall will be affected by too many eddies, which cause much larger vortex structures to be generated than small vortex structures. In the compressible LES study of Mankbadi et al. [17], it is found that SEM could succeed in replicating the Reynolds stress tensor when the radii of the eddies are taken as $0.1\delta$ after a preliminary boundary layer study where the radii were varied, while the moving speed of the eddies is $0.5U_{\infty }$. This also indicated overestimates of the eddy radii in the near-wall region. Detailed comparisons of the radii can be seen in Section 3.
• The velocity of the actual vortex structures should be close to the average velocity of the turbulent boundary layer. The velocity of vortex structures in the near wall region is much lower than that in the outer layer of the boundary layer. Therefore, it is unreasonable to set all the eddy moving velocities as the mainstream velocity $U_{\infty }$. This may be the reason Mankbadi et al. [17] chose $0.5U_{\infty }$ eddy moving velocities.
• The vortex points are evenly distributed in the eddy region at the first step, and the vortex points leaving the outlet plane of the eddy region are regenerated on the inlet plane of the vortex region, which means $x_k^n=-{\sigma }_{max}$. This will cause a problem. After several or more steps, the streamwise distribution of the vortex will be sorted to several points at a time step away, thus losing the randomness of the flow direction, which is shown in Figure 1.

Considering these shortcomings, this paper studies the parameters of a supersonic turbulent boundary layer with Ma = 2.7 to propose a set of easy-to-use SEM applications. Part 2 briefly explains the basic principle of the SEM, Part 3 presents the specific details of the improved method, and finally, in Part 4, the method is applied to generate the turbulent boundary layer at the entrance of a flat plate to demonstrate the effectiveness of the proposed method.

2. Numerical Methods

2.1. NS Code

The simulations are run within an in-house GPU-based 3D unsteady compressible NS equation solver; more details about the code can be found in the literature [19].

A low-dissipation fourth-order hybrid of the central/WENO scheme is used for inviscid fluxes. An energy-preserving central difference scheme is used for the central difference scheme. In terms of the fourth-order energy-conservation scheme, the semidiscretizations of the numerical approximations are:

$$f_{i+1/2}^{EC}=\frac{1}{48}\left[\begin{array}{l}-({\rho }_{i+2}+{\rho }_i)(u_{i+2}+u_i)({\phi }_{i+2}+{\phi }_i)\\ +7({\rho }_{i+1}+{\rho }_i)(u_{i+1}+u_i)({\phi }_{i+1}+{\phi }_i)\\ +7({\rho }_i+{\rho }_i)(u_i+u_i)({\phi }_i+{\phi }_i)\\ -({\rho }_{i-1}+{\rho }_i)(u_{i-1}+u_i)({\phi }_{i-1}+{\phi }_i)\end{array}\right]$$

(1)

where $\phi$ stands for a generic transported scalar property, being unity for the continuity equation or $u_k$ for the momentum equations or $H=\frac{\gamma }{\gamma -1}\frac{p}{\rho }+\frac{u^2}{2}$ for the total energy equation.

Central-upwind fourth-order WENOCU4 scheme for the WENO scheme. It is defined as:

$${\widehat{f}}_{i+1/2}^{WENO}={\displaystyle \sum _{k=0}^{r-1}{\omega }_k{\widehat{f}}_{k,i+1/2},k=0,1,2}$$

(2)

where:

$${\omega }_k=\frac{{\alpha }_k}{{\displaystyle {\sum }_{k=0}^2{\alpha }_k}},{\alpha }_k=d_k\left(C+\frac{{\tau }_3}{{\beta }_k+\epsilon }\right),k=0,1,2$$

(3)

$$d_0=\frac{1}{6},d_1=\frac{2}{3},d_2=\frac{1}{6}$$

(4)

where $\epsilon$ is a small positive value to prevent the denominator from becoming zero and ${\beta }_k$ is the smoothness indicator, which can be given by:

$$\left\{\begin{array}{l}{\beta }_0={\left(-f_{i-1}+f_i\right)}^2\\ {\beta }_1={\left(-f_i+f_{i+1}\right)}^2\\ {\beta }_2={\left(-f_{i+1}+f_{i+2}\right)}^2\\ {\beta }_3=\frac{1}{36}{\left(-2f_{i-1}-3f_i+6f_{i+1}-f_{i-2}\right)}^2\\ \hspace{1em}\hspace{1em}+\frac{13}{12}{\left(f_{i-1}-2f_i+f_{i+1}\right)}^2\\ \hspace{1em}\hspace{1em}+\frac{781}{720}{\left(-f_{i-1}+3f_i-3f_{i+1}+f_{i+2}\right)}^2\end{array}\right.$$

(5)

where ${\widehat{f}}_{k,i+1/2}$ represents the substencil fluxes obtained from the node value $f_i$ as:

$$\left\{\begin{array}{l}{\widehat{f}}_{0,i+1/2}=-\frac{1}{2}{f}_{i-1}+\frac{3}{2}{f}_i\\ {\widehat{f}}_{1,i+1/2}=\frac{1}{2}{f}_i+\frac{1}{2}{f}_{i+1}\\ {\widehat{f}}_{2,i+1/2}=\frac{3}{2}{f}_{i-1}-\frac{1}{2}{f}_{i+2}\\ {\widehat{f}}_{3,i+1/2}=-\frac{1}{12}{f}_{i-1}+\frac{7}{12}{f}_i+\frac{7}{12}{f}_{i+1}-\frac{1}{12}{f}_{i+2}\end{array}\right.$$

(6)

${\tau }_3$ is the reference smoothness indicator as:

$${\tau }_3={\beta }_3-\left(\frac{{\beta }_0+{\beta }_1}{2}\right)$$

(7)

The hybrid algorithm is that the shock-capturing WENOCU4 scheme is only activated near critical cells, while the fourth-order energy conservation scheme is employed for the other flow field with smooth scales. Then, the hybridization is devised as:

$${\widehat{f}}_{i+1/2}=\left(1-\psi \right){\widehat{f}}_{i+1/2}^{EC}+\psi {\widehat{f}}_{i+1/2}^{WENO}$$

(8)

where the weighting function $\psi$ is defined as:

$$\psi =\left\{\begin{array}{l}1, \mathrm{if} \theta >\xi \\ 0, \mathrm{otherwise}\end{array}\right.$$

(9)

where $\xi$ controls the sensitivity to identify the critical regions, which depends on the stability of the program. Generally, we have $\xi$ = 0.98 in our application to simulate the supersonic boundary layer. The shock sensor $\theta$ can be calculated by the nonlinear weight of the substencils. The formula is defined in Equation (10).

$$\theta =\frac{{{\displaystyle {\sum }_{l=0}^3\left|\left({\omega }_l/d_l\right)-1\right|}}^2}{{\left|1/min{d}_m-1\right|}^2+3},l,m=0,1,2$$

(10)

The diffusive term adopts a fourth-order central difference scheme. The time integration scheme uses the third-order explicit Runge–Kutta method.

2.2. Theory of the SEM

Jarrin et al. [12,13] systematically described the implementation of the SEM. The basic theory of the traditional SEM is briefly described here.

First, a synthetic eddy region (ER), shown as Figure 2, a box region symmetrical about the entrance plane, is established near the inflow boundary of the computational domain, as shown in Figure 2. Enough virtual eddies are scattered in this region. The bottom wall of this synthetic region is a nonslip boundary, and the height of the top boundary is greater than the thickness of the turbulent boundary layer. In addition, the thickness of the synthetic region is the double maximum of the radium of the eddies.

The eddies within the ER do not possess a physical vortex structure and are used to generate a random process to simulate the randomness of velocity fluctuations. Furthermore, these eddies have a fixed radius, which is usually comparable to the average grid scale or slightly smaller than the maximum grid scale within the boundary layer. The eddies affect the magnitude of the pulsations on the grid points located within their radii.

The time-averaged flow velocity at each grid point is first determined according to its distance from the wall, and its final velocity is the sum of this time-averaged velocity and the velocity fluctuations. Specifically, the effects of the eddies that can disturb the grid point are superimposed, and this superposition yields a value that quantifies the total disturbance effect, which is eventually fed back to the velocity fluctuations with the help of the known Reynolds stress corresponding to the velocity fluctuations at the grid point.

In general, the eddies are evenly distributed in the synthetic eddy region so that each grid point can be disturbed by a certain number of eddies as long as the number of eddies is guaranteed to be sufficiently large. Figure 3 presents a schematic representation of the movement of a certain eddy, which moves forward with the incoming flow at a certain speed that is often set as a constant in the numerical simulation.

The movement of the eddies is described as follows:

(1)

The eddies are evenly initialized in the synthetic vortex region with a certain radius and carry random information that can influence the generation of turbulent fluctuations at the first iteration. All information about the eddy $k$ can be expressed as $\left(x_k,y_k,z_k,{\epsilon }_k,{\sigma }_k,U_k\right)$, where $\left(x_k,y_k,z_k\right)$ is the initialized location in the synthetic vortex region, $\epsilon =\left(random{(\pm 1)}_x,random{(\pm 1)}_y,random{(\pm 1)}_z\right)$ determines the directions of the velocity fluctuations, and ${\sigma }_k$ and $U_k$ are the radius and the moving velocity of the eddy, respectively. All the eddy information remains unchanged until the eddy regenerates, except $x_k$.

The profile of the radii of the eddies proposed by Jarrin et al. [12] is:

$${\sigma }_1=max\left(min\left[L_t,0.41\delta \right],{\Delta }_{max}\right)$$

(11)

where the integral turbulence length scale is $L_t=\frac{k^{\frac{3}{2}}}{\epsilon }$, $k$ is the turbulence kinetic energy, $\epsilon$ is the dissipation rate, and ${\Delta }_{max}$ is the max grid of the inlet.

The eddies constantly move downstream with the numerical iteration process with the velocity $U_k=U_{\infty }$. Therefore, at each iteration, the new position of eddy $k$ is given by $x_k\left(t+\Delta t\right)=x_k\left(t\right)+U_{\infty }\Delta t$.

(2)

The eddies in motion disturb the grid points within their radii, affecting the generation of velocity fluctuations. The velocity of a grid point on the inlet boundary can be expressed as:

$$u_i={\overline{U}}_i+u_i{}^{\prime }$$

(12)

where $i$ represents the three spatial directions, ${\overline{U}}_i$ is the time-averaged velocity, and $u_i{}^{\prime }$ is the pulsation velocity.

The total number of small vortices is set to $N=max\left(V/{\sigma }^3\right)$, where $V$ is the volume of the eddy region. Then, the pulsation velocity can be calculated by the following equation:

$$u^{\prime }=\frac{1}{\sqrt{N}}{\displaystyle \sum _{k=1}^N\left[a_{ij}{\epsilon }_{k}f(x,y,z,x_k,y_k,z_k,{\sigma }_k)\right]}$$

(13)

where $k$ represents the sequence of the eddies, $a_{ij}$ signifies the Cholesky decomposition values of the turbulent Reynolds stress tensor at the grid points, and the Reynolds stress tensor mentioned herein refers to the matrix of the second-order pulsation velocity:

$$R_{ij}=\left[\begin{array}{ccc}{u}^{\prime }{u}^{\prime }& u^{\prime }{v}^{\prime }& u^{\prime }{w}^{\prime }\\ u^{\prime }{v}^{\prime }& v^{\prime }{v}^{\prime }& v^{\prime }{w}^{\prime }\\ u^{\prime }{w}^{\prime }& v^{\prime }{w}^{\prime }& w^{\prime }{w}^{\prime }\end{array}\right]$$

(14)

$f\left(x,y,z,x_k,y_k,z_k,{\sigma }_k\right)$ is a decay function of eddy $\left(x_k,y_k,z_k\right)$ at grid point $\left(x,y,z\right)$. This function has the following characteristics: if the grid point is within the radius of the eddy, $f>0$; otherwise, $f=0$. $f$ is defined as follows:

$$f\left(x,y,z,x_k,y_k,z_k,{\sigma }_k\right)=\frac{\sqrt{V}}{{\sigma }_k^{1/3}}\psi \left(\frac{x-x_k}{{\sigma }_k}\right)\psi \left(\frac{y-y_k}{{\sigma }_k}\right)\psi \left(\frac{z-z_k}{{\sigma }_k}\right)$$

(15)

where $V_e$ is the volume of the synthetic eddy region and $\psi$ is defined as:

$$\psi \left(\theta \right)=\left\{\begin{array}{l}\sqrt{\frac{3}{2}}\left(1-\left|\theta \right|\right), \left|\theta \right|<1\\ 0,\left|\theta \right|\ge 1\end{array}\right.$$

(16)

where $\theta =\frac{x-x_k}{{\sigma }_k},\frac{y-y_k}{{\sigma }_k},\frac{z-z_k}{{\sigma }_k}$, respectively.

(3)

Once an eddy escapes from the outlet plane of the synthetic eddy region, it is regenerated randomly in the inflow plane synthetic vortex region to keep the total number of eddies constant, which can be given as:

$$\left[\begin{array}{l}{x}_k\left(t\right)\\ y_k\left(t\right)\\ z_k\left(t\right)\end{array}\right]\underset{\Delta t}{\overset{\mathrm{if} x_k\left(t+\Delta t\right)>{\sigma }_{max}}{\to }}\left[\begin{array}{l}{x}_k\left(t+\Delta t\right)=-{\sigma }_{max}\\ y_k\left(t+\Delta t\right)=random(y_k)\\ z_k\left(t+\Delta t\right)=random(z_k)\end{array}\right]$$

(17)

As long as the effects of eddies are similar enough to the natural random process, the fluctuations generated by the eddies are capable of reconstructing the turbulence characteristics of the boundary layer of the flow.

3. Improved SEM

Due to the shortcomings, we modified the original SEM.

• The radius profile of eddies is shown as ${\sigma }_3$ in Figure 4. It is defined as:

  $${\sigma }_3=\left\{\begin{array}{l}max\left(0.41y_w,{\Delta }_{max}\right),y_w<\delta \\ max\left(-0.82y_w+0.41\delta ,{\Delta }_{max}\right),y_w\ge \delta \end{array}\right.$$

  (18)

  where $y_w$ is the distance of the eddies to the wall. The linear distribution along the wall normal direction is adopted. We also tried the velocity distribution of a velocity profile shape, which was slightly different from the linear distribution. Therefore, we use the linear distribution for ease of use. The profile ${\sigma }_2=0.1\delta$ is adopted by Mankbadi et al. [17]. Compared with the profile of the radii of the eddies proposed by Jarrin et al. [12], the profile proposed in this paper reduces the eddy radii near the wall, while it is closer to that proposed by Mankbadi et al. [17].

With the change in radius, the total number of eddies can be calculated as follows: First, an average radius is defined:

$${\sigma }_e=\left(\underset{y=0}{\overset{\delta }{\int }}\sigma \left(y\right)dy\right)/\delta \approx 0.2\delta$$

(19)

Then, the total number of eddies is approximately determined as:

$$N=\frac{V_e}{{\sigma }_e^3}$$

(20)

2.

The moving velocity of the eddies adopts the average velocity of the turbulent boundary layer. However, to avoid the parallel efficiency loss caused by the interpolation of the velocity of the eddies regenerated in each iteration, an approximate expression $U_x={(y_w/\delta )}^{1/7}{U}_{\infty }$ [20] of the velocity profile is adopted. Therefore, the position of each eddy is updated after a time step $\Delta t$

$$x_k\left(t+\Delta t\right)=x_k\left(t\right)+U_k\left(y_k\right)\Delta t$$

(21)

$y_k$ and $z_k$ remain constant before the eddy vanishes.

3.

The new positions of the eddies are set as follows when the streamwise positions of the eddies are greater than ${\sigma }_{max}$

$$\left[\begin{array}{l}{x}_k\left(t\right)\\ y_k\left(t\right)\\ z_k\left(t\right)\end{array}\right]\underset{\Delta t}{\overset{\mathrm{if} x_k\left(t+\Delta t\right)>{\sigma }_{max}}{\to }}\left[\begin{array}{l}{x}_k\left(t+\Delta t\right)=x_k\left(t\right)-2{\sigma }_{max}\\ y_k\left(t+\Delta t\right)=\mathrm{random}(y_k)\\ z_k\left(t+\Delta t\right)=random(z_k)\end{array}\right]$$

(22)

For the temperature and density, the improved SEM correlates them with velocity fluctuations by a strong Reynolds analogy:

$$\frac{T^{\prime }}{\overline{T}}=-\left(\gamma -1\right){\overline{Ma}}^2\frac{u^{\prime }}{\overline{U}}$$

(23)

The pulsation density is then obtained with the ideal gas equation of state ${\rho }^{\prime }/\overline{\rho }=T^{\prime }/\overline{T}$, where $\overline{T}$, $\overline{Ma}$ sand $\overline{\rho }$ are the local time-averaged temperature, Mach number and density, respectively.

4. Application of the Improved SEM to a Supersonic Plate Flow

This paper obtains the results of the DNS [1] of the turbulent boundary layer of a plate with $M{a}_{\infty }=2.52$, from which the Reynolds stress tensor components are extracted, as shown in Figure 5. Since the Reynolds shear stresses $u^{\prime }{v}^{\prime }$, $u^{\prime }{w}^{\prime }$ and $v^{\prime }{w}^{\prime }$ are small compared to the main diagonal elements of the tensor, to conserve memory and reduce the computational overhead, it is assumed that $u^{\prime }{w}^{\prime }$ and $v^{\prime }{w}^{\prime }$ have the same distribution as $u^{\prime }{v}^{\prime }$.

The computational domain is shown in Figure 6. A fine mesh of 473 × 360 × 270 grids is used for the current simulation as shown in

Table 1. Computational grids for the mesh of the channel domain.

$N_x$	$N_y$	$N_z$	Total (Million)	${\Delta }_x^{+}$	${\Delta }_y^{+}$	${\Delta }_z^{+}$
473	360	270	45.98	8.2–13.2	1–12	4.1

. To minimize artificial reflections from the boundaries, large buffer regions (2.7 times the boundary layer thickness, as shown in Figure 6) are introduced on both the side boundaries and the upper boundary. The streamwise and spanwise grids are evenly distributed. In the wall normal direction, the mesh has 21 evenly distributed grid points near the wall as the inner mesh with a resolution of $y^{+}=1$, 191 grid points stretch from the inner mesh to $0.6\delta$, 121 grid points are evenly distributed from $0.6\delta$ to $1.5\delta$, and the rest stretch to the up boundary. The outlet of the computational domain is set as the supersonic outlet, and the thickness of the boundary layer is set to δ = 1.8 mm.

Two simulations are carried out in this paper based on the mesh. The program can run stably when the time step is 4 × 0⁻¹⁰ s, while the corresponding Courant number is approximately 0.5.

The average statistical characters are counted after 300,000 time steps (two overall flow times), and another 150,000 time steps are counted for the statistical characters (one overall flow time). Every step is counted. Both cases are treated with the same methods.

To examine the generation of pulsations at the inlet boundary for the improved SEM, the contours of the flow pulsation velocity at the inlet boundary monitored with the time step are plotted in Figure 7. To clearly distinguish the pulsation velocity direction, red represents a positive pulsation, and blue represents a negative pulsation. The area marked by the black dashed circle in each panel reflects the evolution of a positive pulsation velocity occurring in this series of consecutive moments. The direction and amplitude of the inlet velocity fluctuations are physically realistic based on the continuous forward motion of the eddy. On the other hand, the overall distribution of positive and negative velocity fluctuations in the boundary layer is statistically homogeneous owing to the random dispersion of the eddies, although the characteristic scale of these pulsations varies significantly at different heights.

Figure 8 shows the contours of the transient vortex distribution colored by the local density. The vorticity is concentrated mainly near the wall. The distribution of these vortex contours further reveals the difference between the bottom layer and the outermost layer of the turbulent boundary layer: in the region close to the wall, where the nonslip viscosity of the wall plays a major role, there are many small, closed vortex contours, while the upper region consists mostly of a long strip-like structure split by the strong shear effect of the main flow and the boundary layer. Moreover, the density distribution exhibits a phenomenon where the pulsation density is strongest at the entrance and then gradually weakens downstream. This may be because the pulsation density distribution given by the SEM is not physically realistic near the entrance but is constantly corrected downstream by the flow itself before finally returning to the proper pulsating state.

Figure 9 visualizes the transient dimensionless flow velocity contours and the second invariant of the velocity gradient tensor along the spanwise center section. The second invariant of the velocity gradient tensor is defined as:

$$Q=-\frac{1}{2}\frac{\partial u_i}{\partial x_j}\frac{\partial u_j}{\partial x_i}$$

(24)

which theoretically characterizes the relative equilibrium between the vortex strength and the local shear stress rate. $Q$ provides an intuitive indication of the vortex structure in the flow field. The development of the vortices along the stream is not uniform; for instance, the vortices in the outermost boundary layer go through a process of disappearing and then gradually regenerating. In particular, the velocity distributions of the flow in the boundary layer and the main flow area can be observed from Figure 9. Due to the pulsation of the velocity in the boundary layer, the main flow fluctuates along the upper edge of the boundary layer, and this fluctuation, which originates in the shear region near the inlet between the main flow area and boundary layer, propagates downstream in the main flow.

Figure 10 shows the isosurfaces of the 3D structure of $Q$ colored by the instantaneous streamwise velocity. In both cases, the turbulence structures introduced at the inflow boundary are adjusted to some degree, causing the turbulence to initially decay. This decay is followed by an increase in turbulence as the turbulent boundary layer tends toward equilibrium. The isosurfaces demonstrate that the adjustment distance generated by the traditional SEM is longer than that generated by the improved SEM.

The mean velocity distributions at $x=12\delta$ are shown in Figure 11. The profiles are transformed with

$$u^{+}=\frac{u}{u_{\tau }}$$

(25)

and scaled with $\sqrt{\rho /{\rho }_w}$:

$$u_{vd}^{+}={\displaystyle {\int }_0^{u^{+}}\sqrt{\frac{\rho }{{\rho }_w}}}d{u}^{+}$$

(26)

The density-weighted transformation method was developed by van Driest [21] for compressible fluids. The wall-normal coordinate is transformed according to:

$$y^{+}=y\frac{u_{\tau }}{{\upsilon }_w}$$

(27)

Figure 11 plots the profiles obtained with the traditional SEM and improved SEM and from the DNS by Morgan et al. [1]. The profile acquired with the traditional SEM is almost the same as that with the improved SEM but deviates from the linear sublayer profile and is slightly higher than the log layer profile. A larger velocity is obtained with the traditional SEM along the wall at the edge of the boundary layer due to the underestimation of $u_{\tau }$ with the traditional SEM. The current simulation results are very consistent with those obtained by the DNS conducted by Morgan et al. [1].

The skin friction coefficient is defined as:

$$C_f={\mu }_w\frac{\mathrm{d}{u}_1}{\mathrm{d}{y}_1}/(0.5{\rho }_{\infty }{U}_{\infty }^2)$$

(28)

In the simulation the skin friction coefficient was found to be $C_f\approx$ 0.00273 with the improved SEM and the corresponding Reynolds number based on the momentum thickness ${\mathrm{Re}}_{\theta }$ = 1456. The simulation compares favorably with the experimental results compiled by Coles (1954) and the skin friction correlation given in Bardina, Huang and Coakley (1997) based on the Van Driest II skin friction transformation in Figure 12. The streamwise skin friction coefficients are shown Figure 13. After experiencing a very similar trough, the $C_f$ of the improved SEM rose more rapidly than the traditional SEM. Finally, the restored $C_f$ of the improved SEM is slightly larger, which may be the reason for the lower velocity profile of the outer boundary region.

Further quantitative comparisons can be achieved by calculating the root mean square values of the pulsation velocity at the downstream exit boundary.

The root mean square values of the time-averaged pulsation velocity at $x=12\delta$ calculated by the traditional SEM (total number of small eddies N = 4000) and the improved SEM are compared with the values extracted from the DNS results in Figure 14. For the same number of small eddies, the traditional SEM underestimates the longitudinal pulsation velocity within the boundary layer, especially in the interior of the boundary layer ($y/\delta <0.2$). The improved SEM, on the other hand, restores the DNS results relatively well.

5. Conclusions

In this paper, direct numerical simulations of a supersonic turbulent boundary layer are carried out for a freestream Ma of 2.52. An improved SEM method for generating turbulent boundary layers is proposed. Three corrections have been made on the basis of the traditional SEM:

(1)

The velocities of the eddies are redefined according to the distance of the eddies from the wall in the synthetic eddy region, as the traditional SEM overestimates the eddy radii in the near wall region.

(2)

The moving velocity of the eddies adopts the average velocity of the turbulent boundary layer. An approximate expression, $U_x={(y_w/\delta )}^{1/7}{U}_{\infty }$ [20], of the velocity profile is adopted to avoid velocity interpolation.

(3)

Regenerated streamwise coordinates of the eddies in the synthetic eddy region are also modified to avoid random reduction.

The improved SEM is more effective than the traditional SEM in the simulation of the supersonic boundary layer flow along the wall of a flat plate. The skin friction coefficient of the improved SEM recovers much faster and has lower dimensionless velocity at the outer of the boundary layer than that of the traditional SEM. The improved SEM needs a shorter distance to recover the fully developed turbulence. The final profiles of the mean velocity and the root mean square velocity fluctuations match those provided by Guarini et al. [1] well.