Vortex structure identification for supersonic flows

This work is devoted to the application of vortex structure identification methods for two counter-rotating wingtip vortices in a supersonic flow. The λ2 criterion and the method of maximum vorticity were used for detection vortex core lines and visualization of vortex sheet and rope. Computational model based on the URANS equations with SA turbulence model was applied. Numerical simulations were performed on the hybrid supercomputing system K-60 at the Keldysh Institute of Applied Mathematics RAS using the developed software package ARES for 3D turbulent flows modeling on high performance computing systems.


Introduction
The flight of any aircraft in the atmosphere is accompanied by the formation of vortex structures behind it, in particular, tip vortices at the edges of the wings and other details. The study of the behavior of wake vortex is an important task of modern aerodynamics, both from the point of view of flight safety, and aircraft visibility problem for military flight vehicles.
To date, the main studies have been devoted to tip vortices formed in subsonic and transonic conditions. In recent decades, there has been a steady increase in the number of both experimental and numerical researches of supersonic tip vortices. Their correct numerical modeling usually requires high-resolution grids, which inevitably leads to a large amount of simulated data. For this reason, the use of special methods to identify vortex structures for this class of problems is especially relevant. Extensive reviews have been made on this topic by various authors, including [1,2,3].
In this work we present the results of numerical simulation of two counter-rotating supersonic vortices interaction, the methods of maximum vorticity and the λ 2 criterion were used to identify its core lines and visualization of vortex sheets and ropes. Both of these methods have Galilean invariance and are suitable for analysis of unsteady flows. Computational model based on the unsteady Reynolds averaged Navier-Stokes equations (URANS) equations with Spalart-Allmaras (SA) turbulence model was applied [4]. The numerical simulations were performed on the hybrid supercomputer system K-60 [5] at the Keldysh Institute of Applied Mathematics RAS using the developed software package ARES for 3D turbulent flows modeling on high performance computing systems [6]. The Tecplot software [7] was used to visualize the simulation results.

Model statement
Identification methods were applied to the numerical simulation results of two counter-rotating vortices. The supersonic flow behind two wings -wingtip vortex generators was studied. The general scheme of model is represented in fig. 1. Numerical simulations were carried out in dimensionless statement, per unit length was takeñ L = 1 m. The computational domain dimensions were L x = 0.35, L y = 0.225 and L z = 0.2. The origin of coordinates was located on the common axis of the wings. The axis Ox coincided with direction of free-stream flow, Oz -with common wings axis, Oy was directed from the leeward to windward side of the wings.
The wings were coaxial straight with sharp leading and trailing edges, sharp tip chord and diamond-shaped base. Its were placed at 10 • angle to the free-stream flow and were fixed with a base to the side walls parallel to flow direction. The chord of each wings was d = 0.03, the halfspans were equal to l 1 = 0.075 and l 2 = 0.095 respectively. The distance between wingtip chords was l 3 = 0.03. Thus, the length of the computational domain was 10 wing chords downstream of the common axis of the wings. The Mach number of free-stream flow was M ∞ = 3 and Reynolds number equaled to ReL = 1 × 10 7 .

Computational model
The numerical simulations were performed on 224 computational cores of the hybrid supercomputing system K-60 [5] at the Keldysh Institute of Applied Mathematics RAS using the developed software package ARES for 3D turbulent flows modeling [6].
The core of code based on system of URANS equations with SA turbulence model [4] for compressible flows [8] with Edwards modification [9] for describing a supersonic flow of a perfect viscous compressible fluid. Briefly the system of equations will be presented in the next subsection. The explicit finite volume method based on Godunov-type scheme [10] with the third order reconstruction method (WENO3) was used. A more detailed description of the numerical algorithms is given in [4]. Simulations were carried out on an unstructured mesh with 25 774 200 hexagonal cells. The mesh was refined in the zones of vortex formation and its propagation throughout all simulation domain for better resolution of vortex structures.
Initial conditions for numerical simulations were set to free stream parameters. On the inlet boundary the free stream flow was set with given ratio of the turbulent and dynamical viscosities for SA model (in this work µ t /µ = 1.0). On the vortex generators and on the walls of their attachment the no-slip boundary condition was used. On the other two side surfaces the free stream flow was prescribed. On the outlet boundary supersonic outlet condition was set.

Govening equations
The system of URANS equations in cartesian coordinates (x 1 = x, x 2 = y, x 3 = z) can be written in the following dimensionless form: where q -the vector of conservative variables, f j and g j -the vectors of inviscid and viscous fluxes respectively: Here ρ is the density, u j are the components of the velocity vector u, τ i,j and δ i,j are the components of the viscous stress tensor and Kronecker tensor, E is the total energy of the averaged flow, p and T -thermodynamic pressure and temperature, which are calculated from the perfect gas equation of state: where γ and R are adiabatic index and gas constant for air. The components of the viscous stress tensor and the heat flux vector have the form The effective values of the coefficients of dynamic viscosity and thermal conductivity are determined as where c p is the specific heat capacity coefficient, Pr and Pr t are laminar and turbulent Prandtl numbers, µ 0 is the dynamic viscosity at the reference temperature T 0 , µ t denotes the added turbulent viscosity, which is determined by the modified version of SA turbulence model: A detailed description of the model variables in the equation (6) can be found in [4]. All of the above variables are dimensionless. Nondimensionalization procedure was carried out using free-stream parameters according to the following formulas: where (.) denotes any dimensional value, c -sound of speed. Thus, dimensionless parameters of free-stream flow were equal to ρ ∞ = P ∞ = T ∞ = 1.0 and u 2 ∞ = M 2 γ.

λ 2 criterion
The λ 2 method (or criterion) for vortex identification was proposed in [11]. This criterion is formulated on the search for a pressure minimum across the vortex from the strain-rate transport equation. According to this criterion the vortex flow region is determined based on the analysis of the eigenvalues of the symmetric matrix A = S 2 + Ω 2 , where S and Ω are the strain rate and vorticity tensors of flow: In this case the vortex region is that part of the domain where the second eigenvalue λ 2 (A) < 0 (λ 1 ≥ λ 2 ≥ λ 3 ). This method is quite widespread and is often used in data processing. It should be noted that application of the criterion may be difficult in the case of presence of several individual vortices in the domain. In this situation, additional topological method are required to specify the vortex regions.

Maximum vorticity method
The maximum vorticity method was proposed in [12]. It is based on one of the definitions of vortex flow and consists in detection the local maximum of the vorticity vector modulus |ω| = |∇ × u| in the plane perpendicular to the direction of this vector. This method allows to determine the exact axis of the longitudinal vortex in the case of sufficient resolution of the computational mesh. As well as the λ 2 criterion, the maximum vorticity method requires combined use with topologic approaches in complex cases.

Results
Application of vortex identification methods to above problem statement is shown in fig. 2. The green areas denote vortex flow region bounded by the iso-surfaces of the λ 2 = −5000, red lines are vortex axes by the maximum vorticity method. The iso-surfaces on the wings have been clipped with geometric conditions for better visualization. It can be noted that these methods are higly complementary, and provide a detailed flow pattern. Fig. 3 shows the λ 2 level lines in sections perpendicular to the direction of free-stream flow x = 0.1, 0.2, 0.3, λ 2 = −5000 is separately highlighted in black. Displacement of counterrotating vortices upward (to the leeward side), its divergence at downstream from the wing axis and expansion of the vortex zone are observed. This result correlates with the experimental data of other authors [13].
The special interest is a visualization of the vortex formation region -zone of tip vortex rope, a more detailed conception of which can be obtained by vorticity iso-surfaces. For example, fig. 4 shows the sheet folding from the tip edge of the wing into a rope.

Conclusion
In this work we present the results of numerical simulation of two counter-rotating supersonic vortices interaction, the methods of maximum vorticity and the λ 2 criterion were used to identify its core lines and visualization of vortex sheets and ropes. It can be noted these methods are highly complementary, and provide a detailed flow pattern. The numerical results have demonstrated the robustness of these methods for the case of supersonic wingtip vortices detection.  Figure 2. Superposition of the vortex axes (red lines) found by the maximum vorticity method and the vortex region obtained using the λ 2 criterion.