Numerical simulation of the flow into a circular pipe section

. Computational Fluid dynamics (CFD) is the science that evolves rapidly in numerical solving of fluid motion equations to produce quantitative results and analyses of phenomena encountered in the fluid flow. When properly used, CFD is often ideal for parameterization studies or physical significance investigations of flow that would otherwise be impossible to replicate through theoretical or experimental tests. The aim of this paper is the study of the turbulent airflow and how the vortices formed in turbulent airflow are influenced by the evolution of the hydraulic characteristics of the fluid flow. Unsteady numerical simulation were performed using Reynolds Average Navier-Stokes (RANS) turbulence model adapted to conventional flow into a pipe with variable section which was implemented in the ANSYS FLUENT expert software.


Introduction
This study aims to analyze the vortex structures formed in the aerodynamic tunnel EOL located at the National Aerospace Research Institute (INCAS).For i vortices, a virtual model was created to reproduce the 3D geometry of the tunnel.The airflow, in the virtual model, was described using a mathematical model capable of solving both the velocity and pressure fields in the experimental vein, respectively.
In order to control the phenomenon as well as to reduce possible errors, the study has been divided into two main stages, starting from the realization of a simplified axisymmetric model, and finally creating a three-dimensional complex model.Each of the two steps will in turn include several intermediate phases that will gradually increase the complexity of the model with the addition of new elements or new modelling parameters.By increasing the complexity of the models, the time to solve the equations that describe the turbulent flow also increased.
This study aims to identify the optimal way to place non-invasive measuring instruments for the determination of field velocities in the experimental vein.

The experimental facility
For the reproduction of the flow to different parameters under laboratory conditions, the experimental installation (aerodynamic tunnel) built for testing multiphase flow of air and sand inside of harsh environment, located in the EOL facility of the National Institute of Aerospace Research Elie Carafoli, has been used.The aerodynamic tunnel include a fan with frequency converter, heating battery, experimental vein, openings for dust/sand particles and two cyclone separators.

Fig. 1. Harsh environment testing facility
The central section of the tunnel will be modified to include a modular design, composed of three sections, t

•
A section with 1000 mm length (Test section) made of a transparent tube (to be able to measure with optical instruments) with an inner diameter of 400 mm; • Two identical metal tubes of 1470 mm size with a 400 mm inner diameter.(Central section) The three sections will be assembled with metal flanges and insulated with rubber seals.The top of the diffuser will be provided with an access hatch for inserting the models into the experimental vein.The models will be attached to the diffuser base on a specially built bracket.

Methodology
Turbulent structures involve both a non-stationary movement and a refined degree of computation grid.Compromises must be made to achieve consistent results.
To choose the computational models that characterize the problem as close as possible to reality, the ANSYS Fluent documentation was studied on the hybrid system LES / RANS.After determining the pattern that will define the flow, next step is the definition of the fluid properties.
The calculation models for the turbulence, chosen for the LES / RANS hybrid simulation of the flow, was the DDES with Realizable k-ε.
The DDES (Delayed Detach Eddy Simulation) turbulence model is a RANS/LES model hydride that offers flexibility and a convenient scale to simulate fluid velocity with large Reynolds numbers.
In the DES approach, the unsteady RANS models are employed in the boundary layer, while the LES treatment is applied to the separated regions.The LES region is normally associated with the core turbulent region where large unsteady turbulence scales play a dominant role.In this region, the DES models recover LES-like subgrid models.In the near-wall region, the respective RANS models are recovered.[1] DES models have been specifically designed to address high Reynolds number wall bounded flows, where the cost of a near-wall resolving Large Eddy Simulation would be prohibitive.The difference with the LES model is that it relies only on the required resolution in the boundary layers.The application of DES, however, may still require significant CPU resources and therefore, as a general guideline, it is recommended that the conventional turbulence models employing the Reynolds-averaged approach be used for practical calculations.[1] The DES models, often referred to as the hybrid LES/RANS models combine RANS modeling with LES for applications such as high-Re external aerodynamics simulations.In ANSYS FLUENT, the DES model is based on the one-equation Spalart-Allmaras model, the realizable k- model and the SST k- model.The computational costs, when using the DES models, is less than LES computational costs, but greater than RANS.[1] The k-ε model requires a value of y + parameter that varies between 30 and 300 to produce results with a correct physical significance.This allowed to create an uniform grid with a relatively small dimension of only 5 mm, which was derived into a value of y+~138, for the 2D axisymmetric and y+~168 for the threedimensional model.
The Standard k-ε is a model based on the transport model equations for the turbulent kinetic energy (k) and the dissipation rate (ε).The model used for the initial solution calculation, Realizable k-ε model, differs from the Standard k-ε model in two ways:

•
The Realizable k-ε model contains an alternative formula for turbulent viscosity.

•
The modified transport equation for a dissipation rate ε, was derived from an exact equation for transporting the square mean values of the fluctuations of vorticity.In general, Realizable k-ε is used in favour of the Standard k-ε model because gives correct results for a much larger number of problems [1].
The spectral synthesizer, implemented in ANSYS Fluent, provides an alternate method of generating fluctuations for each velocity component.It is based on the random velocity generation technique originally proposed by [8] and modified by [2].In this method, the components of the fluctuating velocity are calculated by synthesizing a zero divergence speed field by summing a finite number of harmonics Fourier.
Essentially, Kraichnan's method of generating a field of velocity is the following form: where: t is the time, ωn is the frequency, x ≡ (x, y, z) is the spatial coordinate, kn ≡ (kx, n, ky, n, kz, n) is the vector of the form; vn and wn are amplitudes determined by the velocity vectors.Generating a non-stationary flow field using the modified Kraichnan's method [2]: where: l, τ represents the length scale of the turbulent flow and the temporal scale of the turbulence; εijm represents the Levi-Civita tensor; N (M,σ) a normal distribution with the mean value M and the mean square deviation σ;    and n represents a sample of n waveforms and frequencies of the model turbulence spectrum.

Stages of simulations
Virtual tunnel geometry was done using CAD ANSYS SpaceClaim, which is part of an ANSYS 18.2 suite.We have simplified this geometry and introduced it into the ICEM CFD 18.2 program, where we created the appropriate meshing network.Numerical simulations were performed and interpreted in the ANSYS FLUENT 18.2 expert program.
In order to control the phenomenon as well as to reduce possible errors, the study has been chosen to be divided into two main stages.

Axial flow simulation
Phase I: After creating the geometry with SpaceClaim, we have computed the calculation grid using ICEM CFD considering a rectangular cell network of the same size (5mm × 5mm) because the DDES model requires a uniform computation grid to generate consistent results.Finally, the grid included 91,908 cells, 93,310 nodes and 185,217 dots, divided into 8 sections that define the inner fluid.
For this simulation, we used the k-ε turbulence model with Standard Wall functions.The inlet of the tunnel is the velocity inlet type, with the turbulence specified by the intensity of turbulence (TI = 5%) and the hydraulic diameter (Dh = 0.7 m).We chose the outlet type; the pressure of the airflow is variable.The speed at inlet boundary condition is 8.95 m/s.Air was introduced at the entrance to the tunnel ( = 1.225 kg/m 3 and dynamic viscosity =1.7894×10 -5 ).The tunnel walls were set using "no slip" boundary condition.
Regarding the methods of solving, we have chosen the SIMPLE algorithm for speed and pressure coupling.The SIMPLE algorithm uses a relationship between speed and pressure corrections to ensure mass consistency and to obtain the pressure field.
For spatial discretization, we used second order schemes for the first 2500 iterations.We could choose to use a 2nd order scheme from the beginning because the flow inside the simulated model is aligned with the calculation grid, and in this case its use leads to superior cell-type precision through an extension of the Taylor series characteristic of the solution centered on the cell.When choosing for second order schemes, the quantities on the cell faces are calculated using a multidimensional linear reconstruction approach.
The simulation has been left to run either until it stabilizes at a reasonable value or until it reaches a value under the 10 -6 threshold for the residuals.After 2500 iterations, the simulation reached convergence of residuals equal to 10 -14 .
In the second step we changed the SIMPLE algorithm in favor of COUPLED, for a faster convergence and we switched to transient simulation using the DDES -Realizable k-ε hybrid model.The discretization schemes used, were the second order for pressure and magnitude that characterize the turbulence and bounding central differencing schemes for the momentum equations.We modified the inlet condition by adding the spectral synthesizer to generate random perturbations using the turbulence intensity as the input parameter.
To compute the solution, we determined the time steps necessary for a CFL=1 as, Δt = 2×10 -4 s.For a particle to travel through the whole fluid domain, a number of 1332 steps are required.For a coherent analysis, we considered a simulation time to 5328 time steps (simulation time Δt=1.065s),thus a space traveled by a particle equal to 4 times the axial length of the fluid domain.

Three-dimensional flow
Phase II: To approach the three-dimensional model we followed the same modelling and analysis method as the 2D axisymmetric model case.However, due to the high computational effort to obtain the solution, the model was reduced to 1 m (corresponding to the test vein).The inlet boundary condition was reconstructed using the data computed in the 2D axisymmetric simulation.The model was also built in ANSYS Workbench 18.2 using the Space Claim, ICEM CFD, and FLUENT components.
After creating the geometry with the SpaceClaim program, we made the geometric grid using ICEM CFD and created a rectangular cell network with the same hexahedron sizing (5×5×5 mm).
The three-dimensional grid contained 2,099,200 cells, 2,135,625 nodes and 6,333,696 faces, divided into 7 blocks defining the inner fluid.
Initially, for this simulation, we also used the Realizable k-ε turbulence model with Standard Wall Functions.The boundary condition for tunnel's entrance is the velocity inlet type, with the turbulence specified by the intensity of turbulence (TI = 5%) and the hydraulic diameter (Dh = 0.4 m).
The speed on the inlet condition was obtained by extracting a speed profile corresponding to a final step of the 2D axisymmetric simulation and processed using the rotation matrix.The rotation matrix for 2D geometry has the form: The new system coordinates (x ', y') of points (x, y) can be written as: Air was introduced at the tunnel inlet (=1.225kg/m 3 and dynamic viscosity η=1.7894×10 -5 ).The tunnel walls were set using default condition.Outlet condition is outflow type.
The solving methods are identical to those used in the 2D axisymmetric simulation.
The difference between the 2D axisymmetric simulation and the 3D model simulation was the number of iterations at which the solution reaches convergence and the simulation time.
For the 3D model, the simulation reached convergence of residuals equal to 10 -14 after 1500 iterations and the simulation time where the , Δt = 0.16s, where the particle the space equal to 4 times the axial length of the fluid domain, correspond to 800 time steps.

Methods for vortex identification
The identification of vortices is a necessary procedure for understanding the phenomes that occur in complex fluid flows.It allows, for example, being able to track individual vortex in a turbulent flow, from creation and dissipation, and it has a crucial role in the precise understanding of structures existing in the fluid.
The difficulty of finding a universally accepted definition has led to the development of multiple detection methods.These are generally based on the results of the simulations or speed measurements in the experiments, and most often use a field vector of velocities.
Turbulent airflow analysis in our three-dimensional model is structured based on the article published in [3] using computational techniques, visualization methods and vortices identification.
Identification of vortex regions is a complex technique and can be put into practice using different criteria such as local minimum pressure, trajectory and current lines; swirl motions defined by the strain rate tensor [4]; the turbulence intensity of filtered velocity field [5].

Streamline and pathline
Lugt [1979] proposed the use of streamline or closed pathline to detect vortices.Applying this method to the three-dimensional model at two different time steps, we obtain the following results:  This method fails to produce the desired results.An inadequacy of this definition is that a particle cannot undergo a complete rotation motion around the vortex (so there is no closed line of the current line) during the wake of the vortex.This happens when the vortices suffer a transition due to nonlinear processes such as fusion, tearing or disintegration before swirl fluids can exhibit a complete circular motion.[3].

Local minimum pressure
This criterion is based on the following principle: in a vortex, the pressure tends to have a local minimum along the axis generated at the center of a swirl motion when the centrifugal force is balanced by the pressure force [3].This is related to the second condition of [6], namely the pressure in the vortex core is smaller than that near the boundary of vortex.In the Figures 5,6 we can see the implications of this criterion on the section modelled by isosurface with a constant minimum pressure of -10 Pa colored by the swirl motion values.
Bellow we can see the most significant results: This principle is strictly true for steady fluid flows.A defined local minimum pressure may exist in a turbulent motion as can be seen in Figures 5, 6, which does not necessarily involve a vortex.Local minimal pressure is only a consequence of the strain rate of the fluid layers.
The swirl effect of the airflow is caused by the interaction between the axial velocity fields and the velocity fields created on the model's circumference, the pressure playing a crucial role in forming it.While the magnitude of the radial velocity decreases in the axial direction, the radial pressure gradient also decreases.For a rapid rotation identification, we can consider a negative pressure gradient that locally biases the flow, resulting in a flow direction opposite to the average flow direction [7].
We may see this effect by plotting magnitude velocity profiles extracted at different times, separated with a time step equal to Δt= ni+1+0.01s(where ni = time step -0.04s).Fig. 7. Magnitude velocity profiles at different time steps at 0.5m in experimental vein.

Q criterion
Q criterion, by the condition Q > 0 (where Q is the second invariant of velocity gradient tensor ∇), implies a local pressure lower than the ambient flow pressure.By applying Q criterion, we obtained a range of values between {min (-1.6 × 10 16 ), max (4.3 × 10 8 )}.For the three-dimensional model, we used a positive Q value of 10,000 that we modeled as isosurface to determine if the condition imposed by this criterion for vortex identification is met.Increasing Q reduces the complexity of the isosurface, but a higher value The Q criterion is met and this can indicate that there are vortices in the simulated flow, but the use of Q criterion may be inadequate for detecting vortices such as conically symmetric vortex or axial symmetric vortex that passes center of vortex ring.To detect the above mentioned vortices, the use of λ2 criterion may be more precise to generate accurate vortex core geometry.

λ2 criterion
The λ2 criterion starts from the premise that the minimum local pressure is not a sufficient criterion for detecting vortices.Inadequacy of this criterion is generated by the strain rate of a fluid that can create minimal local pressure without generating a vortex.
By Figures 9 to 14 we have pointed out that the eigenvalues, λ1, λ2 and λ3, of symmetrical tensor S 2 + Ω 2 , meet the criteria for the identification of vortices.To identify a vortex, this method imply the following conditions: • Condition I :  1 ,  2 and  3 are eigenvalues of symmetrical tensor S 2 +Ω 2 ; • Condition II :  1 ≥  2 ≥  3 ; • Condition III :  2 < 0 ; In Table 1 we can observe the possible choices of eigenvalues λ1, λ2 and λ3 as well as the difference between the definitions based on the positive value of Q criterion and the negative value of the own value λ2 [3].We plotted the isosurface of the positive value of Q criterion with a constant value of 10,000 to represent correctly the geometry of vortex core, colored by the value range of λ1, λ2 and λ3 in an attempt to validate the concept of vortex identification according to criterion λ2.    1 is validated.To fully validate criterion λ2, we extracted the last range of eigenvalue λ3, using in an identical procedure as in the case of eigenvalues λ1 and λ2.In order to observe from a numerical point of view whether Condition II is met, we decided to structure the maximum and minimum eigenvalues in Table 2. • The condition that λ1, λ2 and λ3, being eigenvalues S 2 + Ω 2 of is fulfilled, • The condition  1 ≥  2 ≥  3 is met by comparing the maximum and the minimum eigenvalues of S 2 + Ω 2 shown in Table 2.
• The condition that λ2 is the second largest value and has a negative value, is met for the simulated model.

Hairpin vortex
Predominantly, hairpin vortices, occur in discussions about the evolution of the turbulent boundary layer [9]; and is conceptualized as a condensed region of rotational motion that has a look that reminds of a "hairpin".Such objects having a rotational movement tend to have "legs" which are in fact two parallel vortices that spins in different direction and they are connected by an arc or horseshoe shape.
To validate the assumption that the hairpin vortex is formed in the boundary layer of the three-dimensional model, we represented a section of the threedimensional model, sectioned in the longitudinal direction to represent the turbulent boundary layer thickness and in the cross direction, with an axial length of 0.5 m, which indicate the downstream side of the pipe.
Bellow we can see the most significant results:  A velocity different from zero induces a cane shape of the vortex resulting in a "leg" stronger than the other, which spins in the opposite direction, around its own parallel axes, and the asymmetry is known to be a probable condition.[10] We may see this effect by representing a cross plane through the structure of a hairpin vortex where the magnitude velocity values and velocity vectors' field are represented in Figure 18.

Conclusion
The study of vortex in turbulent boundary layer is a solid research starting point for several reasons: firstly, in a boundary layer, any vortex with a direction other than the normal wall direction has the potential to function as a "pump" that transports mass and impulse over the mean velocity gradient.Second, the concept of vortex is invaluable as a conceptual for a complex class of threedimensional movements.
Vortices are among the most coherent of turbulent movements and tend to be persistent in the absence of instability.Finally, powerful vortex function as a source for pressure fluctuations because of their low-pressure core and high pressure regions that they induce in nearby fluid streams.
The first step in the turbulence study was testing the various identification criteria known in the literature using appropriate visualization techniques.
The three-dimensional model is a "validation" of the axisymmetric simulation, representing the reality, but also by being more complex, we managed to obtain results that confirm initial theories.
Of the four criteria used, it can be concluded that only the definition based on 2 criterion represents the topology and geometry of the vortex core correctly for the large variety of fluid flow.Although, in most cases, Q criterion and 2 criterion similarly determine the core of a vortex, each of criteria can generate different results for vortices formed in fluid domain.
An expansion of the simulated model and experiments with more advanced equipment is possible, and is required to study at a very different level the phenomena that occur during a turbulent flow.

Fig. 9 .
Fig. 9. Isosurface of the constant positive value of Q colored by the λ1 values -time step Δt= 0.08s.

Fig. 10 .
Fig. 10.Isosurface of the constant positive value of Q colored by the λ1 values -time step Δt= 0.16s.It can be seen that the λ1 eigenvalue range of values is positive, indicating the validation of the first condition of Table 1 for identifying the vortex core.Using the same procedure as for eigenvalues of λ1 we extracted the range of values of eigenvalues of λ2 represented in Figures 11, 12.

Fig. 12 .
Fig. 12. Isosurface of the constant positive value of Q colored by the λ2 values -time step Δt= 0.16s.Eigenvalue λ2 has negative values and thus another condition of Table1is validated.To fully validate criterion λ2, we extracted the last range of eigenvalue λ3, using in an identical procedure as in the case of eigenvalues λ1 and λ2.

Fig. 13 .
Fig. 13.Isosurface of the constant positive value of Q colored by the λ3 values -time step Δt= 0.08s.

Fig. 14 .
Fig. 14.Isosurface of the constant positive value of Q colored by the λ3 values -time step Δt= 0.16s.The value range obtained for λ3, quantified in Figures 13, 14 has negatives values and thus validates the condition that λ1, λ2 and λ3 are eigenvalues of S 2 + Ω 2 (Condition I) and condition that the eigenvalue λ2<0 (Condition III).In order to observe from a numerical point of view whether Condition II is met, we decided to structure the maximum and minimum eigenvalues in Table2.

Fig. 18 .
Fig. 18.Cross plan represented through the structure of a hairpin vortex, time step Δt= 0.144s.