Verification of commercial CFD codes using hydrodynamic instability

Code verification is an essential process in computational analysis. The lack of access to the source code in commercial softwares often makes this process very difficult to be performed. Usually, the use of manufactured solutions is not feasible and analytical solutions of the Navier Stokes Equations are not sufficient because normally many terms of the equations are null. The situation is more troublesome for unsteady flows since it is more difficult to find analytical solutions and errors can accumulate in time. Simulations focused on turbulent fluctuations or sound wave generation associated with aeroacoustics applications are even more important and difficult to perform, because even minute errors can have a large impact on the flow. Here it is proposed a code verification method based on the simulation of a mixing layer, a canonical flow that transitions to turbulence. This flow is linearly unstable and the stability calculation affords a semi analytical solution. The paper shows how this solution can be used for classification of code performance. In the study, several commercial codes have been tested and very different code performances were obtained, both in terms of accuracy and grid requirements.


Introduction
Code verification is an important aspect of code development as it deals with the identification of programming errors [1,2].However, it is the user's responsibility to ensure the correctness of the computational results.Therefore, code verification is often a procedure commonly undertaken by code users too.This is not only to check the correctness of the code, but also to check if the user is using the code correctly.
Several techniques of code verification have been presented in the literature.The method of manufactured solutions [1,3] is regarded as the most complete and accurate.It has been strongly recommended in the context of computational fluid dynamics.This method involves the selection of appropriate functions to represent the velocity and pressure fields (or other variables) which do not satisfy the equations of motion.These functions are substituted in the equations leading to a residue that can be calculated analytically.This residue is programmed into the equations as an analytical source term.The original functions become exact solutions of these equations of motion modified by the addition of the source terms.The production of an analytical exact solution is a key aspect of the method as it is then possible the study of the behavior of the truncation error and determine the numerical order of convergence of the solution.It is desirable that the order of convergence matches the nominal order of the lowest order algorithm used, but this is not always considered a requirement.What is important is that all routines be tested by the manufactured solution, including boundary conditions, and that the rate of convergence is sufficiently high.
Commercial codes are usually very restrictive on what a user can do.The possibility of changing the equations of motion by the introduction of an analytical source term is rarely an available feature.Users have to resort to other types of verification tests, such as the use of analytical solutions of the Navier-Stokes equations, for example.These solutions are obtained under several simplifying hypotheses that cancel many terms of the equation.As a consequence, and despite bringing important information, the test becomes incomplete.The advance of algorithms and computers begin to allow the use of commercial codes for the simulation of at least the large structures of turbulent flows.These flows are much more complex than the analytical solutions available for the Navier-Stokes equations.Tests based on analytical solutions become severely limited.General features of turbulent flows based on theoretical arguments are sometimes invoked to convey an indication of the quality of the code or the numerical solution.These are also limited since they are of statistical nature.For instance, the Kolmogorov 5/3 rule, although important, is based on dimensional analysis only.In fact, several terms of the equations could be missing, or the Reynolds number of the simulation could be wrong, and yet, the Kolmogorov rule would be satisfied.
When the interest is on turbulent flows, hydrodynamics instability can be called upon to help testing codes specially under the restrictive environment offered by some commercial codes.Moreover, large scales of turbulent flows are related to flow instabilities and originate from the growth of very small amplitude disturbances which poses extra challenges to algorithms and implementations.In this regard, this type of test brings information particularly relevant.The current paper shows how such tests can be done.It also shows how several restrictions typical of commercial codes can be circumvented.Results of several commercial codes were used to make the tests more realistic.The specific commercial codes were not revealed as the paper focused on the test and not on the quality of the codes itself.A highly accurate in-house code was also used as a benchmark.

Tests based on Taylor-Green Vortices
As said, analytical solutions of the Navier-Stokes equations can be used to test the accuracy of the codes [4].The Taylor-Green flow (1) is often used for this purpose.It is a fairly complex solution, with variation in space and time, and represents a difficult test.It is however restricted to a two dimensional incompressible flow and cannot be extended to compressible flows.One can think of using this solution to test a compressible code at very low Mach number.Such a test is here presented to illustrate its limitations.One of the codes tested in this paper was used to produce the Taylor-Green flow results.
The code here used required dimensional quantities, so the equations 1 were replaced by a dimensional version u(x, y, t) = U 0 sin(x) cos(y)e −2t v(x, y, t) = U 0 cos(x) sin(y)e −2t (2) where U 0 was taken as 10m/s.The viscosity ν was assumed 1m 2 /s to ensure the viscous terms were not negligible and the base pressure was 101325 P a.In the code, the Mach number used was 0.2, and, consistent with an incompressible solution, ρ was considered uniform and equal to 1.225kg/m 3 .Periodic boundary conditions were applied in a computational domain 0 ≤ x, y < 2π .The initial condition for the simulation was obtained from equation 2 for t = 0.For illustration, figure 1 shows contour plots of u, v, p from equation 2. Figure 2 a shows a comparison between the exact (eq.2) and the numerical solution at (x, y) = (0, π 2 ).The results appear to agree reasonably well.Figure 2 b gives the relative error and shows significant oscillations.These oscillations arise, for example, from the fact that the initial condition used is not an exact solution of equations 2 at M a = 0.2.This inconsistency leads to an acoustic wave, in the case, a stationary wave in view of the periodic boundary conditions.These oscillations decay but take about the same time for the Taylor-Green flow to decay.As a consequence, there is not a time window in which the oscillations are negligible and the fluid is not yet stationary.Under these circumstances, a verification test becomes very limited.A similar situation is observed for the other variables.
Because the oscillations are of a stationary wave character, there are wave nodes, such as (x, y) = (π/2, 0), where the wave oscillations are null.At these points the exact and numerical solution are much closer and it is possible that the truncation error associated with the spatial discretization becomes dominant.Under these circumstances, some grid convergence test may lead to a reasonable estimate of the order of accuracy of the code.The results of a such a test are presented in figure 3.For course grids the decay rate is of second order as expected from the nominal order of accuracy of the code.
For not yet very refined grid, the error builds up significantly.This is consistent with the growth of round off error.Because the equation involves second derivatives, the round off error is expected to increase with grid refinement at a second order rate.It is apparent that, within some limits, the error increases on average at the expected rate.The minimum error attained by the code for this flow for v at (x, y) = ( π 2 , 0) is about 10 -3 .It is expected that the discretization error of the other variables at other positions would be of at least a similar order of magnitude.Indeed, studies performed but not given here show exactly that.This error level is a relatively high and is thought to be associated with the fact that only single precision calculations are supported by this code.
The Taylor-Green tests were to a large extent consistent with the expectations and con-veyed important information about the code accuracy and correctness.Nevertheless, clearly the test suffered from serious limitations and a very selective approach of what to analise had to be exercised in order to extract some meaningful information.To what degree this approach also covered code misimplementation or misuse is uncertain.Moreover, the procedure did not test the compressibility and three-dimensionality, and does not indicate the consequences of the high minimum error obtained.Furthermore, if a code does not provide the possibility that an initial condition is given as a generic function of space, tests based on analytical solution may prove to be impossible to perform.

Brief introduction to hydrodynamics instability theory
A linear hydrodynamic stability analysis, as in solid mechanics systems, is usually performed via the normal mode decomposition method.The approach consists of investigating the response of the flow to small-amplitude disturbances.
This theory is well established and the results have been verified by experiments and computational simulations in amongst many fundamental flows.Hydrodynamics instability can be called to help in the verification of CFD codes.
In any situation, when a basic flow can be established, for instance, a parallel flow, the Navier-Stokes equations can be written in terms of disturbance variables and linearized to study small-amplitude disturbances developing around this basic flow.Thus, the flow field can be decomposed according to v(x, y, t) = ṽ(x, y, t) p(x, y, t) = p(y) + p(x, y, t) where u, v represent the velocity components in x− and y−directions, p is the pressure and the tilde indicates disturbance.The stability equation resulting from substitution of infinitesimal two-dimensional modal perturbations into the incopressible Navier-Stokes equations and linearisation about a parallel mean state is the well-known Orr-Sommerfeld equation [5,6]: ω is a complex frequency.A positive value of the amplification rate ω i indicates growth and the mode is temporally unstable, while ω i < 0 indicates damping, and the mode is stable.The theory can be applied to compressible flows.Likewise viscous effects can also be considered, which in the incompressible case leads to the Orr-Sommerfeld equations.The equations are more complicated, but the framework is the same.Prediction from of the linear stability theory for different Mach numbers and Reynolds numbers for different flows can be readily found in various books.Here to keep the analysis simpler, only the inviscid case is considered.

Code test based on Linear Stability Theory
In this section a code test procedure based on linear stability theory is presented.A very accurate 6th order in house code [7,8] was used as an example to indicate how accurate the numerical solution can be.Moreover, as will be seen, the complete access to the source code permits the implementation of source terms in the flow equations which enable a more accurate comparison between the numerical and the theoretical solutions.Among other limitations, commercial codes usually do not allow the implementation of such source terms, but there are ways of circumventing these difficulties.These are discussed in the next section.
To simulate linear stability theory, a base flow has to be chosen.Here a hyperbolic tangent profile was considered, as shown in figure 4.This is a canonical flow that has a well documented instability theory.This instability is associated with the formation of large coherent structures in turbulent shear flows.
The velocity scaling of this flow is usually given by the free-stream velocity while the length scale is given by the so called vorticity thickness (δ w ), which, as represented figure 4 is associated with the extrapolation of the derivative at the hyperbolic tangent profile inflexional point.In figure 4, both these scalings are 1 and the profile is given by ū = tanh 2y. ( The linear stability as described in the previous section considers a steady base flow.In the numerical solution, however, the base flow widens with time owing to transport effects in the vertical direction.Forcing terms can be added to the governing equations to cancel the vertical transport effects.These terms are ∂ 2 ū ∂y 2 , ∂ ū2 ∂y and ∂ 2 T ∂y 2 , where T denotes temperature.According to the linear stability theory, the streamwise boundary condition was periodicity.The streamwise domain was one disturbance wavelength.The vertical boundary condition is meant to represent a homogeneous boundary condition at infinity.Several boundary conditions can be used for this purpose.They vary essentially in terms of efficiency, meaning that the most efficient is the one that represents in the smallest vertical domain and effectively infinity domain.In the current tests free slip was used at the horizontal boundaries.
The initial condition is composed of the base flow plus the disturbance.The initial disturbance is periodic in the streamwise direction.As long as the disturbance is sufficiently small, in principle any shape could be used in the vertical direction, but it is interesting to use a shape close to the eigenfunction of the stability problem.In the current application, the initial disturbance was ṽ(x, y) = Aαe −σy 2 cos αx.
and were utilized to amplitude of perturbation A = 10 −4 m/s and σ = 2 for the Gaussian curvature.
The Mach number used in the simulation were 0.01, 0.4 and 0.8 at Reynolds number 500 .It is not necessary to use this wavenumber, in fact it is desirable to use a range of wavenumbers.It is only important that the modes are sufficiently unstable.
Figure 5 shows time evolution of vorticity.Initially on base flow vorticity is observed.As time increase, a vorticity oscillation is seen, the disturbance flow.At yet later times a vortex is formed which then decays owing to viscous effects.
The amplitude can be measured from any variable anywhere in the flow, but v at (x, y) = (π/α, 0) is a convenient choice for the initial condition here used.Initially the disturbance amplitude displays some irregular oscillations with are associated with an initial transient.During this transient the stable modes composing the disturbance decay while the single unstable mode grows and dominates the disturbance flow.This transient can be substantially reduced if initial disturbance shape in the vertical direction is very similar to the unstable eigenfunction.After this transient the disturbance exhibits the exponential growth as predicted by the linear theory.At later times, where the disturbance attain amplitudes of about 1% of the base flow, the phenomenon saturates.This stage is associated with the vortex formation in figure 5.It is apparent that the stage associated with the exponential growth displays amplitudes that are too small to be observed in figure 5.
The stage that is important for the numerical test here described is that of the exponential growth.Indeed, figure 6 shows that an excellent agreement between the numerical   7. Comparison between theoretical results and grid independent numerical solutions for several disturbance streamwise wavenumber using house code.Lines are the theoretical solutions, while the symbols are the numerical solution.
solution and the theory.
Prior to perform comparisons such that of figure 6, it is essential to verify that the solution is independent of both the computational grid and domain.The grid required was relatively course in view of the high accuracy of the code used.Similar tests were performed to check that the solution was independent of the vertical domain.
Tests can be performed at other wavelengths, and it is advisable to do that.Figure 7 shows a comparison of the linear growth obtained numerically and theoretically for several wavenumbers and Mach numbers.For these cases, the theoretical results cannot be obtained from equation 7, but are available at reference [9].In figure 7 the agreement with theory is also very good.

Strategies for circumventing some commercial code limitations
The test described in the previous section can provide important information regarding of the performance, including commercial CFD codes.Nevertheless, because commercial codes are usually not very flexible in several strategies may have to be used in order to reproduce the stability theory results.In the results presented in this paper, such strategies were actually used and the results then indicate what can be achieved with their use.
First of all, commercial codes rarely provide the means for the inclusion of source terms, similar to those used in the previous section to prevent the base flow from evolving in time.If that is the case, one can run the solution at very high Reynolds number, which controls the widening of the base flow.What is important is that the base flow remains constant in the time scale that the disturbance evolves.
From evolution of the base flow at a high Reynolds number simulated by a commercial code 1, the base flow does not change up to time t = 50, while the linear evolution of the disturbance takes only about t = 35.
For some commercial codes an initial condition cannot be given by a generic function of y.However, usually at least one can establish initial conditions from a combination of uniform streams.In that case a hyperbolic tangent profile can be produce, from a initial condition like From this initial condition, the diffusive effects in the y-direction will eventually produce a hyperbolic tangent profile with the desired length scale.It is important however to use a low Reynolds number to ensure that the flow is stable and that disturbance do not grow.The low Reynolds number has the extra advantage of producing the desired profile more quickly.The hyperbolic tangent profile so obtained can then be directly used as the base flow for a high Reynolds number simulation.
If the commercial code allows the introduction of a generic initial condition, an initial disturbances like equation 9-10 can be use.Otherwise, the stability can be triggered from the discretization errors that are always present.In this case, limited control over the disturbance wavelength is exercised and the disturbance that eventually dominates corresponds to the most unstable wavelength.For instance, for a large streamwise computational domain, it is expected that the most unstable spanwise wavenumber prevails.If the initial disturbances, that is, the discretization error, is sufficiently small, there will be a time range in which the most unstable mode is sufficiently larger than the other modes for them to have a negligible effect and yet the dominant mode be sufficiently small to satisfy linear theory.With double precision accuracy it is usually possible to find a sufficiently long time window with these characteristics.Yet, some control over the disturbance wavelength can the exercised by choosing an appropriate streamwise computational domain.In fact, if the streamwise computational domain is smaller than the most unstable wavelength, obviously this mode will not arise.Under these circumstances, the most unstable mode will have a wavelength corresponding to the streamwise computational domain.Results obtained under these circumstances are also discussed in the next section.

Hydrodynamic stability test of commercial CFD codes.
The object of the paper is not to establish the performance of specific commercial codes, but rather, explain how the test can be performed.In addition, it is important for the reader to be able to access the kind of information that the test can provide for existing commercial codes.This section is then presented in order to provide this information.The results were chosen to illustrate the different types of code behavior that were observed in the several tests performed for several codes.Results from three codes are presented, here refer to as code 1 to 3. In the several simulations, care was taken that numerical results were independent of the vertical domains used.
Figure 8 shows the amplitude evolution for different space and time resolution for code 1.The Reynolds number was 10 5 and the Mach number was 0.003.For this code, a generic initial condition could be introduced.Therefore, the base flow was given as equation 8 and the initial disturbance as equation 9-10.The streamwise wavenumber was α = 0.89.In the picture, the initial transient was removed from the plots.Results for grid III and time step of 0.25 are fairly independent of the resolution.It is seen that the grid independent results are closest to the theoretical solution.The results are not as good as that of figure 6 most probably in view that the in house code was tested to be 6th order accurate in space, while the current code is nominally second order only.The grid for convergence is also significantly more refined than that of figure 6, which can also be attributed to the lower order of accuracy.Nevertheless the numerical solution is quite good, if one considers that this error is at the disturbance level and that in the linear regime the disturbances are a small fraction of the base flow.
Similar results for code 2 are given in figure 9.The Reynolds number was 10 5 and the Mach number was 0.003.This code did not allow the use of a generic initial condition, and the hyperbolic tangent profile had to be produced as described in the previous section.Likewise, an initial disturbance also could not be introduced in a controlled way.The instability was triggered by the numerical error.The streamwise computational domain used corresponded to the streamwise wavenumber of the most unstable mode.The circumvention strategies used were successful and an exponential growth of disturbances was obtained along a wide time window.Again, the initial transient was removed from the plots.The grid convergent results achieved a level of accuracy somewhat similar to that of code 1, but required a substantially more refined grid.The test then indicated that a similar level of accuracy can be achieved, but at a higher computational cost.Because the results from code 2 eventually converge, there is an indication that code 2 has a higher truncation error than code 1, despite the fact that for both codes, second order accuracy is claimed.Figure 10 compares theoretical growth rates and numerical ones obtained by code 2 for several streamwise wavelength.Since code 2 did not allow full control of the disturbance, the test was limited to wavenumber higher than the most unstable one.The control of the disturbance wavelength could only be made via the streamwise computational domain, as discussed in the previous section.Despite these difficulties, the level of agreement between theory and computation is similar to that obtained for the most unstable mode indicating that the circumvention strategy used worked well.
Results for code 3 are given in figure 11 at M a = 0.4.Here again the initial condition could be fully controlled.The base flow and disturbance used was the same as the other codes.Results indicate that the solution becomes independent of the grid for a relatively course (similar to that of code 1), but the numerical solution failed to agree with the theoretical one.Figure 12 shows a comparison between theoretical results and grid independent numerical ones at a number of streamwise wavenumbers for code 3. The results so not agree with the theory for all the wavenumber range tested.Nevertheless the most unstable wavenumber as predicted numerically agree fairly well with theory.
For code 3 the error does not seem to be of a truncation nature because the results do converge for a course grid.Another important source of numerical error is round off and   indeed there is some further indication that this is the error origin.First of all, for this code it was not possible to reproduce the linear stability theory for disturbances as small as those used for code 1.Moreover, while the other codes used double precision accuracy, code 3 only provided single precision accuracy.If round off error were the dominant cause of error for this code, better results would be obtained for solutions of larger magnitude as compared to the disturbance flow.Indeed this is the case.Figure 3 was produced using results from code 3.They showed that for the large Taylor-Green vortices the numerical solution agrees fairly well with the analytical solution.Indeed the amplitude evolution in figure 11 never quite achieved an exponential growth.It does seem that the disagreement observed in figure 11 stems from the fact that at the amplitude levels simulated were too large and the disturbances were already saturating under the influence of the nonlinear effects.This conjecture is also consistent with results from figure 12.
In it interesting to mention in passing that the linear stability test also indicates the grid resolution requirements for the proper simulation of a vortex in a free shear flow.This can be an important information for Direct Numerical Simulation (DNS), Large Eddy Simulation (LES) or detached eddy simulation (DES).If the vortices are large, the Taylor-Green vortices can give an indication of the level of refinement needed, but if the initial evolution of the vortices is required this is better estimated by the linear stability theory.

Conclusion
Verification tests is a necessary step of code development and, indeed, code use.When it is possible to use the technique of manufactured solution, this should be used.In commercial CFD code this is rarely an option.Use can be made of analytical solutions of the transport equations, which can lead to important information such as the order of accuracy of the numerical solution.Test using analytical solutions are nevertheless limited.On the one hand, analytical solutions are available only for very simple flows, often restricted to two dimensional incompressible conditions.On the other hand, tests based on the linear stability theory can be more demanding.In the paper its was shown that a code that performed reasonably well for the analytical Taylor-Green flow, exhibited a more limited performance in the linear stability test.Yet, codes for which the an initial condition cannot be given as a general function of x, y and z, may prevent the used of tests based on analytical solutions.
The paper shows that linear stability theory can be used to classify CFD code performance.This test can be used for fairly generic types of flows, including three-dimensionality and compressibility.A number of strategies to overcome the commercial code restrictions were discussed and examples of their use was given.
For illustration, results of linear stability theory tests on a number of commercial codes were presented.The tests were able to classify the codes with regard to solution accuracy.They were also able to compare code efficiency, that is, grid resolution requirements for a certain level of accuracy.The test was also able to suggest the origin of the dominant discretization error, namely, the truncation or round off error.

Figure 1 .
Figure 1.Contour plots of a Taylor-Green flow.

Figure 2 .
Figure 2. a) Comparison of the exact and numerical solution for v. b) Relative error between exact and numerical solution for v.

Figure 3 .
Figure3.Estimation of the order of accuracy of the code based on a Taylor-Green flow solution for v.For reference, the symbols represent absolute error, the black line represents a second order of accuracy decay and the red line represents a second order growth.

Figure 5 .
Figure 5.Time evolution of the flow shown by snap shots of the vorticity at selected times.

Figure 6
Figure 6.a) Time evolution of the disturbance amplitude obtained by direct numerical simulation for v and comparison with the linear stability theory prediction in logarithmic scale.Line is the theoretical solution, while the symbols are the numerical solution.

Figure
Figure 7.Comparison between theoretical results and grid independent numerical solutions for several disturbance streamwise wavenumber using house code.Lines are the theoretical solutions, while the symbols are the numerical solution.

Figure 8 .
Figure 8.Comparison between theoretical and numerical solutions for a number of computational grids with 50, 100 and 200 points per wavelength, for code 1.

Figure 9 .Figure 10 .
Figure 9.Comparison between theoretical and numerical solutions for a number of computational grids for code 2. The streamwise wavenumber was α = 0.89.

Figure 11 .
Figure 11.Comparison between theoretical and numerical solutions for a number of computational grids for code 3.

Figure 12 .
Figure12.Comparison between theoretical results and grid independent numerical solutions for several disturbance streamwise wavenumber using code 3.