CFD Online Discussion Forums - Looping over all faces within a boundary condition (turbulent inlet development)

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Looping over all faces within a boundary condition (turbulent inlet development) (https://www.cfd-online.com/Forums/openfoam-programming-development/152019-looping-over-all-faces-within-boundary-condition-turbulent-inlet-development.html)

Looping over all faces within a boundary condition (turbulent inlet development)

Background

I have coded up the work of Billson 2003 (turbulent inlet) within openFoam with a lot of help from the material found on this webpage:

http://www.tfd.chalmers.se/~lada/pro.../proright.html

However as some of you may know when a turbulent inlet is applied to an incompressible simulation the pressure field is complete nonsense (this is also a problem with openFoam's supplied turbulent inlet when used with incompressible solvers). As it is the pressure field I am interested in I am now trying to implement a mass flux correction to solve the problem. For those interested the paper I have used as a reference is Kim 2013 "divergence free turbulence inflow conditions for large eddy simulations with incompressible flow solvers".

Please be aware the code below is not mathematically correct, I just want to get it running first before dealing with the finer mathematical points.

What I have been doing is first summing up the normal velocity components using the following:

Code:

forAll(patchField , facei)

        {

                normaliseSum(patchField[facei],uNormalFaceSum);



                faceCounter = facei;

        }

The function being called is:

Code:

template<class Type>

void turbStochInletFvPatchField<Type>::normaliseSum(const vector& u,scalar& uNormalFaceSum)

{

//This is temporary until I code in a general case to find the normal component of velocity to the turbulent inlet:

        scalar uNormal = u[0]; 





//Then sum the normal components, passing by reference to the function ensures the scalar uNormalFaceSum is updated:

        uNormalFaceSum += uNormal;

}

Then I had hoped to use the summed turbulent field and the bulk flow field that openFoam is expecting to normalise the turbulent field, thus resolving the mass imbalance which is causing the pressure fluctuations for an incompressible solver.

Code:

scalar uncorrectedBulkFlowField = uNormalFaceSum / faceCounter;



scalar massFluxCorrector = meanU / uncorrectedBulkFlowField;



patchField = massFluxCorrector * patchField;

However this doesn't work because I am looping over a patchField. Basically if you change any variable in the loop it appears to change the variable type so it has all the properties of the patchField. This change of variable is then passed on to variables outside the loop meaning the scalar correction factor I was hoping to find (by dividing the expected mass flux from the expected bulk field by the actual mass flux) ends up becoming a scalar / patchField hybrid, annoyingly this can be added, subtracted or multiplied with scalars (and itself) but not divided which is what I need to do.

Question
Does anyone know how to loop over patch faces without looping over patchField whilst within a boundary condition? I don't want to have to call the name of the patch I am working on either as that will reduce the portability of the code. I may also later need help on how to then extract the velocity field from the patch.

This will also help another problem I have, because of the patchField loop the function normaliseSum also has to be defined for:

Code:

        void normaliseSum(scalar,scalar&);

        void normaliseSum(sphericalTensor,scalar&);

        void normaliseSum(symmTensor,scalar&);

        void normaliseSum(tensor,scalar&);

Despite the fact I am only using:

Code:

void normaliseSum(const vector&,scalar&);

Things I have tried to get the code working
Changing the function to a scalar function then returning a scalar value, doesn't work the patchField loop still changes the variables type.

Trying to redefine the variables outside the loop by renaming them as something else. Doesn't work the hybrid variable is passed to the new one.

Defining the variables as scalars within the loop, won't compile.

Using a for loop, rather than forAll, for loop doesn't work.

Thank you for any help.

Just in case anyone else stumbles across the problem.

I have managed to semi solve the problem. I was utterly wrong about the cause though. What was happening is that using the ForAll loop spreads the scalar variable across multiple processors, so dividing it doesn't work (other simple operations were fine, but not complex ones like log etc)

I added:

Code:

reduce(uNormalFaceSum, sumOp<scalar>());

after the loop.

This has created a new problem in that it now sums the scalar value over each processor on the inlet face (surprise surprise, I told it to).For example if I decompose it onto four processors, of which the inlet is over two processor domains instead of an answer of say 10, it doubles and gives me the answer 20.

But I'm sure I can sort that out tomorrow morning. I will update once it is done. Although if anyone else wants to comment (on ways around the inefficiency of looping over the face on multiple processors) feel free to.

Hi Mat,

Did you manage to sort out your boundary condition? I'm having a problem with the superfluous pressure waves coming in after I've used turbulentInlet BC and need the pressure data.

Sorry for the slightly slow reply, I forgot my login details...

The problem is with the underlying physics of the turbulent inlet. Under a laminar flow condition at the inlet the velocity flow field is divergence free i.e.

div(V) = 0

However, this condition does not hold for the OpenFOAM turbulent inlet, the velocity field across the inlet is not divergence free. This means, in simplified terms that the PISO/SIMPLE /etc algorithm has to produce extreme pressure fluctuations to support the non-divergence velocity field at the inlet, which are then transported downstream to the rest of the domain.

To get round this you'll have to use a different turbulent inlet which is divergence free (or closer to being divergence free). For larger inlet sizes with small scale turbulence fluctuations the turbulent inlet of the Rostock OpenFOAM extension can be used. I believe it was coded based on the work of Kornev 2007, Method of random spots for generation of synthetic inhomogeneuos turbulent fields with prescribed autocorrelation functions.

http://www.lemos.uni-rostock.de/en/d.../cfd-software/

I extracted the vortons library and turbulent inlet and compiled and used them separately, rather than installing all parts of the Rostock extension, mainly because I had no use for the other parts and I also altered the Rostock turbulent inlet and didn't want to recompile the entire Rostock library for one change in the inlet code.

For smaller inlets with larger turbulent fluctuations I applied a mass flux correction to the Rostock inlet. See Poletto 2013, A new divergence free synthetic eddy method for the reproduction of inlet flow conditions for les and Kim 2013 Divergence free turbulence inflow conditions for large eddy simulations with incompressible flow solvers, for examples. I am not sure who came up with the idea first.

This is required as the Rostock inlet does not produce a totally divergence free turbulence field for all turbulent fluctuations over all inlet sizes. It works fine for larger inlets with small or large turbulent fluctuations but the pressure issue returns for very large fluctuations (turbulence intensity greater than 10) for "small" inlets. I had issues once the inlet became smaller than 2m.

In addition the closer to the inlet you want to collect pressure data the greater the pressure fluctuations become. Again this issue is worse for smaller inlet sizes and larger turbulent fluctuations but if you are collecting data very close to the inlet you might want to apply a mass flux correction regardless.

The reason for this issue is the Rostock inlet does not produce a mathematically totally divergence free turbulence field for all turbulence fluctuations for all inlet sizes. If you did want to pursue coding an inlet which was divergence free for all fluctuations/ inlet sizes I would look at the work of Yu 2014, A fully divergence free method for generation of inhomogeneous and anisotropic turbulence with large spatial variation.

I hope that helps and good luck in solving your problem!

Quote:

Originally Posted by Sevex (Post 543167)

Does anyone know how to loop over patch faces without looping over patchField whilst within a boundary condition? I don't want to have to call the name of the patch I am working on either as that will reduce the portability of the code.

Hi there,

I am not sure if this is what you are looking for but you can loop over the patch faces using the following code within the BC:

Code:

const fvMesh& mesh     = patch().boundaryMesh().mesh();

const faceList& faces  = mesh.faces();

const label& startFace = patch().patch().start();

const label    endFace = startFace + size() - 1 ;



for(i=startFace;i<=endFace;i++)

{

        // your code

}

I have implemented it successfully in one of my custom BCs. I don't know if it is the most elegant way but it works.

Quote:

Originally Posted by Sevex (Post 543167)

I may also later need help on how to then extract the velocity field from the patch.

For this you can use the code that is implemented in many BCs:

Code:

const fvPatchVectorField& Up =

        patch().lookupPatchField<volVectorField, vector>("U");

All the best!

Matthew, thank you for such a descriptive answer. I read it in another forum later and the new OF (ver 1606) seems to have implemented a boundary condition that solves this issue. I'm also trying to workout my way for a divergence free boundary condition too.