[shchepan]

Niklas, you were right, this check does not work. But I cannot understand why.

                vector Sf = Sfbc[j]; //area of the face    
                vector pHit = position_ + lam*(endPosition - position_);
                //pHit is the radius vector of hit point
                scalar S = 0.0;            
 
                forAll(facePoints, fi)
                {
                    label fi1 = (fi + 1) % facePoints.size();
                    vector p1 = points[facePoints[fi1]];
                    vector p0 = points[facePoints[fi]];
                    //lHit1 and lHit2 are vectors between face points and hit point
                    vector lHit1 = p1 - pHit; 
                    vector lHit2 = p0 - pHit;
                    
                    //I calculate the area of triangle between lHit1 and lHit2
                    scalar S += mag(lHit1 ^ lHit2) / 2.0;
                }

                //pHit can be within edges only if sum of the areas of triangles are equal to area of the face 
                if (S == mag(Sf))
                {
                    faces.append(facei);
                }

S is more than mag(Sf) for all faces!

[niklas]

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)

    const faceList& faces = mesh.faces();
    const labelList& f = faces[facei];
    label nPoints = f.size();

    for(label i=0; i<nPoints; i++)
    {
        // construct the face normal of the triangle containing points i and j=i+1
        label j = (i + 1) % nPoints;
        const vector& pi = points[f[i]];
        const vector& pj = points[f[j]];

        vector line1 = pj - pi;
        vector line2 = Cf - pj;
        vector line3 = pi - Cf;
        vector n = line2 ^ line3;
        n /= mag(n);

        // calculate the distance-fraction to the triangle-plane
        scalar lambdaNominator = (Cf - from) & n;
        scalar lambdaDenominator = (to - from) & n;
        // check if trajectory is parallel to face
        if (mag(lambdaDenominator) < 1.0e-30)
        {
            if (lambdaDenominator < 0.0)
            {
                lambdaDenominator = -1.0e-30;
            }
            else
            {
                lambdaDenominator = 1.0e-30;
            }
        }

        scalar lam =  lambdaNominator/lambdaDenominator;

        // check if the hit is within the triangle edges
        if ( (lam >= 0.0) && (lam <= 1.0))
	{
            vector px = from + lam*(to -from);

            // calculate the outward pointing triangle edge normals
            vector normal1 = line3 - ( (line3 & line1)/(line1 & line1) )*line1;
            vector normal2 = line1 - ( (line1 & line2)/(line2 & line2) )*line2;
            vector normal3 = line2 - ( (line2 & line3)/(line3 & line3) )*line3;
            normal1 /= mag(normal1);
            normal2 /= mag(normal2);
            normal3 /= mag(normal3);

            // calculate the distance to the edges
            scalar d1 = (px - pi) & normal1;
            scalar d2 = (px - pj) & normal2;
            scalar d3 = (px - Cf) & normal3;
 
            // check if point is inside triangle and allow for small overlap on edges
	    // inside means that all d's, 1-3, are negative.
            // d2 and d3 are internal edges, so tolerance can be increase further for those if needed
            if ( (d1 <= myTol) && (d2 <= myTol) && (d3 <= myTol))
	    {
	        hitFace = true;
            }
        }
    }

[shchepan]

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!)

[niklas]

eehhh....
example:

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

px - pi = (1 0 0)

how can this vector be from px to pi?

[shchepan]

Originally Posted by niklas

eehhh....

I am really sorry, I am little bit tired(

[shchepan]

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

[ambros]

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

[niklas]

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.

[francescomarra]

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

[white_friday]

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

[Kabra]

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?

[Kabra]

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?