CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM (http://www.cfd-online.com/Forums/openfoam/)
-   -   Find a particle's distance to a wall (http://www.cfd-online.com/Forums/openfoam/76042-find-particles-distance-wall.html)

gruber May 12, 2010 10:23

Find a particle's distance to a wall
 
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 May 14, 2010 03:31

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 February 14, 2012 09:34

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 February 14, 2012 09:43

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

michaelb February 14, 2012 09:50

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 February 14, 2012 10:08

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 February 14, 2012 13:33

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 February 14, 2012 16:34

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.

Code:


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

Code:

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


michaelb February 15, 2012 06:59

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 ?

Code:

#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 February 15, 2012 07:22

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


All times are GMT -4. The time now is 00:51.