[CedricVH]

There is a bug in simulations with lagrangian particles (basicKinematicParcels in a kinematicCloud) using wall interactions in OpenFOAM 1.6.x. I have defined in kinematicCloudProperties that particles hitting the wall, should stick to it:

StandardWallInteractionCoeffs
{
    type        stick;
}

In StandardWallInteraction.C, in the correct(...) function it is found that this means that, if a particle hits the wall:

• The particle is kept in the system
• U of the particle is set to 0 ( U = vector::zero; )

This function is called by the hitPatch(...) function in KinematicParcel.C. However, it seems that this function is never called!

I am using particles with a diameter of 9e-6 which is much smaller as my grid size, so it is conform with the assumptions of OpenFOAM. I have done a simulation without flow and only gravity working on the particles (which I set to (0 -1000 0) to speed it up). I injected the particles somewhere and let them fall to a bottom wall. A video of my simulation is seen here: http://nl.tinypic.com/player.php?v=qp2auu&s=5

Instead of hitting the wall, the particles float above it with a distance of 9.4955e-06. This distance is bigger than the particle diameter, so it is not the distance from wall to particle centre. The U is also not becoming (0 0 0), but has a value that is much bigger which means that U = vector::zero is not called.

You can also see in the video that the particles behave strangely when hitting the bottom, the particles are grouping together. This is not a bug in the post processing software as these results are the same in Paraview and TecPlot.

I have attached my case with this post.

[niklas]

which code are you using?

[CedricVH]

I have made a custom solver (pimpleLagrangianFoam), but it happens with every solver that implements kinematicCloud, such as rhoPisoTwinParcelFoam.

[niklas]

If I cant reproduce the result, its gonna be hard to fix it.
Can you send me the code (rhoPisoTwinParcelFoam is new to me, which version of foam is that in?)

niklas dot nordin at nequam dot se

N

[CedricVH]

rhoPisoTwinParcelFoam is in OpenFOAM 1.6.x, but it is not in the src directory. You can find the solver in tutorials/lagrangian/rhoPisoTwinParcelFoam. Before you can use it, you have to compile it with wmake.

In the same directory, you have the case simplifiedSiwek. Open the file constant kinematicCloud1Properties and change the StandardWallInteractionCoeffs from rebound to stick. Then just run rhoPisoTwinParcelFoam on this case.

At the last time step, most particles should be stuck to the bottom wall (Y=0). However, when you check the timestep/lagrangian/kinematicCloud1 folder, you can see that their Y value is a little bit bigger than 0 and that their U vector is not (0 0 0) as you should expect.

[CedricVH]

When doing a git pull today, I saw that the developers had made a change in the surfaceInterpolation.C file and I was already hoping that this would have solved my problem

However, this is not the case (and yes, I have recompiled everything) Is there anyone able to reproduce this problem with the testcase of the rhoPisoTwinParcelFoam solver?

You should first update your OpenFOAM 1.6.x installation with git pull to make this solver work correctly (bugreport of Alberto Passalacqua).

[niklas]

open
src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

and find the line that reads

        dt *= p.trackToFace(p.position() + dt*U_, td);

change it to

	if (mag(U_) > VSMALL)
	{
            dt *= p.trackToFace(p.position() + dt*U_, td);
	}

then lookup the calcVelocity routine and add this (after vector Unew=...)

    vector Unew = Ures.value();

    label facei = this->face();
    if (facei != -1)
    {
        bool internalFace = this->cloud_.internalFace(facei);
	if (!internalFace)
	{
            label patchi = this->patch(facei);
            const polyPatch& patch = mesh.boundaryMesh()[patchi];
	    bool keep = true;
            td.cloud().patchInteraction().correct(patch, facei, keep, Unew);
	}
    }

and type wmake libso in src/lagrangian/intermediate

The first modification is because when you set the particle velocity to zero on the face.
At the next trackToFace it will not move, but it will reach its end position, and this will make
the algorithm believe that it has left the face and no longer belong to it.

the second modification is just to make sure that the forces acting on the particle
does not move it, but keeps it there.

[CedricVH]

Wow, thank you very much Niklas! It works perfectly!

As a super moderator, do you have any contacts with OpenCFD? It would be nice if this bugfix could be implemented in the offcial git version of OpenFOAM 1.6.x!

[niklas]

Great.

Im sure Graham will see this when he gets back from his course and come up with
a better solution

Im not 100% sure how the addition of
td.cloud().patchInteraction().correct(patch, facei, keep, Unew);

in calcVelocity will influence particles that just bounce on the wall
(if it will do the interaction twice)
in which case another solution has to be implemented,
but that requires interactionType knowledge from outside the wallInteraction class,
so I'd wait for an official fix.

[CedricVH]

I don't think that the second part of your code works properly as facei always stays -1. When I run a case and a particle hits the wall, it is not associated with a boundary face.

When I call this code on a cloud:

forAllConstIter(passiveParticleCloud, aCertainCloud, iter)
{
    label partPatch = iter().patch(iter().face());
    Info << "Particle position: " << iter().position() << endl;
    Info << "Particle face: " << iter().face() << endl;
    Info << "Particle patch: " << partPatch << endl;
}

It gives me these results:

Particle position: (0.0227693 4.5e-06 0.05)
Particle face: -1
Particle patch: -1
Particle position: (0.0388337 4.5e-06 0.05)
Particle face: -1
Particle patch: -1
...

The positions are correct (y = 4.5e-6 = radius of particle) and also U = (0 0 0), but the particles are not associated with a boundary. The function iter().onBoundary() also returns false.

I think that these functions work with the particle centre and do not take it's size into account. Maybe a fix in the official git release?

[niklas]

did you do a restart in ascii mode?

ascii mode does not write face info, so a restart removes particles from faces.

[CedricVH]

I don't know what you mean. I just started the case from scratch. The writeFormat in my controlDict file is set to ascii, but I don't know if that has anything to do with it.

[niklas]

ok. so that isnt the problem then.

It's always good to run lagrangian simulations in binary mode.

Then I dont know what's goinf on and I cant really say whats wrong (if anything) since I
dont know where in the code you print statements are.
If you put an info statement inside the
if (facei != -1)
{

does it write something then?

[CedricVH]

The code I posted was in a very simple (and not very correct) functionObject. However, changing the writeFormat in my controlDict file to binary solved the problem!

This is a very strange behaviour and it should be warned for. It's not correct that changing the output format gives different results.

Thank you very much for your help! I really appreciate it!

[niklas]

I agree that it can appear as strange behaviour.

The reason facei (if I remember correctly) is not passed on in ascii mode is
because if you save it in ascii, you have a finite representation of the position.
and this is equivalent of moving it.

It can therefore be abit dangerous to do restarts when you have lagrangian data, since
the position that was written no longer have to be inside the cell (or on a face).

[gmacpherson]

Originally Posted by niklas

Great.

Im sure Graham will see this when he gets back from his course and come up with
a better solution

Im not 100% sure how the addition of
td.cloud().patchInteraction().correct(patch, facei, keep, Unew);

in calcVelocity will influence particles that just bounce on the wall
(if it will do the interaction twice)
in which case another solution has to be implemented,
but that requires interactionType knowledge from outside the wallInteraction class,
so I'd wait for an official fix.

Hi Guys,

Andy has put a fix for this into our internal development version of the particle tracking, but it isn't backwards compatible with current cases. Therefore, this fix will come out with the next release.

Thanks for reporting it, best regards,

Graham

[LaSerpe]

Hi, I'm having similar problems running icoUncoupledKinematicParcelFoam on version 2.3.

the case is a water inpinjment study for a 3d wing: when I run the lagrangian tracking I get a certain number of unactive particles, but when I look at the positions file most of them seem to be stuck in a cell which not own a boundary face.
I've tried to change mesh refinement, I switched off wall function and change integration time step but nothing worked.


Giulio

[wyldckat]

Greetings Giulio,

For some reason I've had your post on my to-do list, to keep in mind to have a look into this as soon as I could... probably because I was curious...

Problem is that the description you've given isn't enough for me to try and replicate this issue with one of OpenFOAM's tutorials, nor am I able to assess if this was a bug that has already been solved in OpenFOAM 2.3.x or not.

Therefore, can you provide more information on how to reproduce the same error you're getting? Or at least an example case where this issue can be reproduced?

Best regards,
Bruno

[LaSerpe]

Dear Bruno,
here I attach 2 picture showing the problem, the first picture is the visualization of the particles after the simulation run, you can see that at the leading edge of the wing the distribution in homogeneous.

The second instead is showing you the problem: I wrote a simple routine which reads which particles is active and which is not, then it takes the label of the cell from the positions file and compare that label with those enlisted in the owner file, starting from line = NInternalFaces, so I will find the label of the boundary cell where my particle is stuck.
This produces an inhomogeneous distribution which is not correct.
So at the end it turns out that some particles are considered unactive but their position is not near the wall (I verified this by reading and jumping from file to file)

This seems to be a problem related to the 3D case, because in our 2D simulation all worked fine

Sorry, I'm not working on this anymore, so I've to search and try to figure out where I saved all the data on my hard drive if you need meshes or dictionaries I used.

Best regards,
Giulio

[wyldckat]

Hi Giulio,

Thanks for the feedback! Well, from your description and images, I can try to create a case to attempt and reproduce the same issue. But if my guess is correct, this is one of those hard cases to reproduce the same error.
My guess is that the particles that were deactivated got welded into their location because of two possible scenarios:

1. The mesh had singularities in the location where the particles got stuck. In other words, there were cell vertexes exactly where the particles were going to hit the wall and therefore got stuck at those vertexes.
2. And/or the particles got stuck because at the time they hit that location, they were heavily packed and surrounded by other particles, leading to the algorithm to mistakenly freeze them there.

If this is the case, it's hard to reproduce, since it:

• Either requires making the mesh be correctly aligned with the flow and how to make the particles hit said location.
• Or requires throwing a lot of particles at the problem... which would be fun

If you can share a test case, it would be great.
If not, I'll try to reproduce this sometime in August.

Best regards,
Bruno