[gruber]

Hello,

I am implementing an algorithm for particle movement. Therefore, I need to calculate the distance between a particle and the wall in the particle's velocity direction (if there is a wall).

I looked for the solution in the Particle and Cloud files but I couldn't find it.

Thank you for your help

[gruber]

Maybe I should rephrase my question. How can three different points of one boundary patch be read out to describe a plane equation? And how can this be done for different boundary patches?

Thank you for your help

[michaelb]

Hi Michael,

I'm interested in doing the same thing ie finding particle distance to the wall. Did you come up with a solution to your problem ? I would be thankful if you could share any hints !

Best regards,

Michael

[niklas]

Its quite easy, but
can I first ask why you need it?

[michaelb]

hi Niklas,

thanks for adressing the problem.
I'm glad you're saying that it is easy, although I'm not sure it will be for me.
I need to get the particle distance to the wall to get some statistics about particle repartition in a complex domain.
Have you done that before? What would be the starting point ? I was thinking about the wallDist function but it spits out a volScalarField which isn't very suitable here...

Cheers,

Michael

[niklas]

ok, but just so you know it is an extremely time-consuming operation that you want to perform

first. the easy part

the general formula for the distance between a point and a plane/face

x = particle position
C = face centre (or arbitrary point on a face)
n = wall/face normal, where n must be normalized, i.e | n | = 1

distance = | n .dot. ( x - C) |

next is the time-consuming part.
in order to find the closest distance to a wall you need to extract all the wall-faces and calculate
the distance to each face and pick the closest one.

if you are doing this in parallel it will be slightly trickier, since you need to get the wall-faces on the
other processors and check those as well.

foam already have access functions that give you access to patches and their face normals and face centres,
which you should find quite easy with a little searching.

[michaelb]

thanks Niklas,

yep first part was clear,

the second part will require me to dive into the code, which I've pretty much avoided until now... I'll give it a try !

If you are used to deal with lagrangian particle simulation, would you please have a look here and advise ? :
http://www.cfd-online.com/Forums/ope...tml#post344331

Thanks

Michael

[niklas]

I realized that the way I described above won't actually work.
Its probably better to just check the distance to all the points on the face.
Still, the code to do it might be of some help.

const polyBoundaryMesh& patches = mesh.boundaryMesh();

forAllConstIter
(
    typename Cloud<SprayParcel<ReactingParcel<ThermoParcel<KinematicParcel<particle> > > > >, 
    parcels, 
    iter
)
{
    // get the position for the particle                                                                                                                                            
    const vector& position = iter().position();

    scalar distance = GREAT;
    // loop over all patches                                                                                                                                                        
    forAll(patches, patchI)
    {
        if (isA<wallPolyPatch>(patches[patchI]))
        {
            const polyPatch& patch = patches[patchI];
            const vectorField& n = patch.faceNormals();
            const vectorField& C = patch.faceCentres();
            forAll(C, i)
            {
                vector normal = n[i]/mag(n[i]);
                //scalar di = mag(position - C[i]); // distance to face centre
                scalar di = mag( normal & ( position - C[i] ) );
                distance = min(distance, di);
            }
        }
    }
    Info << "distance to wall is " << distance << endl;
}

The reason it won't work is because of this...
in the 'figure' below, it will believe its closer to wall A than wall B

                                      wall A
                             --------------
         x                  |
                             |
                             |  wall B

[michaelb]

Hello Niklas,

I don't have such sharp edges in my domain, so I guess it will be ok for me...
I cannot get the following compiling though : I 've a problem with the wallDist[k] initialization.
wallDist[] here is the IOField I would like to write in timestep/lagrangian/kinematicCloud/wallDist.

Do you have a suggestion on how to do that correctly ?

#include "argList.H"
#include "Cloud.H"
#include "IOdictionary.H"
#include "fvMesh.H"
#include "Time.H"
#include "timeSelector.H"
#include "OFstream.H"
#include "passiveParticleCloud.H"
#include "IOPosition.H"

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    timeSelector::addOptions();
    #include "addRegionOption.H"

    #include "setRootCase.H"

    #include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
    #include "createNamedMesh.H"
    #include "createFields.H"

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "Scanning times to determine distance to wall for cloud " << cloudName
        << nl << endl;



    const polyBoundaryMesh& patches = mesh.boundaryMesh();
    
    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);
        Info<< "Time = " << runTime.timeName() << endl;

        Info<< "    Reading particle positions" << endl;
        passiveParticleCloud myCloud(mesh, cloudName);
        label np = myCloud.size();
        Info<< "    Read " << returnReduce(np, sumOp<label>())
            << " particles" << endl;
            
        IOField<label> wallDist(myCloud.fieldIOobject("wallDist", IOobject::NO_READ), np);
        label k = 0;
        forAllConstIter(passiveParticleCloud, myCloud, iter)
        {
            // get the position for the particle                                                                                                                                            
            const vector& position = iter().position();

            wallDist[k] = GREAT;  //  <----- BAD LINE
            //loop over all patches                                                                                                                      
            forAll(patches, patchI)
            {
                if (isA<wallPolyPatch>(patches[patchI]))  
                {
                    const polyPatch& patch = patches[patchI];
                    const vectorField& n = patch.faceNormals();
                    const vectorField& C = patch.faceCentres();
                    forAll(C, i)
                        {
                            vector normal = n[i]/mag(n[i]);
                            scalar di = mag( normal & ( position - C[i] ) );
                            wallDist[k] = min(wallDist[k], di);
                        }
                }
            }
            //Info << "Distance to wall is " << wallDist[k] << endl;
            k++;
        }
        wallDist.write();
    }
    return 0;
}

Cheers,

Michael

[michaelb]

Got it !
The error wasn't where I expected it to be. I had to replace :
IOField<label> wallDist(myCloud.fieldIOobject("wallDist", IOobject::NO_READ), np);
by
IOField<scalar> wallDist(myCloud.fieldIOobject("wallDist", IOobject::NO_READ), np);

cheers,

Michael