[ville]

Hi,
A proposed, new approach for discretizing the convection terms
for the Navier-Stokes equations is given in the paper below based on
implementation of the new SSD scheme in the OpenFOAM environment
together with the RK4 scheme and projection method:

Vuorinen V., Schlatter P., Boersma B., Larmi M., and Fuchs L., A Scale-Selective, Low-
Dissipative Discretization Scheme for the Navier-Stokes Equation,
to appear in Computers and Fluids.

• http://dx.doi.org/10.1016/j.compfluid.2012.09.022

The paper shows an "easy trick" to reduce numerical diffusion due to flux limiting
or upwinding. Normally, in unstructured codes, the diffusive error is O(dx)

and it makes the overall accuracy also O(dx) for turbulent flows. The SSD

scheme is able to maintain O(dx^2) as the overall accuracy and it makes the

diffusive error O(dx^3). This is shown in the paper. In case you want more

info or want to implement the scheme in e.g. PISO solver you can ask me here
for more info.



Best,
Ville

%%%%%%%%%%%%%%%%%%%%%%

Brief instructions for implementation at top-level code are provided below.
We use it mostly to stabilize implicit LES computations but it is generally

a good candidate for explicit LES simulations as it is very non-dissipative.



About implementing SSD-scheme (Vuorinen et al., to appear in
Comp.&Fluids, 2012) to OpenFOAM top-level code

1) A Laplacian high-pass filter can be defined using the function
Up = -(Delta/(2*pi))^2 laplacian(U), where cutoff wavenumber is
kcut = 2pi/Delta. Setting Delta = Nyqvist wave gives kcut = pi/dx which
is recommended for SSD. We implement the filter here in conservative/self-
adjoint form (see the publication).
Create fields for a filter (e.g. in createFields.H)


surfaceScalarField DxDyDz
(
IOobject
(
"DxDyDz",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionSet(0,2,0,0,0,0,0)
);

DxDyDz = Foam:ow(mesh.surfaceInterpolation::deltaCoeffs() ,-2.0);

scalar pi = 3.141592653589793;

DxDyDz = (filterScope*filterScope)*DxDyDz/(4.0*pi*pi);

%%%%%%%%%%%%%%%%%%%%%%%%%%

%additionally implement high-pass filtered fields for
%the quantities you are solving e.g. velocity and set the BCs
%in /0/ directory to be e.g. 0 for the fluctuation (I used 0 typically)
% For example for the momentum

volVectorField rhoUprime
(
IOobject
(
"rhoUprime",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
rho*U
);
rhoUprime.correctBoundaryConditions();


2) Now things are settled up and the simulation begins normally.
You have to compute the fluctuation before solving the NS-equations.
For example, what I do in an explicit density based RK4 code:


rhoprime = fvc::laplacian(DxDyDz,rho); rhoprime.correctBoundaryConditions();
rhoUprime = fvc::laplacian(DxDyDz,rhoU); rhoUprime.correctBoundaryConditions();
rhoEprime = fvc::laplacian(DxDyDz,rhoE); rhoEprime.correctBoundaryConditions();
pprime = fvc::laplacian(DxDyDz,p); pprime.correctBoundaryConditions();


3) The SSD scheme is designed and validated for FULLY explicit time
integration schemes. Normally the convection term in e.g. momentum eq
would be written as fvc::div(phiv, rhoU).
The convection term is re-written as

fvc::div(phiv, (rhoU+rhoUprime)) - fvc::div(phiv, (rhoUprime))

Note that if you discretize both terms with the same scheme
you are not doing anything different from fvc::div(phiv, rhoU) i.e.
fvc::div(phiv, rhoU) = fvc::div(phiv, (rhoU+rhoUprime)) - fvc::div(phiv, (rhoUprime))

The trick is to choose "linear" for the first term (which is smooth
because rhoU+rhoUprime is smooth; the fluctuation is taken away;
i.e. we low pass filter by writing rhoU+rhoUprime)
and "blended" for the latter term which adds upwinding to it.

That's it!

4) FOR IMPLICIT formulation we tried it as well with good experiences.

Now you have probably terms such as fvm::div(phiv, rhoU) .
Do not do anything to them BUT discretize them with "linear".
You have to add an EXPLICIT source term to the RHS of the equations
which is fvc::div(phiv, rhoUprimeLinear) - fvc::div(phiv, rhoUprimeBlended).
Here rhoUprimeLinear and rhoUprimeBlended are the same field
computed in the same way as rhoUprime i.e. with the HP-filter.
You just need to copy that to two fields to separate the terms.
As you might have guessed the "Linear" is discretized with "linear"
and "Blended" with "blended". Perhaps you need to change the
signs + -> - but that works, we tested it.

What effectively happens is that you get a sink of momentum on the RHS
of your eqs which is essentially of form ~|U|f''''(x) i.e. contains fourth
derivatives. You see this easily in the case of linear convection eq
from Taylor series expansion and the facts do not change overhere.
More info in the paper.

[FMDenaro]

Many thanks, I just downloaded the paper for a reading

[francesco_capuano]

Thanks, very interesting. Filtering and numerics have been seldom taken into account simultaneously in the past, at least to my knowledge. I recall a similar attempt in the paper

Bose, Moin & You, "Grid-independent large-eddy simulation using explicit filtering", Physics of Fluids, 2010

although their approach is somewhat simpler. They explicitly filter the velocity field (using a commutative filter) and then just take care that the resolved scales are not affected by the numerical error (simply checking the modified wavenumber diagram, for a given scheme).

Is there any chance that your new solver / scheme would be included within a future OpenFOAM release?

[ville]

Hi,
Just a link to the final version of the article:
http://dx.doi.org/10.1016/j.compfluid.2012.09.022

Thanks for the hint on the previous paper as well.

The scheme is so easy to implement in the top-level code that
I mostly see it as two steps:

0) add the lines shown above which create a surfaceScalarField that is used
in the hp-filter to take into account grid spacing variation
1) add a .H file where you evaluate the hp-filtered fields
2) add a source term to the governing equations.

That's why I'd see it as a 5 minute modification of any existing solver
to implement.

From my behalf the RK4-solver could well become a part of the OF distribution (extended) in the future. Would you have some instructions of the proper actions to take?

Best,
Ville

[francesco_capuano]

Dear Ville,

thanks for your reply, I'll try to follow your instructions and modify an existing OF solver.

Concerning the OF-ext distribution, unfortunately I have no idea about the procedure to follow, but I'm pretty sure that other members could help you, especially in the OpenFOAM subsection of the forum.

Regards,
Francesco

[FMDenaro]

I just had a look to the paper, what I see is that the definition of the fluctuation u' is actually the recovering only of the resolved frequency smoothed by a top-hat like filter. I mean, in no way you can define a fluctuation velocity possessing wavenumbers outside the Nyquist frequency. You just use a deconvolution-like approach to define the sub-grid scale resolution.
It is correct?

[ville]

Hi,
I'm not doing deconvolution as we are limited to the grid scales i.e. k<=kNyq.
So the argument is that we are not resolving well the waves
at the grid scale ~2*dx. Thereby, it makes sense to
split the resolved scales into smooth and fluctuating parts using the filter.
The upwinding acts a dissipation model for the resolved
part of the fluctuation (i.e. models the unresolved effects).

Regarding what is really physical is a harder question.
Here, the stabilization was carried out so that the diffusive
error becomes of "higher order" i.e. acts as a high order filter
which works technically and stabilizes enough. As the paper showed, the SSD scheme
is quite non-dissipative and it is even less dissipative than UW3.

Any further comments? It is interesting to here the opinions here!

-Ville

[FMDenaro]

Originally Posted by ville

Hi,
I'm not doing deconvolution as we are limited to the grid scales i.e. k<=kNyq.
So the argument is that we are not resolving well the waves
at the grid scale ~2*dx. Thereby, it makes sense to
split the resolved scales into smooth and fluctuating parts using the filter.
The upwinding acts a dissipation model for the resolved
part of the fluctuation (i.e. models the unresolved effects).

Regarding what is really physical is a harder question.
Here, the stabilization was carried out so that the diffusive
error becomes of "higher order" i.e. acts as a high order filter
which works technically and stabilizes enough. As the paper showed, the SSD scheme
is quite non-dissipative and it is even less dissipative than UW3.

Any further comments? It is interesting to here the opinions here!

-Ville

Maybe I have not understood ... it seems to me that you are using the classical definition of a differential filter (top-hat), that is you integrate the Taylor expansion around a volume that if (and only if) is centered around x
leads to
u_bar(x)=u(x) + (h^2/24) Lap(u(x)) +.... = (I + A^2 + ... )u(x)

Then, you truncate and defines
u'(x) = - (h^2/24) Lap(u(x))

This is a well known regularization procedure (Leray) that is a deconvolution procedure since you are formally inverting the truncated top-hat operator (I + A^2 + ... )^-1 = (I-A^2 + ....).
This is if I have well understood....
I have many references on such procedure

[ville]

Hi,
Seems like we are not talking about the same thing although I agree with
you that a Laplacian-term is indeed one of the first terms in the expansion
when inverting a filter.. I'd say that using a Laplacian filter in CFD and OpenFOAM
is quite a standard way of separating scales. I prefer the formalism
where the Laplacian filter is defined using (1/kcut)^2 pre-factor because
it relates the filter with the wave length very clearly.
See the cited reference in Phys.Fluids by Stoltz et al. for a similar approach.
One could use some other pre-factor(cutoff) equally well and even argue
a better choice for sure. Here it does not really matter too much
because the truncation error is eventually controllable with the blending factor.
So as long as one separates the scales in a meaningful way one can
control the dissipation with alpha.

Are you familiar with so called implicit LES? I have mostly in mind ILES
when writing the paper. NS-eqs are of form NS(u,p) = 0 and the normal LES
filtered eqs of for NS(uf,pf) = SGSterms, where uf and pf represent the resolved numerical
solution. In ILES the SGSterms are modeled using a truncation error of the convection term.
Here I showed how the diffusive error can be made Odiff(dx) -> Odiff(dx^3)
using this numerical trick.

Deconvolution approaches belong to their own category. I don't really prefer
them because one tries to invert a filter using the van Cittert procedure and
"gains" SGS terms represented on the grid scale. So what are the SGS
terms represented on the grid scale is another
question. Anyways, even in deconvolution methods damping terms are added
and I think there is not very much difference with the SSD scheme.

-V

[FMDenaro]

Ok, I think that the issues could be seen in a unified way...
For example, see Eq.s (6) to (9) in
http://onlinelibrary.wiley.com/doi/1...O;2-9/abstract

as well as Eq.s (6) to (9) in
http://onlinelibrary.wiley.com/doi/1...d.278/abstract

I think that is possible to give a common framework
Regards
Filippo

[ville]

Thanks for the links, I'll have a look at them on monday as I don't have access to
them right now.
Best,
Ville

[pankajhrajput]

Originally Posted by ville

DxDyDz = (filterScope*filterScope)*DxDyDz/(4.0*pi*pi);

Hi Ville

I am trying to code your scheme into a compressible OpenFOAM solver which I intend to use for a low dissipation LES for a transonic jet (something similar to your work "Large-eddy simulation of highly underexpanded transient gas jets" Phys. of Fluids). However, I am not able to understand how and where you specified the 'filterScope' term. Is it a parameter that you have to specify beforehand or is it computed during run-time?

Thanks
Pankaj