Direct Numerical Simulation Code Validation for Compressible Shear Flows Using Linear Stability Theory

1.Universidade de São Paulo – Instituto de Ciências Matemáticas e de Computação – Departamento de Matemática Aplicada e Estatística – São Carlos/SP – Brazil. 2.Departamento de Ciência e Tecnologia Aeroespacial – Instituto de Aeronáutica e Espaço – Divisão de Propulsão Aeronáutica – São José dos Campos/SP – Brazil. Correspondence author: Leandro Franco de Souza | Universidade de São Paulo Instituto de Ciências Matemáticas e de Computação – Departamento de Matemática Aplicada e Estatística | Av. do Trabalhador São-Carlense, 400 – Centro | Cx. Postal: 668, CEP: 13560-970 – São Carlos/SP – Brazil | E-mail: lefraso@icmc.usp.br received: Dec. 12, 2016 | accepted: Jun. 06, 2017 section editor: John Cater ABSTRACT: In order to simulate compressible shear fl ow stability and aeroacoustic problems, a numerical code must be able to capture how a basefl ow behaves when submitted to small disturbances. If the disturbances are amplifi ed, the fl ow is unstable. The linear stability theory (LST) provides a framework to obtain information about the growth rate in relation to the excitation frequency for a given basefl ow. A linear direct numerical simulation (DNS) should capture the same growth rate as the LST, providing a severe test for the code. In the present study, DNS simulations of a two-dimensional compressible mixing layer and of a two-dimensional compressible plane jet are performed. Disturbances are introduced at the domain infl ow and spatial growth rates obtained with a DNS code are compared with growth rates obtained from LST analyses, for each basefl ow, in order to verify and validate the DNS code. The good comparison between DNS simulations and LST results indicates that the code is able to simulate compressible fl ow problems and it is possible to use it to perform numerical simulation of instability and aeroacoustic problems.


INTRODUCTION
Th ere is an increasing concern about noise-related health problems in many engineering areas.Some noise sources are due to turbulence in a fl ow fi eld and aeroacoustics is the area responsible for the investigation of noise generation and propagation (aerodynamic noise).Th e problem of noise generated by fl ow has some specifi c characteristics, such as multiple frequencies, multiple scales and propagation to long distances with very low dissipation (Tam 2004).To deal with these characteristics, a numerical simulation must have a transient numerical scheme, very low numerical dissipation and dispersion to preserve the propagation of acoustic waves and non-refl ecting boundary conditions (Tam 2004).
Two main strategies are frequently used to simulate aeroacoustic problems.Th e fi rst one is based on acoustic analogies (Lighthill 1952;Lighthill 1954;Curle 1955;Ffowcs-Williams and Hawkings 1969), where noise sources are obtained from computational fl uid dynamics (CFD) and its propagation is calculated solving an inhomogeneous wave equation for pressure, subjected to these 02/13 sources.The second approach uses a high fidelity simulation capable of computing directly the acoustic waves together with the fluid flow, which is called direct aeroacoustic simulation.
In order to perform direct aeroacoustic simulation, a numerical code was developed within our research group (Souza et al. 2005).It uses parallel processing with domain decomposition to perform the calculations in a feasible time.High order temporal and spatial discretization schemes were adopted to minimize the dispersion and dissipation phenomena on the simulation of the flow field and the acoustic waves.A series of strategies, such as filtering and mesh stretching as well as characteristic boundary conditions were implemented in order to obtain a proper solution of the aeroacoustic problems.
However, before using a direct numerical simulation (DNS) code for aeroacoustic predictions, it is a good practice to observe the code behavior when submitted to unsteady, small disturbances.The linear stability theory (LST) provides a framework to obtain information about the amplification rate in relation to the excitation frequency and corresponding eigenfunctions.The results from LST can then be used for comparison with numerical simulations.
In the present study, results are presented for DNS of a two-dimensional compressible mixing layer and two-dimensional compressible plane jet.Disturbances were introduced at the domain inflow and the code verification was done by comparing the spatial amplification rate obtained by the DNS code with the amplification rates obtained from a LST analysis of each baseflow.The good comparison obtained indicated that the code satisfactorily simulates compressible flow problems and can be used to perform direct aeroacoustic simulations.

FORMULATION AND NUMERICAL SCHEME
The code was developed to solve the two-dimensional unsteady compressible Navier-Stokes equations together with the continuity and energy equations in dimensionless form.A Cartesian reference frame (x, y) was adopted, where the x axis is aligned with the streamwise direction and the y axis with the normal direction.The corresponding velocity components are u and v.In conservative form, the solution vector is given by: A perfect gas relation for pressure is assumed to close the formulation .In the conservative form, the unknowns are the density ρ, the momentum densities ρu, ρv and the total energy E, given by: where c v is the specific heat at constant volume and T is the temperature.All equations can be written in vector notation as: where the flux vectors F and G are given by: where p is the normalized pressure.The normal stress components are given by: and the shear stress component is: The heat fluxes, considering the Fourier law for heat conduction, are given by: The dimensionless quantities are the Reynolds number Re, the Prandtl number Pr and the Mach number Ma.The dimensionless fluid properties are the viscosity µ, the thermal conductivity ϑ and the heat capacity ratio κ.All fluid properties are normalized by reference velocity u ~∞ , temperature T ~∞, density ρ ~∞, and viscosity μ ~∞, where the superscript ~ represents dimensional values and the subscript ∞ corresponding to free-stream values.Time integration is performed using a 4 th order, 4-steps Runge-Kutta scheme.Spatial discretization in the streamwise x and normal y directions uses a 6 th order compact finite difference method and the generated tridiagonal equation systems are solved using the Thomas algorithm.
A non-uniform grid is used.An interest region with uniform grid is delimited.Then a grid stretching of approximately 1% is applied in both directions that, together with spatial low-pass filtering, create a damping zone.Disturbances become increasingly badly resolved and are dissipated as they propagate through the damping zone before they reach the outflow boundaries.
The code was parallelized through domain decomposition in both directions using distributed-memory parallelization.A message passing interface (MPI) library is used for inter-node communications.A grid transformation in the (x, y) plane is used to map the physical grid in an equidistant (ξ, η) computational grid.The metrics required for the derivative calculations are presented below.The first derivatives are given by: where the metrics are given by: LST analysis was used to obtain the frequency with spatial growth rate and the corresponding eigenfunctions are introduced as disturbances at the simulation inflow boundary.

RESULTS AND DISCUSSIONS
Compressible Two-dimensional mixing layer A compressible mixing layer consists of two streams with different velocities flowing parallel to one another, so that the resulting streamwise velocity profile assumes an "S" shape.Due to the inflectional point on the velocity profile, the flow is unstable and transition to turbulence may occur due to the generation and growth of vortical structures.For the present verification, the flow configuration has been closely matched to the case investigated by Babucke et al. (2008).
The Mach numbers of the upper and lower streams are Ma 1 = 0.50 and Ma 2 = 0.25, respectively.The Reynolds number is 500, based on the vorticity thickness at the inflow, which is also used to normalize length scales in x and y directions.As in Babucke et al. 2008, the inflow (x 0 = 30) is chosen in such a way that the vorticity thickness is: The inflow streamwise velocity, density and temperature profiles are shown in Fig. 1.These profiles are obtained by solving the steady compressible two-dimensional boundary layer equations.
The eigenfunctions used to introduce the disturbances at the DNS inflow boundary were obtained by linear stability analysis.The dimensionless fundamental frequency is ω = 0.62930.The corresponding eigenfunction amplitude and phase are shown in Fig. 2.This frequency and three subharmonics were used to disturb the flow.The maximum amplitude of the dimensionless fundamental frequency is 2 × 10 -3 , while for the other subharmonics the maximum amplitude is 1 × 10 -3 .All amplitudes are normalized by the corresponding maximum streamwise velocity amplitude.Furthermore, a phase shift of Δθ = -0.028was introduced for the first subharmonic, Δθ = -0.141for the second and Δθ = -0.391for the third, as used by Babucke et al. (2008).
Th e computational mesh, with 2500 × 850 points in the x and y directions, respectively, is presented in Fig. 3, showing the domain regions and boundary conditions.In the longitudinal direction, the mesh has constant increment ∆x = 0.157 up to x = 300, forming the region of interest.Further downstream, the mesh is stretched with a rate of approximately 1.0%, forming a damping zone.In the normal direction, the mesh is stretched from the center to the boundaries, with the lowest increment of ∆y = -0.15 and the highest increment of ∆y = 1.06.Tests with a coarse and fi ner mesh indicate that the adopted mesh is suitable for the simulations carried out here.
At the infl ow, the undisturbed boundary conditions were defi ned by a null v velocity; the u velocity, temperature and density profi les are presented in Fig. 1.Th e total energy E is calculated by Eq. 2. Th e disturbances were introduced at each time step of the Runge-Kutta method.Furthermore, characteristic boundary conditions are used at the infl ow to avoid outgoing wave refl ections back to the computational domain.Th e characteristic variable c, which corresponds to the outgoing waves at the infl ow boundary, can be obtained from the basefl ow and fl ow fl uctuations by:  Once c (j=1) is calculated, the fl ow fl uctuations at the infl ow boundary can be calculated as: and these values are added to the infl ow boundary variables.
To minimize refl ections at the freestream (upper and lower boundary) and at the outfl ow boundaries, caused by normal and oblique acoustic waves (Babucke et al. 2008), a combination of grid stretching and spatial low-pass fi ltering is applied in the damping zone, forcing the fl ow variables to a steady-state solution.Th us, the disturbances become increasingly badly resolved as they propagate through this region and by applying a spatial fi lter, the perturbations are substantially dissipated before they reach the boundaries.Th e computational domain was decomposed into 2 parts in the y direction and 7 parts in the x direction, adopting 14 processing elements to perform the calculations.
Th e chosen time step was Δt = 2π/ω•n steps , where n steps is the number of points per wavelength of the fundamental frequency.In the present case n steps = 752 was adopted, corresponding to a time step of Δt = 1.33 × 10 -2 .Seventy-six periods of the fundamental frequency were simulated and the last eight were chosen for post-processing analysis.By considering 76 periods, we ensure that the results are free of numerical transients.
Th e maximum amplitude of the normal velocity component v, for each streamwise x position, shows how the disturbances behave (Fig. 4).In the initial part of the domain, the amplitudes grow exponentially until non-linear eff ects are observed, causing a saturation of the maximum amplitude of v.

(19) (20)
The spatial growth rate α i is calculated by Eq. 21.It can be compared with results from LST, as shown in Fig. 5. Table 1 presents the difference between DNS and LST results, where there is an initial portion of the domain with large error due to the receptivity region.Further downstream, where the transient of all three modes has vanished, the differences decrease as shown in Table 1.Although α i is a very sensitive value, mainly for the subharmonics, the mean values of the DNS corresponded well to those predicted by LST.Thus, the good agreement between simulation and LST serves as verification of the implemented DNS computational code for small disturbances propagation.Besides the spatial amplification rate comparison, the spanwise vorticity Ω z was compared with the results presented by Babucke et al. (2008), as shown in Fig. 6.It is possible to see the formation of the first Kelvin-Helmholtz vortices which will pair with neighboring ones to form new structures downstream.A good qualitative agreement was observed. (21)

08/13
Compressible Two-dimensional plane JeT Plane jets are a typical shear flow present in several applications, such as in combustion and propulsion.There is also a fast growing interest in aerodynamic noise generated by jets in the aircraft industry.In this flow, a central core with high velocity mixes with a parallel stream with lower velocity.In the interface between both flow streams, there are high gradients of the flow properties, making the flow unstable to disturbances.
In the present investigation, a two-dimensional plane jet was studied.As in Reichert and Biringen (1997), the inflow velocity profile is:  U ∞ was U j /U ∞ = 1.67 for all cases.The jet is perfectly expanded and isothermal such that the temperature, density and pressure profiles are kept constant.The adopted Prandtl number was Pr = 0.71.Since the DNS results are compared to inviscid linear stability analysis, the DNS simulations are done with a high Reynolds number (Re = 10000), thus the inertia effects are much higher than the viscous effects.
A computational Cartesian grid with 1400 × 800 nodes in x and y directions was used in all simulations.In the longitudinal direction, the grid is uniform with ∆x = 0.25 up to x ≈ 250 and then stretched at a 1.0% rate.In the normal direction, there is also a uniform grid with ∆y = 0.106 up to y ≈ ±5.20 and then stretched at a 1.0% rate.The uniform spacing zone corresponds to the interest zone, while the stretched grid region corresponds to the damping zone.
At the inflow, the undisturbed boundary conditions are given by Eq. 22 for the u velocity component and null for the v velocity component.The non-dimensional density and temperature are equal to one, and the total energy E is calculated by Eq. 2. As in the case of the two-dimensional mixing layer, the disturbances were introduced at each step of the Runge-Kutta method and the necessary care was taken to avoid numerical contamination of the solution.
The chosen time step was the same as the one previously considered Δt = 2π/ω•n steps , where n steps is the number of points per wavelength of the fundamental frequency.The number of time steps used for the fundamental frequency was n steps = 752.In each case, 76 periods were simulated for the fundamental frequency, where the last eight were chosen for post-processing analysis.An inviscid LST analysis was performed for each case, considering both sinuous and varicose modes.The simulated cases, their correspondent spatial amplification rate and frequency are presented in Tables 2 and 3, for sinuous and varicose modes, respectively.The frequencies were chosen close to the maximum amplification rate according to LST.The corresponding jet, freestream and convective Mach numbers (Ma j , Ma ∞ and Ma c ) are also shown in these tables.The convective Mach number is given by:  (23)

10/13
The variation of the spatial amplification rate with Ma j is shown for the sinuous mode in Fig. 8 and for the varicose mode in Fig. 9. Increasing the Mach number, there is a decrease in the amplification rate value, and the flow becomes more stable.The varicose mode is more sensitive to compressibility effects and the growth rates decay faster with increasing Mach numbers.The decrease of the amplification rate with compressibility and consequent mixing of the flow is a known effect.This effect is of great importance in combustion flows and can help in noise control.Figures 10 and 11 show the spanwise vorticity for the sinuous and varicose modes, where the vortex generation pattern is anti-symmetrical for the sinuous mode and symmetrical for the varicose mode.For both modes, increasing the Mach number is followed by an increase in the undisturbed jet core length and a lower spreading of the vorticity thickness, since the flow is more stable as predicted by the LST analysis.
The maximum amplitude of the normal velocity v along the y direction, as a function of the streamwise x position (|v| max ), is presented in Fig. 12 for the sinuous mode, and in Fig. 13 for the varicose mode.An initial receptivity region can be noticed  After the linear growth region, there is a region where the non-linear effects starts to become significant, where the LST theory is no longer valid.Since DNS results have these three distinct regions, where the linear region may be contaminated by the overlap of the initial transient and the nonlinear regions, it was decided to compare the maximum amplitude of the normal velocity with the amplitude 12/13 predicted by the LST theory, instead of the amplification rates.A good agreement was observed between LST and DNS results and this can be used to consider the DNS code verified.

Figure 1 .
Figure 1.Infl ow profi les for u, p and T for the 2D mixing layer.

Figure 4 .
Figure 4. Maximum amplitude of v for each x position.Comparison with results from Babucke et al. (2008).

Figure 5 .
Figure 5. Spatial amplification rate in the streamwise direction.Comparison between LST and DNS for the 2D mixing layer.

Figure 7 .
Figure 7. u velocity profile for the 2D plane jet.

Figure 8 .
Figure 8. Spatial amplification rates for the 2D plane jet, sinuous mode.

Figure 9 .
Figure 9. Spatial amplification rates for the 2D plane jet, varicose mode.

Figure 10 .
Figure 10.Spanwise vorticity for the 2D plane jet, sinuous mode.Contour levels from -0.5 to 0.5.Solid lines for positive values and dashed lines for negative values.

Figure 11 .
Figure 11.Spanwise vorticity for the 2D plane jet, varicose mode.Contour levels from -0.5 to 0.5.Solid lines for positive values and dashed lines for negative values.

Figure 13 .
Figure 13.Maximum streamwise amplitude variation of the normal velocity component for the varicose mode.

Table 1 .
Spatial amplification rate comparison between DNS and LST for the 2D mixing layer.

Table 2 .
Mach numbers and spatial amplification rates of LST simulations.Sinuous mode.

Table 3 .
Mach numbers and spatial amplification rates of LST simulations.Varicose mode.