[zx-81]

Hi all,

I have developed a new incompressible fluid solver targeted for doing very high-resolution simulations of turbulent flows on standard workstations.

The following video shows a turbulent jet, simulated with one billion degrees of freedom (1536x768x768 grid) at about 3 minutes processing time per frame on a standard intel 16 core work station. (watch video in HD if possible)

https://youtu.be/DjjsjYLg89U

Another example shows a turbulent 3D flow (passive marker concentration) over a square cylinder (Grid resolution 1024 in longest direction)

https://www.youtube.com/watch?v=1NC6I6W0mQs

The kinetic energy spectrum (computed over 2D slices) shows that all scales down to the nyquist limit of the grid are resolved (solid line is the -5/3 Kolmogorov slope)

https://www.youtube.com/watch?v=0R59Cq9oVqM

Note the small 'pile-up' of energy at the highest-frequency end of the spectrum (nyquist-limit of the grid) which is due to the simulation using nominally zero physical viscosity. It also shows that the solver is practically free of numerical dissipation.

With first results looking promising so far, I would now like to do some 'serious' verification and benchmarking in order to compare the accuracy and performance of the solver to established standard solutions.

Large parts of the solver are still in 'research & devlopment mode' and I keep trying out new ideas to further improve the algorithm and code implementation. So I would like to be able to readily verify and quantify the degree of improvemnt of newly implemented modifications.

It is straight forward to do such testing for single components of the solver (i.e. the pressure poisson solver for which one can simply compute the error residual and number of iterations over a set of synthetic test cases as measures of accuracy and performance).

But unfortunately I have not yet found a satisfactory way to do this kind of accuracy verificationfor the whole simulation of high-resolution turbulent flows.

I read the standard papers and links on this site about CFD Verification & Validation but those do not seem quite exactly address my problem.
The problem is that for fully turbulent flows (i.e. deterministc chaos) there is no fixed 'point-to-point' reference solution to compare the simulated flow solution against.

Of course one can compute averaged/integrated quantities (drag coefficients etc.) but this does not tell me about the quality of the highly resolved fine structure of the turbulent flow field. (But the quality of this fine structure is exactly what I a most interested in)
For example there could be numerical artefacts like 'checker board oscillations' present in the fine structure of the solution but the integrated/averaged measures of quality could possibly still look ok.

So I would like to ask this community for any further ideas. My background is not classical CFD but software engineering and computer graphics so maybe I'm missing some input from CFD experts on this.

Any advice or comment will be greatly appreciated !

[arjun]

When we wrote our DNS-LES solver in cylindrical coordinates, we used smooth sphere cases to validate results.

So I would just try to see if I could predict drag crisis properly.

[FMDenaro]

Originally Posted by zx-81

Hi all,

I have developed a new incompressible fluid solver targeted for doing very high-resolution simulations of turbulent flows on standard workstations.

The following video shows a turbulent jet, simulated with one billion degrees of freedom (1536x768x768 grid) at about 3 minutes processing time per frame on a standard intel 16 core work station. (watch video in HD if possible)

https://youtu.be/DjjsjYLg89U

Another example shows a turbulent 3D flow (passive marker concentration) over a square cylinder (Grid resolution 1024 in longest direction)

https://www.youtube.com/watch?v=1NC6I6W0mQs

The kinetic energy spectrum (computed over 2D slices) shows that all scales down to the nyquist limit of the grid are resolved (solid line is the -5/3 Kolmogorov slope)

https://www.youtube.com/watch?v=0R59Cq9oVqM

Note the small 'pile-up' of energy at the highest-frequency end of the spectrum (nyquist-limit of the grid) which is due to the simulation using nominally zero physical viscosity. It also shows that the solver is practically free of numerical dissipation.

With first results looking promising so far, I would now like to do some 'serious' verification and benchmarking in order to compare the accuracy and performance of the solver to established standard solutions.

Large parts of the solver are still in 'research & devlopment mode' and I keep trying out new ideas to further improve the algorithm and code implementation. So I would like to be able to readily verify and quantify the degree of improvemnt of newly implemented modifications.

It is straight forward to do such testing for single components of the solver (i.e. the pressure poisson solver for which one can simply compute the error residual and number of iterations over a set of synthetic test cases as measures of accuracy and performance).

But unfortunately I have not yet found a satisfactory way to do this kind of accuracy verificationfor the whole simulation of high-resolution turbulent flows.

I read the standard papers and links on this site about CFD Verification & Validation but those do not seem quite exactly address my problem.
The problem is that for fully turbulent flows (i.e. deterministc chaos) there is no fixed 'point-to-point' reference solution to compare the simulated flow solution against.

Of course one can compute averaged/integrated quantities (drag coefficients etc.) but this does not tell me about the quality of the highly resolved fine structure of the turbulent flow field. (But the quality of this fine structure is exactly what I a most interested in)
For example there could be numerical artefacts like 'checker board oscillations' present in the fine structure of the solution but the integrated/averaged measures of quality could possibly still look ok.

So I would like to ask this community for any further ideas. My background is not classical CFD but software engineering and computer graphics so maybe I'm missing some input from CFD experts on this.

Any advice or comment will be greatly appreciated !

Hello,
If you want to seriously validate your code, first you need to address the kind of formulation you are using to focus on the paremeters more relevants. For example, looking at your energy spectra I exclude you are running a DNS case as the dissipative part of the spectrum is totally missing. Therefore, if you are using LES you should detail the filter size and the SGS model you are using.

Homogneous-isotropic case is the classical test-case.

The first sever complex-flow test-case I would suggest is the plane channel flow, many DNS database at several Re_tau numbers are available. The backward-facing step is further available.

Let us know

[zx-81]

Originally Posted by FMDenaro

Hello,
For example, looking at your energy spectra I exclude you are running a DNS case as the dissipative part of the spectrum is totally missing.

This test case was simply run with zero physical viscosity to verify that the solver is free of numerical dissipation (Any remaining energy attenuation below the Kolmogorov line would be due to numerical viscosity)
It could of course be turned into a DNS any time by specifying a non-zero viscosity.

The goal of the solver is to resolve as many wavenumbers as possible within the inertial subrange (this is where he turbulence fine structure develops) down to the nyquist scale of the underlying grid. (The dissipative part of the spectrum is less intersting because it can be relatively easily simulated by, sloppily speaking "just blurring away the fine structure" generated in the inertial subrange).

So once that "resolving power" of the solver is verified, it could then be used as a proven basis solver engine for DNS, LES, implicit LES etc.

[FMDenaro]

Originally Posted by zx-81

This test case was simply run with zero physical viscosity to verify that the solver is free of numerical dissipation (Any remaining energy attenuation below the Kolmogorov line would be due to numerical viscosity)
It could of course be turned into a DNS any time by specifying a non-zero viscosity.

The goal of the solver is to resolve as many wavenumbers as possible within the inertial subrange (this is where he turbulence fine structure develops) down to the nyquist scale of the underlying grid. (The dissipative part of the spectrum is less intersting because it can be relatively easily simulated by, sloppily speaking "just blurring away the fine structure" generated in the inertial subrange).

So once that "resolving power" of the solver is verified, it could then be used as a proven basis solver engine for DNS, LES, implicit LES etc.

The issue is much more complex ... in your case, due to the inviscid flow assumption the cascade is terminated at the level of the Nyquist cut-off by the numerics but the shape of the spectrum close to the cut-off can deviate from the inertial range (and however, the inertial range is a good reference only in case of homogenoeus flows) not only for the numerical dissipation. You are running a case in which you have unresolved terms, in this sense you have a LES solution. But I don't know (you did not give us any detail) how your code works (FD,FV,FEM,SM) and if you have an explicit SGS model. Therefore, you must be aware that you have an "imperfect" solution as the resolved large scales interact with the unresolved large scale across the numerical cut-off. We call that "no-model LES".
I strongly suggest to run the 3D periodic (forced) homogeneos turbulence and plot the total kinetic energy versus time to see if it remains constant.

[zx-81]

Thanks for the suggestion. However, as far as I know the drag crisis effect will show only for very large Reynolds numbers > 1e5 which is too large for my DNS solver which should be able to resolve a maximum of about Re ~ 1e4.

Anyway. Your response highlights the problem I'm having with such kind of tests that will tell me if the algorithm/code works correctly in predicting "total"/average kind of effects.

However, such tests (drag coefficient etc.) does not tell me if the 3D fine structure of the turbulent velocity field is simulated correctly. For example, there could be numerical problems leading, e.g. to small checkerboard-like high-frequency distortions but global/integrated values such as drag coefficients etc. could still look correct.

So it seems that all I can do for now is to 'just look' at high-resolution 3D renderings of the velocity field or other quantities (vorticity, marker concentration) and see if the turbulent fine structure "looks plausible".

Of course there cannot be any kind of 'exact reference solution' for turbulent velocity fields, because (due to determinitic chaos), each individual realization will be totally different in it's fine structure.

So how can one verify that the instantaneous fine structure of the velocity field is resolved correctly then ?

[FMDenaro]

No there is no meaning in comparing pointwise/istantaneous values. You must compare only the statistics

[zx-81]

Thank you for the reply.

Originally Posted by FMDenaro

But I don't know (you did not give us any detail) how your code works (FD,FV,FEM,SM) and if you have an explicit SGS model.

The code is based on a standard fractional-step finite-difference scheme.
In order to be able to perform hi-resolution 'quasi-DNS' simulations on normal work stations, the code does contain numerous performance optimizations such as using a co-located grid with approximate pressure correction and optionally some 'heuristic sub-grid turbulence' approximations, such as the vorticity confinement method (VC) and some other 'energy re-injection' heuristics that I'm trying out.

So in order to quickly try out new ideas for performance and/or accuracy improvement of the code, it would be nice to have a suitable 'measure of quality' of a fluid solver.

I already use some of the suggested verification methods based on global statistics of the flow field. Those are great in telling if the code is correct over all. But I'm still missing some measure of 'the degree of quality' or 'resolving power' of the whole flow solver.

For example, it would be very interesting to know how much the approximate pressure correction of the co-located grid degrades the quality of the flow field compared to the more exact staggered grid approach. (I implemented both methods and there is no noticeable difference both visually and in global statistics)

As another example, I once had a coding bug that caused unnecessarily high numerical diffusion in the otherwise correct scheme. As the bug only affected the diffusion, most global statistics would still look ok. The problem showed up only after visual inspection of the velocity field which 'looked too blurry'.

Another related problem is how to verify the order of convergence for the whole solver. This is straight forward to do for single components, i.e. the pressure poisson solve where on can compare a reference solution against successively refined grids to check the rate of spatial convergence.
But how could this be done for the whole solver (lacking a point-to-point reference because of the turbulent random nature of the solution)

Anyway: I have now found a good introductory document about CFD verification http://sokocalo.engr.ucdavis.edu/~je...ano-Hirsch.pdf
that I'm studying right now. It contain some interesting ideas such as the "manufactured solution" approach that I think worth trying out.

[FMDenaro]

Maybe you would not like what I am going to write but that is my idea of developing, validating and using a code for turbulence.

First, with regards to the DNS, there is no necessity to have a high order discretization. In principle, you can perform an accurate DNS even with a first order scheme. That is due to the fact that the DNS requires that the Nyquist and the Kolmogorov frequencies are equal. At this level of resolution, the presence of a numerical diffusion lies within the range of the physical diffusion.
Much more relevant is the accuracy in LES but I do not know if you are interested in this formulation.
Therefore, I think that the high accuracy in a new code requires only the use of analytical test, this way you can validate the real accuracy order.
If you are interested, you can read Sec.7 here:
https://www.researchgate.net/publica...taggered_grids

the paper focusing on the projection method on colocated grid.
After the validation versus analytical test, the simulation of turbulence requires however the use of the statistical analysis. If you are interested in the verification of an unsteady case in turbulence you can try to replicate the flow problem described in Sec.4 here https://www.researchgate.net/publica...physical_flows

[arjun]

In addition

https://people.sc.fsu.edu/~jburkardt..._3d_exact.html

I can share the c++ code that you can use. I used that to benchmark and code FVUS/Wildkatze flow solver.

Originally Posted by FMDenaro

Maybe you would not like what I am going to write but that is my idea of developing, validating and using a code for turbulence.

First, with regards to the DNS, there is no necessity to have a high order discretization. In principle, you can perform an accurate DNS even with a first order scheme. That is due to the fact that the DNS requires that the Nyquist and the Kolmogorov frequencies are equal. At this level of resolution, the presence of a numerical diffusion lies within the range of the physical diffusion.
Much more relevant is the accuracy in LES but I do not know if you are interested in this formulation.
Therefore, I think that the high accuracy in a new code requires only the use of analytical test, this way you can validate the real accuracy order.
If you are interested, you can read Sec.7 here:
https://www.researchgate.net/publica...taggered_grids

the paper focusing on the projection method on colocated grid.
After the validation versus analytical test, the simulation of turbulence requires however the use of the statistical analysis. If you are interested in the verification of an unsteady case in turbulence you can try to replicate the flow problem described in Sec.4 here https://www.researchgate.net/publica...physical_flows

[Martin Hegedus]

Not sure of the ins and outs here, but one possibility is to run it with a slip boundary condition. The solution should then become irrotational.