Hiting the wall by particle
 

Hello!
I have a problem with the particles tracking. The problem is to simulate the hit of the particle to the wall (the boundary of the mesh) correctly.
There is solidParticle library in the OpenFOAM. The function trackToFace in solidParticle.C should take into account the size of the particle in the case of collision. But if the cell near the wall is smaller then the size of the particle doesn't work! In this case particle can approach the wall closer than the radius of the particle.

Can you tell me how can I solve this problem?

yup, that is correct.
I think this is a deliberate limitation, because what you want to do is very expensive.

Also, particles can only belong to 1 cell at a time, so how do you treat this?
This should also be considered, but again, extremely expensive to solve.

N

in more detail, what Niklas is mentioning - the lagrangian particle tracking implemented in FOAM is sitting on the assumption, that the particle is smaller than the cell. Otherwise you would need to integrate e.g. the forces on particle around the particle surface as it expands over more cells. New formulation of coupling to the continuous phase would be needed. The continuous phase would not be continuous any more -- force transfer across particle.... You would need to use VOF, or 6DOF or some other treatment...

matej

Niklas and Matej, thank you very much for answer!
Matej, you are right that in the case of big particle we should take into account a lot of forces (for example lift force, that takes into account the circulation of the flow around the particle). I have already implemented this and some other forces in the solidParticle library.
The most problem is in the region close to the wall. We should make the size of the mesh elements near the wall much smaller than in the middle zone to simulate the hydrodynamic correctly, especially in the case of EM governed flow (my case). Thus appears the described problem of the collision with the wall.

The way how to solve this problem that I want to suggest is so. I propose to make a virtual wall in the distance of the radius of the particle from the real wall. Then make the size of the particle zero for the function wallImpactDistance. Thus we will indicate the collision when the center of the particle reach this virtual wall.

Now I can not understand how to implement the check of the collision with this virtual wall in the function trackToFace. If you can give me some advise in this question, I'll be very thankful!

This is a really difficult problem.

Let me outline how I would do it and where I see the challenges. Maybe you have a better idea, in which case I would love to hear it.

Since we are assuming that wall-interaction can occur even though we are not in a cell close to a wall, we have to move our wall-faces inwards according to the wallImpactDistance and check every wall-face before we move the particle to see if we will hit a wall face first.

This operation is pretty straight-forward: we just add 1 extra check to the original algorithm
before we move the particle.

Something like this -
1. original algorithm.
2. check wall faces.
3. move particle.

However...point 2 in the list above is trickier than one would think
Maybe a simple example can serve as an illustration to the problem.

Assume we have 2 walls as below, that we move inwards and
the particle wants to move from A to B.

All faces on wall 1 will intersect with the trajectory A-B, but none will actually be hit.
Hence, when you check for wall-interaction you have to also check if the intersection-point
is within the wall-face.
And, since we have to take into account the size of the particle and move the wall-faces, we
are creating small holes in the mesh through which the particle can escape.

So in short, wall-interaction has to take into account if the interacting point is within the face (very expensive check) and it's also important when you move the wall-faces that you
dont create holes in the mesh. This implies that you cant just move the faces, you have to move the points that constitutes the faces.

I hope it is clear now why I find this a bit difficult.

N

Query
 

Hey Guys,

May be you all can help me, I have been trying to play with Lagrangian particle solver, I had two questions in this regard

a) Is there some restriction on this solver, when we run it in parallel.

b) I have even tried running it on my desktop, every thing keeps on going well for certain time steps but eventually it hangs with message on the screen as follows

Time = 0.41

Courant Number mean: 0.0123124 max: 0.234255
Moving Particles
53 Particles moved. 3 walls hit. 0 particles left the model.
DILUPBiCG: Solving for Ux, Initial residual = 0.000270292, Final residual = 8.83653e-08, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.000453004, Final residual = 9.93462e-07, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.00126948, Final residual = 1.39543e-06, No Iterations 1
DICPCG: Solving for p, Initial residual = 0.0139312, Final residual = 0.00119803, No Iterations 3
DICPCG: Solving for p, Initial residual = 0.00276161, Final residual = 0.00026989, No Iterations 56
time step continuity errors : sum local = 6.67705e-08, global = -5.95145e-09, cumulative = 8.18682e-07
DICPCG: Solving for p, Initial residual = 0.00173553, Final residual = 0.000171299, No Iterations 8
DICPCG: Solving for p, Initial residual = 0.00126418, Final residual = 0.000115941, No Iterations 8
time step continuity errors : sum local = 2.8673e-08, global = -6.02904e-09, cumulative = 8.12653e-07
ExecutionTime = 597.83 s ClockTime = 598 s

Time = 0.42

Courant Number mean: 0.0123111 max: 0.23442
Moving Particles


It just does't proceed beyond this
All other residuals seems to be under control and courant number is well below zero. I am also attaching the file cloudProperties, Any help is enormously appreciated.

Thanks,
Sumit

Niklas,
you are right, there are a lot of problems in your example of 2 walls and obtuse angle between them. Let us initially consider simpler case. Let we have some simple geometry, for example cylinder (there are no obtuse angles). In this case we can build the virtual wall in the distance of particle radius from the real wall.
You propose to apply some function to indicate heating this virtual wall before the move function. But now algorithm check heating with the wall inside of move function. There are trackToFace function that indicate the collision and move the particle only to the face of cell or wall. Thus we should implement this new check into the trackToFace function, but not outside of the move function, as you suggest.
Am I right?

Sumit,
where are some problems with parallel calculation when you are working with the big cloud (huge number of particles). I have not solve this problem yet.

yes, you are right, but that's what I ment.
it will all be in the trackToFace routine. It has to be.

trackToFace checks all faces in the cell to see if the particle can move to its end position
or if it should move the particle to the intersecting face.

What I meant was that before we move it (which is actually done inside the trackToFace routine) we have to add an additional check for all, and I do mean all, wall faces before we can go ahead and move it.

And I really dont recommend that you implement something that will work only on a special type
of geometry.

Niklas,
I am very thankful for your help and ideas!
But I want to ask you some more technical questions.

There is such code in Particle.C:
00307 // change cell
00308 if (internalFace) // Internal face
00309 {
...
00326 }
00327 else
00328 {
00329 ParticleType& p = static_cast<ParticleType&>(*this);
...
00341 if (!p.hitPatch(patch, td, patchi))
00342 {
...
00371 else if (isA<wallPolyPatch>(patch))
00372 {
00373 p.hitWallPatch
00374 (
00375 static_cast<const wallPolyPatch&>(patch), td
00376 );
00377 }
...
00382 }
00383 }

Here is indicated the collision with the wall and applied the hitWallPatch method. But I cannot understand how to add the new virtual wall to the list of faces and to put flag that it is the wall, not cell face.
Can you help me in this question?

As a next step I'll think how to patch up the possible holes in the virtual wall.

That part of the code is just responsible for the wall-treatment and will remain unchanged.

What you want to do is to edit the Particle.C-file and the trackToFace-routine.
just stick an include-file there and add the necessary checks and the rest should work,
since boundary-faces are moved according wallImpactDistance.
Just remove the check in the lambda-routine that wallImpactDistance is shorter than cell-size also of course.

Hitting a face is signalled with facei_ not being equal to -1.
So when you have moved the particle to the virtual wall, setting facei_ accordingly,
the tracking will think it is on a wall and apply the boundary treatment.

Niklas,
I understand, that in myWallCheck.H I should check the wallImpactDistance in direction of normal vector of face for all faces. If this wallImpactDistance for the faces[j] is less than d_/2.0 then we should redefine this faces[j] as a boundary face. And all other algorithm will work.
But
1) how can I get the normal vector of faces[j]?
2) how can I mark this faces[j] as a boundary face?

Thank you a lot, Niklas!

P.S. I think I have found answer on the 1st question:
mesh.Sf()[facei_]/mesh.magSf()[facei_]
Am I right?

correct, and if the face is on a boundary the normal is pointing outwards.

You dont have to mark a face as a boundary face. FOAM knows if the face is a internal face or not.
You just have to loop over every patch and check if its a wall or not, then
you can just grab the faces of an entire patch and loop over those.
not sure of the exact syntax, but to grab them should be something like this
const surfaceVectorField& Sfbi = mesh.Sf().boundaryField()[patchi];
or
const faceList& bcfaces = mesh.cells().boundaryMesh()[patchi] alternatively
mesh.boundaryMesh()[patchi].cells()
(not sure how it is)
Also, remember that cells() returns the indeces to the faces, while faces() returns the indeces to the points.
if it doesnt work try looking in the wallFunctions directory for the turbulenceModels
should be something useful there
If you cant figure it out I will try to find out the exact expression.

Sorry, Niklas, I don't understand your algorithm:
I just don't understand.

What I proposed:
Thus I suggest to deceive the FOAM and tell it that boundary face is not real boundary face, but the face of the cell, when the particle in this cell is closer to the real wall than d_/2.0. If it is possible, of course.

Niklas, would you be so kind to explain me the algorithm you suggest and/or estimate my one?

I am away at a conference at the moment, but when I get back
I will write you a piece of code that I think will make it clear.

N

Thank you very much, Niklas! I'll be waiting for code.
Have a good time on the conference!

On a train back home, lets hope the shaky connection works :)
Here's how you access all the stuff you need.
I'm leaving it up to you to fill out the details, i.e, how you check if
the face-hit is within the face-edges

myWallCheck.H

Niklas, thank you. But I want to ask you to explain some places in your code to understand it:

myWallCheck.H
And general question. I am not sure that I understand your idea correct. So you suggest:
1) to remove the check that the distance between the cell center and face center is larger than wallImpactDistance in lambda function in ParticleI.H
2) to add this code in Particle.C
After this modification will been made, FOAM will check if the particle hit the wall anyway, but not only when the particle is in the boundary cell.

Have I understood your idea correct?

1. I thought that was obvious. It checks if the patch is a wall or not.
2. You will (probably ) need them later when you move the face centre point, to check if the face hit is
withing the face-edges
3. yes.
4. nope. the lambda calculation assumes that the face is infinite, since it only has a face centre point
and a face normal. I would suggest that you decompose the face into triangles using the face-center
and the points on the edges and then check if the face hit is within one of these triangles.

general.
1. yes.
2. yes.

and if the face is hit you need to append the face index to the faces-variable in the trackToFace-routine

I have written the code to check if pHit is within the edges of the face.
But I am not sure that I know how to add this face to faces. Is there any method to add the new face of thirst of all I should increase the size of the faces massive for 1 element and them fill this element this the new face?

faces.append(facei)

...I seriously doubt that check though.
I think you should come up with something else

Niklas, you were right, this check does not work. But I cannot understand why.
S is more than mag(Sf) for all faces!

comparing scalars with == is always a bad idea.

if you are going to do that use mag(scalar1 - scalar2) < eps instead,
or (mag(area1/area2 - 1.0)) < eps to get relative difference (even better here),
but apart from that you are decomposing the area using the point and if a face is slightly warped
(which it always is numerically when you have more than 3 points)
you will get different areas depending on where the point is.
I would not use that area-comparison method, but if you want to get it working on simple
cases where you know the faces arent so bad just increase the tolerance.

Im currently working on a new tracking algorithm and what Im doing there is to decompose the face in
triangles, using its face center and face points, and then check if the face hit is within one of these
triangles, the algorithm for that looks like this

(I just copy/paste'd this so variable names might differ, but you should get the idea I hope)

Your method:
normal1 is outward pointing triangle edge normal.
vector (px - pi) is vector FROM px TO pi.
Thus if px is within the triangle then vector (px - pi) will have the same direction as normal1. Thus (px-pi) & normal1 will be positive (but you've written that it should be negative!)

eehhh....
example:

px = (1 0 0)
pi = (0 0 0)

px - pi = (1 0 0)

how can this vector be from px to pi?

I am really sorry, I am little bit tired(

Niklas, thank you very much. It seems that treated solver works correct. I am really very thankful!

Hi Niklas and Shchepan, I'm also working with particles that near wall have the radius bigger than the cell size and than I have the problem on the rebound.
I've implemented the idea shown in this thread and I've solved almost all the problems on the rebound. I still have some problem on few particles.
I'm thinking in other ways to solve the problem on the hitting of the wall by particles and I have a question for you:
I'm not sure if I have understood properly the logic behind the particle tracking algorithm implemented in OpenFOAM but, in a situation of thetraedric mesh (like the 2D sketch attached below), is it possible that particles with center position in a cell that not have a wall face could not hit the wall properly? I mean although their radius length is smaller than the characteristic cell size?

Thanks for the answer.
Fiorenzo

your problem is the same as shchepan's.
Since the particle only can belong to 1 cell at a time there is a chance the particle can 'escape'.
When your particle can extend into neighboring cells you have to check all boundary faces for wall interaction, otherwise you will make assumptions that can cause problem.

Dear Niklas,

thank you for the answer to Fiorenzo. We are working together on this problem. We would highlight a situation that could cause the algorithm to fail. We have the feeling that this situation has not been considered before. In the attached image, we illustrate two particles hitting the walls starting from different cells but following exactly the same trajectory. The starting position of the particles is indicated in black, the hit in red. Particles have dimension smaller than the cell size, indicated with a dashed line. Both cells include the hitting position of the particle centre. It is not a matter of the dimension of the particle, it depends upon the fact that the second cell does not share any face plane with the boundary patch plane. Of course the smaller the particle dimension, the smaller the probability that this situation happens.

Your suggestion to include in the hitting test also the neighbour cells should fix this problem but I do not know how much difficult could be to implement such test.

Any comment on this will greatly help us to better understand the algorithm implemented in OpenFOAM.

Thank you in advance,

Franco

dear niklas:
I'm a fresh for the particle-flow interaction, and I need to make some treatment depends on the angle that a particle hit on the wall.
I need to consider the heat exchanges between a particle and the wall the particle hits, and the thermal resistance has some relationship with the hit agle,so and I wonder how can i get the expression about this "angle"?

looking forward for your reply

yours sincerely

Bella

Hi there,
I would like to throw something into this rather old topic. I had very similar problems with wall-collisions and particle sizes like they are depicted here (I am using OF2.3.0 with swak4foam).
Due to a necessarily fine mesh 80% of my particles are bigger than the cells. So what I do is:

1. Check if particle Volume is bigger than cell volume.
2. If it's not, continue with original code.
3. If it is, check how far particle extends into neighbour cells.
4. Repeat process for neighbours neighbours and so on until the new neighbours don't fullfill the "inside particle" criterium anymore (e.g. dist(cell.centre(), p.centre())).
5. Do the originial code for all cells considered inside the particle.

Unfortunatelly this doesn't work so well as it seems that there is a problem with the hitFace() function when the considered cell isn't cellI_ (the particles cell) anymore.


Therefore I simplified it all and simply deactivate the particles if a cell considered inside has any wall-faces.


Maybe someone has some ideas on this topic?

Just to push this topic a little bit and maybe get some new ideas.
Here are two plots for the filtration of particles with the original code and with my additions.

As can be seen the modifications improve the filtration behaviour as the filtration rate doesn't decrease for bigger particles anymore. However, the transition regime between the original code and the modified one remains a problem as there is no interpolation scheme for particles which are in their dimensions close to the cells dimensions.


Does anybody have an idea on how this could be improved?