[MartinaFormichetti]

Hi everyone!

I am trying to implement a dUdy gradient in a wall function, I want this calculated by sampling U2(x) two vertical rows of cells away from the wall at y2 and U1(x) at the row of cells next to the wall at y1. I am not an expert in coding but I read on https://jibranhaider.com/blog/mesh-i...n-in-openfoam/ that I should use these mesh connectivities to access mesh information relative to the boundary patch. I just don't understand the usage of those or their definition very well. Based on other wall functions I tried the following:

const label patchi = patch().index();
    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
    (
        IOobject::groupName
        (
            turbulenceModel::propertiesName,
            internalField().group()
        )
    );
    // Access Ux and y
    const volVectorField& U = turbModel.U();
    const volScalarField& Ux = U.component(0); // Assuming Ux is along x-direction
    
    const fvPatchScalarField& Uw = Ux.boundaryField()[patchi];
    const scalarField U1(Uw.patchInternalField());
    const scalarField U2(U1.patchInternalField());

    const scalarField& y = turbModel.y()[patchi];
    const scalarField y1(y.patchInternalField());
    const scalarField y2(y1.patchInternalField());
    
    const scalarField& dUdy = mag((U2-U1)/(y2-y1));

I understand why it doesn't work because the scalarField class doesn't have a member named patchInternalField() (this is also the error I got), it is just to make it easier to understand what I want to do.

I also tried some variations of the following code, all without success because I don't understand how these connectivities and indices work:

 const label patchi = patch().index(); // Index of the wall patch

    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>(
        IOobject::groupName(
            turbulenceModel::propertiesName,
            internalField().group()
        )
    );

    // Access Ux and y
    const volVectorField& U = turbModel.U();
    const volScalarField& Ux = U.component(0); // Assuming Ux is along x-direction
    const fvMesh& mesh = U.mesh();
    const vectorField& cellCentres = mesh.cellCentres();
    const scalarField& y = cellCentres.component(1); // Assuming y-coordinate is along y-direction

    // Calculate U1 and U2
    const scalarField U1 = Ux[patchi + 1]; // Ux at the wall patch (patchi)
    const scalarField U2 = Ux[patchi + 2]; // Ux one cell centroid above the wall patch

    // Calculate y1 and y2
    const scalarField y1 = y[patch + 1]; // y at the wall patch (patchi)
    const scalarField y2 = y[patchi + 2]; // y one cell centroid above the wall patch

    // Calculate the gradient component d(Ux)/dy
    const scalarField dUdy = mag((U2 - U1) / (y2 - y1));

Any help, tips or sources to read more about this and understand it better would be highly appreciated.
Thank you in advance to all who will reply!

Best,
Martina

[Tobermory]

Dear Martina

First of all, I assume that this is a low Re number wall function, otherwise the wall shear stress is not directly connected to the velocity gradient at the wall (if this is gibberish to you, then you need to research nutkWallFunction and wall functions in general).

Now for your question ... it's a far from trivial thing you are trying to do. Let me muddy the waters a little: imagine a patch face on the wall; you can find the adjacent cell fairly easily, and if you know which face is on the top of this face, you could use that to find the cell 2 away from the wall. Sounds good, eh? But the wall adjacent cell probably has 6 faces (or maybe more if its a poly, or maybe less if its a tet!), so how are you going to choose the right face? If it's a tet or a pyramid, what then?

Ok - simplify - assume you have a hex-only mesh, then you could cycle through each face of the wall adjacent cell and then find the face that has the largest wall-normal distance ... and then use that to find the cell 2 away from the wall. Phew - that will take some careful coding.

Now what happens at a corner?

Overall, I have possibly not been terribly helpful here, and please don't let me put you off, but I just want to make sure that you are aware of the magnitude of the journey you are embarking on ... good luck.

[MartinaFormichetti]

Originally Posted by Tobermory

Dear Martina

First of all, I assume that this is a low Re number wall function, otherwise the wall shear stress is not directly connected to the velocity gradient at the wall (if this is gibberish to you, then you need to research nutkWallFunction and wall functions in general).

Now for your question ... it's a far from trivial thing you are trying to do. Let me muddy the waters a little: imagine a patch face on the wall; you can find the adjacent cell fairly easily, and if you know which face is on the top of this face, you could use that to find the cell 2 away from the wall. Sounds good, eh? But the wall adjacent cell probably has 6 faces (or maybe more if its a poly, or maybe less if its a tet!), so how are you going to choose the right face? If it's a tet or a pyramid, what then?

Ok - simplify - assume you have a hex-only mesh, then you could cycle through each face of the wall adjacent cell and then find the face that has the largest wall-normal distance ... and then use that to find the cell 2 away from the wall. Phew - that will take some careful coding.

Now what happens at a corner?

Overall, I have possibly not been terribly helpful here, and please don't let me put you off, but I just want to make sure that you are aware of the magnitude of the journey you are embarking on ... good luck.

Hi Tobermory,

I am quite familiar with the physics and mathematics of the problem I am dealing with.

I also know that it is a challenge to do what I am trying to do with computing the gradient using data points from two rows of cells above the surface and the cells just above the surface (since I have been trying for the past week or so), which is the reason why I posted on this forum since I cannot find any material that I fully understand on the mesh connectivities or the way openfoam uses indexes to distinguish between cells or loop over them.

Apart from your comment explaining why it is difficult, do you know of any sources I could look into to get more info on this and do this?

Thank you for your help.
Martina

[Tobermory]

Happy to hear! You'd be surprised how many innocent posts are raised without the poster realising what they are embarking on!

Since what you are trying to do is at such a fundamental level, I doubt that you will find a reference that will hold your hand and walk you through it. But it's pretty simple I think - it will still require some leg work from your side, but you will learn a lot in the process. The Doxygen pages are invaluable.

As a starter, to find out how you can access the wall adjacent cell look in epsilonWallFunctionFvPatchScalarField::calculate (https://cpp.openfoam.org/v8/epsilonW...8C_source.html), where you will see an example of this.

Now to get the next cell out, you need to research how to address the faces of a given cell - start looking at fvMesh, and drill down the derived classes and you'll find member function faces(). That might be useful. Once you have the faces, you are almost there, albeit as I discussed you need to decide which face and cater for edge effects like corners.

Good luck!

[MartinaFormichetti]

Originally Posted by Tobermory

Happy to hear! You'd be surprised how many innocent posts are raised without the poster realising what they are embarking on!

Since what you are trying to do is at such a fundamental level, I doubt that you will find a reference that will hold your hand and walk you through it. But it's pretty simple I think - it will still require some leg work from your side, but you will learn a lot in the process. The Doxygen pages are invaluable.

As a starter, to find out how you can access the wall adjacent cell look in epsilonWallFunctionFvPatchScalarField::calculate (https://cpp.openfoam.org/v8/epsilonW...8C_source.html), where you will see an example of this.

Now to get the next cell out, you need to research how to address the faces of a given cell - start looking at fvMesh, and drill down the derived classes and you'll find member function faces(). That might be useful. Once you have the faces, you are almost there, albeit as I discussed you need to decide which face and cater for edge effects like corners.

Good luck!

I see, I have been looking for a long time at the nutkWallFunction - https://www.openfoam.com/documentati...8C_source.html - I think the way they get the data at the neighbouring cells is quite similar to the one you suggested with data at the patch given by .boundaryField()[patchi] and the neighbouring cells given by .internalField().

My problem so far has been, as you also mentioned, getting the data at the cells above this first row of cells. I looked it up online and I only found this page - https://jibranhaider.com/blog/mesh-i...n-in-openfoam/ - mentioning the mesh connectivities/mesh coordinates but it doesn't say much. I wasn't looking for a guide to do this as I assume not many people have attempted this before but at least a manual page or something that gives some more information on these. Do you know if I can look this up in more detail anywhere?

[Tobermory]

I wasn't looking for a guide to do this as I assume not many people have attempted this before but at least a manual page or something that gives some more information on these. Do you know if I can look this up in more detail anywhere?

Nope - I'm afraid that there's little official document of the code basics - the Developers assume that you will either spend the time getting your head around the whole code (which is the best way, if you are a long term user), or you will pay them to develop a bespoke tool/solver. Fair enough I guess - it is freeware, and they need to make a living.

There are some other sources online, training material and papers/Masters theses, that you'll find from Google searches - Chalmers, Wolf Dynamics etc. But I've already told you where to look in the code ...