[CFDanielGER]

Hello everybody,

I am currently trying to implement a new hybrid turbulence model in OpenFOAM (An Enhanced Version of DES with Rapid Transition from RANS to LES in Separated Flows) but unfortunately got stuck.

At the moment I have two questions.

1)
For calculating the new length scale for the model I first need to access for each cell the cell vertices and the cell center and build the corresponding vectors. Then I need to do calculations with these vectors and the vorticity vector. I guess I need a for-loop for this. How can I do this in OpenFOAM in an efficient way (accessing the center and vertices of the cell and do the calculations)?

2)
For another quantity I need to build the average over the cell and its closest neighbours. How can I access the closest neighbours of a cell in OpenFOAM? This function should also work in parallel computing (If I have to access a neighbour cell from another subdomain for example).

I really appreciate your help. If you need additional information please ask.

Many thanks for your help!
Daniel

[clapointe]

I think I can help with the 2nd question, at least. You can access cell neighbours with cellCells(). I remember seeing/using something like that long ago... here's a good thread with straightforward usage : How to identify cell neighbours.

As for the first question, I'm not quite sure. I imagine that you could get a list of points associated with each cell in a similar way to the way you get the cell neighbours. A quick search and it looks like you can make a list of points, and use indices obtained from cellPoints to pull the coordinates of specific points from that list : Point labels associated with a cell. The cell coordinates are mesh.C() if I remember correctly.

You might also check out what's done in for the IDDESDelta -- it may provide some clues/be a good reference for your implementation.

Hope this helps,
Caelan

[CFDanielGER]

Thank you very much for your help.

I decided now to first implement an easier version of the new length scale which is based on deltax, deltay and deltaz. I saw from the maxDeltaxyz.C that the maximum delta is already calculated. I also think that it is used in the cubeRootVolDelta.C file.
My question is, do I have to code and calculate the delta in each direction from the beginning or is it also possible to import it?

Many thanks for you help.

Best regards,
CFDaniel

[clapointe]

I imagine that if your delta is sufficiently similar to an existing one, it'd be easiest to adapt existing code to suit your needs. Building off the preexisting code structure will also help ensure that your DES model will be able to access the new delta as well.

Caelan

[CFDanielGER]

What I am currently doing is to adapt the SpalartAllmarasDDES.C file since for calculating the new delta I also need some flow properties and it seemed "easier" to work directly in the turbulence model file than working in a new delta file (I am very unexperienced by coding in OpenFOAM but I got some tipps that working directly in the turbulence model file may be more convenient since I can access easier the flow variables I need).

I am now wondering if I can access the deltas in each direction or if I have to calculate them by myself? For example when I want to access the viscosity or velocity I just write this->nu() or this->U(). Now I am wondering if the same is possible with the deltax/deltay/deltaz.

[clapointe]

In les turb models, the delta is accessed in a similar way :

this->delta()

To start, you could probably assume (at the expense of generality) the use of one delta type. Then its calculation is already coded -- you'd probably just need to access it (like above) and then augment the delta in your les model.

Caelan

[CFDanielGER]

Yes, I also saw this. I think this delta then is dependent on the input you select in the turbulenceProperties file. Since there is only maxDeltaxyz and cubeRootVolDelt available I guess I then have to implement some line of codes to get deltax, deltay and deltaz?

Probably something like this (found in the IDDESDelta.C) but for each direction in x, y and z:

    scalarField& faceToFacenMax = tfaceToFacenMax.ref().primitiveFieldRef();

    const cellList& cells = mesh.cells();
    const vectorField& faceCentres = mesh.faceCentres();

    forAll(cells, celli)
    {
        scalar maxDelta = 0.0;
        const labelList& cFaces = cells[celli];
        const vector nci = n[celli];

        forAll(cFaces, cFacei)
        {
            label facei = cFaces[cFacei];
            const point& fci = faceCentres[facei];

            forAll(cFaces, cFacej)
            {
                label facej = cFaces[cFacej];
                const point& fcj = faceCentres[facej];
                scalar ndfc = nci & (fcj - fci);

                if (ndfc > maxDelta)
                {
                    maxDelta = ndfc;
                }
            }
        }

        faceToFacenMax[celli] = maxDelta;
    }