Implementing a new LES Model in OpenFoam
 

Hello Guys,

I have a tricky problem with the implementation of a new LES Model to OpenFoam. The new LES-Model is more a vegetation model and consists of two parts. First an additional drag term for ne NS-equation. Second a transport equation for k (unresolved turbulent kinetic energy).
So I set up a ne solver, named woodles and added to the Ueqn the following (red):

fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ sgsModel->divDevBeff(U)
+ fvm::Sp(cd*mag(fvc::Sp(lad,U)),U)
);
lad is a volScalarField and read in as IODictionary from constant directory. cd is a scalar read in from the transportProperties dictionary.
Than I added a new class for my LES Model derived from the LESModel class.
Made my changes for k and nuSgs ran into a little problem.
I need cd and lad in my LESModel too. But I dont want to read in my lad dictonary and the transportProperties dictionary twice (reducing file IO). Is there a way to pipe the values of cd and lad to my SGS Model without changes in the GenEddyVisc and LESModel class?

kind regards,
Fabian

I think you would want to do it the other way around. The constructor for your new SGS model class reads in the cd and lad values from files; you then provide additional access functions in the new class giving access to these fields. When the code calls the model generically (sgsModel->divDefBeff(U)) it behaves like an LESModel since it is derived from that and implements the interface, but you can cast it and access the new functions that way (a bit messy though).

Alternatively - why not subsume your term fvm::Sp(cd*mag(fvc::Sp(lad,U)),)) into divDevBeff(U)? Just overload divDevBeff in the new class to return the whole expression.

Gavin

Thx for your reply. I tested it before to read in the fields in my new model, but for a reason I forgot, this didn't work. So I decided to put it into my solver. Overloading the divDevBeff function isn't the best way in my opinion. I don't want to "hack" my foam. The idea to provide a new access function for the force term seems to be a very good way. But I am not sure how to enable my solver to access the new function. Using sgsModel->newfunction(U) will not work for my understanding because the solver do not know during compilation which model I want to use. So I have to access it directly: newModel::newfunction(U) right? This looks like a hack to ...

kind regards,
fabian

Disclaimer; this is a bit fiddly and is stretching the limit of my knowledge of C++ programming, certainly without spending 1/2 an hour reading Stroustrup. One of the more experienced programmers will undoubtably correct me on this though!

I _think_ that you would need to do some sort of cast operation to ensure that the compiler knows that you are treating sgsModel as a newModel. However this is ugly and will break if you try and run it with any other SGS model. Alternatively you could probably rewrite the whole solver to remove the run time selection, in which case sgsModel can be made a newModel throughout. However I really don't see why you shouldn't overload divDevBeff. What you are trying to achieve is a modification of the SGS term (I think) and so it should be in there. It doesn't really "hack" the code - this is the correct way of providing this (IMHO). Also it shouldn't be anything more involved than providing a wrapper around the existing function which calls the underlying function, adds on your new term, and spits out the whole matrix.

Gavin

Allright, after thinking about your suggestion, I have to agree. The only way to tell the solver about some new functions is to add something to LESModel class or write my own base class for vegetation models. This is to time expensive and I would prefer to expand the divDevBeff function. But after a first try, a problem arise. I would like to write out the term fvm::Sp(cd*mag(fvc::Sp(lad,U)),U) (is stands for a drag force) but I dont have an idea how to achieve this. And I remember, this was the problem last time too. In my solver I put this:

volVectorField Fveg
(
IOobject
(
"Fveg",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
cd*lad*mag(U)*U
);
runTime.write();

after the end of the PISO loop and this works. But at the moment I dont have an idea how this will work in my LESModel.

I am still not very happy with the solution to put it into the diffDevBeff function. For me it is still a hack and it restricts further changes. I tried today to construct a new LESModel class and derive my own model from this new LESModel class. But there came up some problems. I think this is a little bit connected to my less knowlegde of C++, Fortran is more straight away but I will try it again. May be somebody tried this before and give me some hints? Or there is another way to enable my solver to access some additional functions in my own LES model? Nobody has an good solution?


Fabian

So I will answer my self :-D I copied the LESModel class to my OpenFoam user dir and modified this class. Now it provides an additional drag function neccesary for vegetation models. My vegetation model is now derived from my new extended base class and is added to my own solver. So everything works very good and there is no hack (in my opinion).

Fabian