The influence of the numerical solver selection on the nozzle impulse flow simulation results

Numerical simulations are currently used for different applications in a various fields of science. Certain solutions are not as obvious as the others while the results can give very valuable conclusions. Computational Fluid Dynamics (CFD) is one of the tools that can be used to solve different problems related with the mass and heat transfer. Nowadays it is already known that the impulse flow simulation allows to determine pressure pulsation attenuation parameters by a given geometry. However, the nozzle shape optimization method strongly depends on the numerical results obtained from the impulse flow simulation. In commercial CFD software Ansys-Fluent the obtained results depends strongly on the chosen numerical methods, especially the spatial discretization method. This is the reason to use other software as a benchmark. Alternative software FlowVision was used to perform the impulse flow simulation for the same geometries to compare the results. As there is a different problem definition in both systems the calculations, accuracy and results differ from each other. The paper describes the numerical differences between solvers. Article contains discussion about obtained results and includes hints how to avoid mistakes when user change software, especially in solving unusual CFD problems..


Introduction
Numerical calculations used to simulate reality in the modern world are a typical engineering tool used to improve design process. Users often don't have knowledge about the calculations performed by the software. In the case of well-known problems, the experienced user is able to determine the value of the obtained results. Interpretation problems increase when a new problem is solved using numerical methods.
Using CFD software numerical calculations can give, for the same problem, different results which depends strongly on the solver algorithms. Most common problems described in the literature are related with the appropriate selection of the closing models. However, this is a matter related with the simulated phenomena physics and is not related to the calculation method. Numerical problems are related with the solution methods such as a spatial discretization schemes, time integration schemes etc. In the paper [1] the differences between implicit and explicit time integration schemes are described for the turbulent incompressible flows. Influence of the spatial discretization methods on the numerical results for different calculations are described in [2,3].
Verification of newly introduced CFD solvers using the established commercial software is very desirable and the results of such comparisons are regularly presented at international conferences [4]. Modelling the impulse flow, damping is not only the effect of calculated physical phenomena but is also a numerical dissipation result. Authors in [5] describes the interpretation problems resulting from numerical dissipation. This paper presents impulse flow damping comparison for different shapes using two solvers with the same flow parameters, boundary conditions and additional models (turbulence, wall functions etc.).

Solvers
In the both analysed systems different solvers can be used. In the Ansys-Fluent software there are two main solvers: Pressure-and Density-Based. FlowVision have three different approaches from which the preferred is the AST "implicit-new" solver. This approach also define the time and spatial discretization schemes. In this paper the impulse flow simulations for two solvers are compared: Ansys-Fluent pressure-based solver and FlowVision AST "implicit new" solver. Algorithms of selected solvers are presented in the figures 1. and 2. Linearization of the governing equations using both softwares use the implicit approach. In the FlowVision the "implicit" solver developed by Aksenov [6] is used.

Ansys-Fluent
In the Ansys-Fluent pressure-based segregated solver the individual governing equations for the solution variables are solved one after another [7]. This approach makes this solver memory-efficient, although achieving the convergence is more time-consuming. Because the impulse flow is not a complicated state of flow the SIMPLE pressure-velocity coupling algorithm is selected. Calculating the variables the QUICK spatial discretization scheme is used for all values. In the literature it is described that QUICK scheme gives most satisfying results for the nozzle impulse flow damping [2].
where Sm is an additional mass source, Sh is an additional volumetric heat source and F is an additional force defined by the user.

Flow Vision
The FlowVision AST "implicit-new" solver use velocity-pressure split method, where the Navier-Stokes equations are solved using conservative velocities at a given time step. Variables like density, pressure, temperature or velocity are stored in the cell centers, although the mass flow rates are determined at the cell faces. This is a new type of the algebraic solver, it is a combination of the Aggregation AMG, selective AMG and RParFBBS solvers [8]. The P-V split approach is widely described in the paper [6]. More details of the general solver algorithm and the discretization schemes are the company industrial secret. Continuity equation in the FlowVision takes the same form as in the Fluent software. However, momentum (4) and energy (5) equations are treated in a slightly different way [8].
General momentum equation takes the form: As can be seen there are additional expressions which include to the equation isotropic (D) and anisotropic ( ) resistance. Additional forces are: Felforce exerted by the electrostatic field, Fuser-additional user defined volume force and F which is calculated using equation (5) : Using FlowVision user can choose one of three energy equations implemented in the software : -equation for total enthalpy, -equation for thermodynamic enthalpy, -simplified equation for thermodynamic enthalpy, in which the convective term is absent.
For slow gas flows it's recommended to use the simpler energy equation (6), which is formulated via the thermodynamic enthalpy [8]: where Jq is a specific heat flux and Q letter variables are additional volume source of energy due to viscous dissipation (Qvis), radiation (Qrad) and user sources (Quser).

Numerical Model
Numerical analysis of impulse flow were performed for three different geometries. As a basic shape the straight pipe is selected, the Venturi orifice represents geometry with linear neck-shape and as a most complicated geometry the pipe with hyperboloid shape neck is chosen. Geometries are defined as a surface models for the axisymmetric Fluent calculations as it is presented in the figure 3. In the FlowVision software the STL geometry file must be imported into the program because there is no geometry design system. To perform axisymmetric calculations model must be 3D element created from profile turned by, at last, two degrees as it is presented in the fig. 4.  Investigated geometries for both simulations are designed in the ANSYS software.

Mesh generation
Mesh generation conducted in the Ansys-Fluent is well known process. User can define global parameters for the model discretization or can define local grid parameters. Software generate grid that meet all defined conditions. For QUICK discretization scheme the grid have to be Cartesian or equivalent (the subsequent cells must be specified). Mesh detail for the Venturi orifice is presented in the fig. 5. Mesh generation philosophy in the FlowVision software is very different than in a classical mesh systems. Program calculate the dimensions of the computational domain and creates box-shaped catesian grid around it. Later, only this part of the initial grid, which includes the model geometry is calculated. For the pipe element with hyperboloidal neck the initial grid detail is presented in the fig. 6. For the initial grid user can define only number of cells in each direction. In presented case 2 cells are defined in the depth direction.
FlowVision mesh automatic adaptation module is very sophisticated and allows to create very detailed grid without user experience. Simplifying: the grid adaptation is a consequence of the cell splitting and merging. In automatic mesh refinement the cell is divided into four smallest cells and the process is conducted due the appropriate dimension of the smallest cell is obtained as it is presented in the figure 7 [8].   Sliced cells must be adapted as they have different shapes. Adaptations define the splitting and merging parameters of those cells: -that cell belong to the specified object volume, -that cell is adjacent to the specified object surface, -that cell is adjacent to a surface that has certain boundary conditions defined.
Adaptation also impacts the surface or the volume of the object. The grid refinement is performed mainly by two types of adaptation. Basic types are splitting and merging, but program can use also additional adaptation types. One of the most important is the "adaptation to solution" which allows to adapt the grid in the object area where big gradients of some variable occurs. However the program has a limitation, especially for axisymmetic flow-the computational domain and the initial grid should be oriented parallel. This could be a problem if the model is not parallel to the coordinate system of the FlowVision workspace.

Boundary Conditions
Nozzle impulse flow in the numerical calculations is defined as a one-time step excitation with known value. In this case, the model consist of 4 boundary conditions for each software as it is presented in the fig. 9. Main difference between boundary conditions using Fluent or FlowVision is that in Fluent the axisymmetric calculations are performed on a "flat" model with one boundary defined as an axis and using the FlowVision there is a 3D turned model and the axisymmetric case is obrained by defining the symmetry BC on the side surfaces. Differences in the boundary conditions in the two compared softwares are neglible so they are not described in detail in the paper.

Turbulence model
Numerical simulations are performed using basic k-ε turbulence model. Ansys Fluent has three k-ε models: -Standard, -RNG, -Realizable. FlowVision software have implemented four models named: -KES (Standard), -KEAKN (Abe, Kondoh, Nagano), -KEFV (FlowVision), -KEQ (Quadratic), where the KES model is the default turbulence model used in this software. In the presented simulations the standard k-ε model, which is a default turbulence model in both system, is used. Software developers describes this models with some differences. As the paper have to concentrate on the grid and solver influence on the results, due to lack of space, the turbulence models are not investigated in details.

Wall functions
For some calculations proper choose of the wall functions model have significant influence on the final result. In the presented case the same wall functions models are defined. In the FlowVision two models WFFV and WFS can be chosen for the k-ε turbulence model. In the Ansys-Fluent software there is a more wall functions calculations methods defined. Except the standard method software have defined four more wall function calculation methods. In the presented investigations the standard wall functions in both softwares are used. In the FlowVision standard wall functions model is presented by the WFS model.

Results
In the presented investigations the working fluid is air modelled as an ideal gas. The environmental pressure is equal to 1 atmosphere and the fluid temperature is defined as a 300 Kelvin. In the presented case 20000 time steps were performed which gives the simulation time equal to 0,04 s. Calculations performed in both systems takes a similar time on the default calculation speed-up settings. As a simulation results the spatially averaged value of the mass flow on the outlet is investigated. For the impulse flow simulation such a result is a damped, oscillating function. In the figures 10., 11. and 12. the comparison between signals obtained in Ansys-Fluent and FlowVision software for three different shapes are presented.  As it can be seen in the figures the flow rate signal is damped in a different way in both software's. Ansys-Fluent software suppress the impulse flow much more than FlowVision software and this applies to all nozzles. In this case, can be concluded that the geometry shape and mesh definition don't have a significant influence on the flow damping. Analysing the graphs can be seen that the amplitudes of the forthcoming peaks are constantly falling in the Ansys-Fluent calculations and, that have rather pulsation-characteristic using FlowVision software. It is also visible that the flow is more intensive in the Ansys-Fluent software as the obtained signal have much more sharp peaks. As the pipe lengths, mass flow impulse value, fluid parameters: temperature and pressure are the same, very interesting is that the signal period is different in both systems. It means that the flow speed is different. It is possible that it's influenced by differences in the outflow boundary, minor gas parameter differences due to the accuracy, cell size differences, numerical dissipation and/or the model viscosity importance, which seems to be higher in the Ansys-Fluent software. In the figure 13., the detailed first amplitude of the twin hyperboloid nozzle flow signal for both investigations is presented. In the figures can be seen global similarities and local differences. For both softwares the maximum value is preceded by the function oscillation which consist of two peaks. In the FlowVision software the first mass flow rate value peak is bigger and the second one is smaller. The graph is smooth and free of sudden value drops, while in the Ansys-Fluent results can be seen this substantial behavior. Another difference which is not clearly visible in this case, is that the mass flow rate value rises faster for the FlowVision calculations, and falls slower, when in the Ansys-Fluent software is exactly the opposite. This is more visible for the empty pipe and Venturi orifice comparison of the first amplitude presented in the figures 14. and 15.  To get the answer how and why presented differences exist more detailed investigations of the numerical model are needed. In order to evaluate some differences between the obtained results, the pressure distribution maps study can be performed as a preliminary calculations assessment. In the figures 16.-18. the comparison of the pressure contours for investigated elements are presented. As there is no need to make quantitative comparison the pressure values are not presented to obtain better quality of the figures. Qualitative comparison allows to asses how much the pressure distribution differs in both cases and how it can affect the results. In the figures red colour represents higher pressure and blue colour represents lower pressure. For the empty pipe flow it's very significant that the under-and over-pressure occurs on different sides of the moving disturbance. In the figure 16. the disturbance moves backward after it is reflected from the outlet and it is moving into the left side. It seems that in the Ansys-Fluent the gas in expanded behind the disturbance when in the FlowVision gas is sucked thorught the outlet to follow disturbance. This is another indication that, to better understand the differences in obtained results, special attention should be paid to the outflow boundary condition. For complicated shapes the flow is not as simply as for empty pipe which means that analysis of the pressure contours cannot give much more information. In the figs. 17. and 18. Can be seen that the distribution of high and low pressure is in some places more and in another places less similar.  Due to lack of space in this paper the FFT analysis is not presented but it should be known that for the impulse excitation (dirac delta in a mathematic sense) the Fast Fourier Transform analysis has a fundamental meaning.

Conclusions and discussion
Using the impulse flow simulations as a tool to estimate the damping parameters of the different nozzles shapes requires a suitably prepared numerical model. Previous studies have shown that the impulse damping in the numerical simulations using Ansys-Fluent software overestimate the damping in the empty pipe and underestimate the damping in other nozzles. This attempt, to use another software for such a non-standard simulations, allows to be aware how much the results depends on the solver, and minor differences in the calculation algorithms. This paper presents differences and similarities between Ansys-Fluent and FlowVision software without assessment, which software is better, because in this case such an opinion is not possible.
Presented FlowVision software can be distinguished by a very interesting and effective grid generation philosophy which is often a weakest element in the other software. One of the conclusions of this paper is also that using commercial software, operator is not able to determine in 100% how the problem is resolved, especially if the problem is unusual. The author also want to leave the reader with an assessment whether one or another software will be better to calculate his specific problem and what is important, to pay attention to, when changing the calculation software.