[pablitobass]

Dear All,

I have been working for a while on the implementation of a new FLux Limiter in OpneFOAM, and after some attempts it seems that I found they way to do that. I have compiled it and it seems that it works properly (at least it allowed me to get reasonable results). However, I have a doubt, tthus I would ask if someone could be so kind to help me.

I know that a flux limiter can be seen as a function that allows a "blended" scheme to behave like an UD scheme or like a HR scheme in different region of the flow field depending on the particular conditions in the flow field itself.
In principle (I guess), the HR scheme could be whatever you want, but in my case I need to use the QUICK scheme of Leonard with a particular limiter.

Is it possible to specify this scheme in the field "interpolationScheme" into the "fvScheme" of the case? (I tried but it seems that it does not work. Probably I used wrong settings)

Does anyone know what is the default HR scheme? Maybe the linear?

Thank you very much in avantage!

Paolo

[Tushar@cfd]

I think GAMMA scheme is mostly used for OpenFOAM work

-
Best Luck!

[pablitobass]

Hi Tushar,

Thanks a lot for the reply!

In practice, my problem is the following. The flux at the cell face can be calculated as:

(phi)f = (phi)_UD + psi(r)*[(phi)_HR + (phi)_UD]

where (phi)_HR is the flux at the face "f" evaluated whit a Higher Order scheme, which in my case should be the QUICK scheme, and psi(r) is the flux limiter that I have added on my version of OpenFOAM.

by considering what you said, I am not allowed to use a QUICK scheme for the interpolation at the cell face? However I am a bit confised. In fact, I had have a read to to an article in OpenFOAM Wiki, and I found the this:

The surfaceInterpolation classes implemented in the finite volume library perform interpolation from volume fields to face fields, a critical calculation in the discretization process for the finite volume method, particularly for co-located meshes, as is used by OpenFOAM. OpenFOAM implements dozens of different schemes for this operation. This kind of interpolation is denoted:


Notation used by Marupio in this wiki


1 Schemes The schemes are selected in the fvSchemes file, and loaded using runTime selection. There are more than fifty schemes available in OpenFOAM-extend, including:

• QUICK - quadratic upwind interpolation;
• linear - central differencing scheme;
• upwind - upwind differencing scheme; and
• skewCorrected - skew upwind scheme.

and seems that would be possible to select the QUICK scheme, but I do not understand how to do that.

Does anyone have any suggestion?


Thanks a lot!


Paolo

[Tushar@cfd]

Refer the following

http://www.openfoam.org/docs/user/fvSchemes.php

check out section 4.4.5, you can apply the above said scheme as:

div(phi,U)        Gauss QUICK;

I think this will work

For more information refer the cfd-online post (I guess you know this):

http://www.cfd-online.com/Forums/ope...nder-hood.html

-
Best Luck!

[pablitobass]

Hi Tushar,

Thanks again for helping me! I am afraid that when you specify "Gauss QUICK" on the fvSchemes dictionary, you are actually specifying the QUICK limiter. I have been looking for a long time the QUICK scheme inside OpenFOAM and the only files I found about it are actually inside the following foleder:

src/fiiteVolume/interpolation/surfaceInterpolation/limitedSchemes/QUICK

on which there are two files: QUICK.H and QUICK.C

I am pretty sure that those files define the QUICK flux limiter and not the quadratic interpolation scheme. In fact as you can see here:

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.

OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.

Class
Foam::QUICKLimiter

Description
Class with limiter function which returns the limiter for the
quadratic-upwind differencing scheme.

Note that the weighting factors are not bounded between upwind and
central-differencing, some downwind contribution is possible although
the interpolate is limited to be between the upwind and downwind cell
values.

Used in conjunction with the template class LimitedScheme.

SourceFiles
QUICK.C

\*---------------------------------------------------------------------------*/

#ifndef QUICK_H
#define QUICK_H

#include "vector.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
Class QUICKLimiter Declaration
\*---------------------------------------------------------------------------*/

template<class LimiterFunc>
class QUICKLimiter
:
public LimiterFunc
{

public:

QUICKLimiter(Istream&)
{}

scalar limiter
(
const scalar cdWeight,
const scalar faceFlux,
const typename LimiterFunc:hiType& phiP,
const typename LimiterFunc:hiType& phiN,
const typename LimiterFunc::gradPhiType& gradcP,
const typename LimiterFunc::gradPhiType& gradcN,
const vector& d
) const
{
scalar phiCD = cdWeight*phiP + (1 - cdWeight)*phiN;

scalar phiU, phif;

if (faceFlux > 0)
{
phiU = phiP;
phif = 0.5*(phiCD + phiP + (1 - cdWeight)*(d & gradcP));
}
else
{
phiU = phiN;
phif = 0.5*(phiCD + phiN - cdWeight*(d & gradcN));
}

// Calculate the effective limiter for the QUICK interpolation
scalar QLimiter = (phif - phiU)/stabilise(phiCD - phiU, SMALL);

// Limit the limiter between upwind and downwind
return max(min(QLimiter, 2), 0);
}
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************** *********************** //

See also here for reference:

http://en.wikipedia.org/wiki/Flux_limiter

My goal is to specify a limiter in this way:

PHP Code:

div(phi,U)        Gauss My_Limiter; 


and use the quadratic interpolation for the face flux as a Higher Order scheme.

Thanks a lot for your time and your kindness!

Best Wishes,
Paolo

[crazzy.pirate43]

Dear pablitobass,

If you please, I'm trying to use the blended scheme but every time iput into the interpolation scheme I get eror so is there is any specific method to write it? and what is the relation between it and the divschemes?

Thanks in advance

[pablitobass]

Hi crazzy.pirate,

I am afraid I did not really undestand your question. Could you be a bit more precise?

However, to my knowledge, when you define the entries in the divergence scheme you have two choices;

1. divergence of a non-advective term,
2. divergence of an advective term

The way they are defined is different. In the first case, if you have any vector or tensor field F, you have to use the following sintax

div(F) Gauss linear;

which means that you are using a Gauss integration and a linear interpolation scheme to evaluate the field F at the cell faces.

In the second case, considering once again the field F, your entries will be

div(phi, F) Gauss <selected_interpolation_scheme>;

In this case you have a variety of options to chose which includes blended and limited schemes.

Hope it helps.

Paolo

[arjun]

Originally Posted by pablitobass

Hi Tushar,

Thanks a lot for the reply!

In practice, my problem is the following. The flux at the cell face can be calculated as:

(phi)f = (phi)_UD + psi(r)*[(phi)_HR + (phi)_UD]

where (phi)_HR is the flux at the face "f" evaluated whit a Higher Order scheme, which in my case should be the QUICK scheme, and psi(r) is the flux limiter that I have added on my version of OpenFOAM.


Paolo

I think this is a typo


(phi)f = (phi)_UD + psi(r)*[(phi)_HR - (phi)_UD]

[pablitobass]

Originally Posted by arjun

I think this is a typo


(phi)f = (phi)_UD + psi(r)*[(phi)_HR - (phi)_UD]

Might be that I have messed up something, arjun, but I can't really see what. If I consider that the limiter is zero

psi(r)=0

the previous formula returns the flux calculated with the upwind scheme. On the contrary, when

psi(r)=1

it returns the flux calculated just with the higher order scheme.

If you think there is an error, colud you please point it out? Thanks!

Paolo

[arjun]

Originally Posted by pablitobass

Might be that I have messed up something, arjun, but I can't really see what. If I consider that the limiter is zero

psi(r)=0

the previous formula returns the flux calculated with the upwind scheme. On the contrary, when

psi(r)=1

it returns the flux calculated just with the higher order scheme.

If you think there is an error, colud you please point it out? Thanks!

Paolo

Isn't it supposed to be like that that 0 for upwind and 1 for HR and intermediate values for blend.

When I wrote i though limiter to be this blending. This is why I pointed out.

If you have other behavior in mind that i mistook what you wrote.

[crazzy.pirate43]

Dear Pablitobass,

Thanks a lot for quick reply. The problem that in my case there is no non-advective terms but however, is it right that all the non-advective terms take (Gauss linear) as interpolation scheme.
On the other hand, for the advective terms the selected interpolation scheme wich you refered to has to be the same like that one in the interpolation scheme itself or not?

Thanks for your help