CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Post-Processing (https://www.cfd-online.com/Forums/openfoam-post-processing/)
-   -   Lagrangian data - particle velocity post processing (https://www.cfd-online.com/Forums/openfoam-post-processing/222030-lagrangian-data-particle-velocity-post-processing.html)

Ake November 8, 2019 10:14
Lagrangian data - particle velocity post processing
 
Hi everyone,


i am doing LES and i just added lagrangian particles to pimpleFoam (openfoam 7). But now i am struggeling with the post processing, which isn't as easy as I thought/expected.:confused: :(

I am able to see the lagrangian data in paraView. But i have no idea how to sample the particle velocity to get (for example) an axial particle velocity profile. Additionally i need the average of the particle velocity, what could be done easily in the controlDict if i have a field. So what i need is a field of the particle velocity with that i am able to create the time average.


In the Cloudfunction "VoidFraction" (VoidFraction.C) there is lagrangian data "written in a field" i guess.
--> theta[p.cell()] += dt*p.nParticle()*p.volume();

So i have tried to do it the same way and edit the VoidFraction cloudfunction with this or something similar
--> Uparticle[p.cell()] += dt*mag(p.U());
but this didnt work properly. The solver compiles but the magnitude of the velocity of the created field is much too high (about 10⁹). I think there is something wrong in the transfer of the "point" velocity to the "cell" velocity.




Does anyone have any idea how this could work?

I'm not sure this is the right way. Maybe someone has a better idea to sample the particle velocity at any position? Is there an existing CloudFunction which can be used for that?


Thanks in advance

fresty November 11, 2019 05:01
Have you tried using 'calculate' to acquire p.U(component) and 'integrate variables' in ParaView?

Ake November 13, 2019 06:35
Quote:
Originally Posted by fresty (Post 749436)
Have you tried using 'calculate' to acquire p.U(component) and 'integrate variables' in ParaView?
Thanks for your answer. I havent tried it that way. But finally, I solved my problem myself. Here is my solution for all those who are interested.



I created a new vector field 'UPar' in the file KinematicCloudI.H similar to the given scalar field 'theta'. Additionally i created a scalar field of the particles sum 'SumPar' (=parcels sum). I need this field because it is possible that one cell contains more than one particle and then the averaged velocity is calculated.


Code:
template<class CloudType>
inline const Foam::tmp<Foam::volVectorField>
Foam::KinematicCloud<CloudType>::UPar() const
{
    tmp<volVectorField> tUPar
    (
        volVectorField::New
        (
            this->name() + ":UPar",
            mesh_,
            dimensionedVector(dimVelocity, vector(0,0,0)),
            extrapolatedCalculatedFvPatchVectorField::typeName
        )
    );


    volVectorField& UPar = tUPar.ref();
    forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
    {
        const parcelType& p = iter();
        const label celli = p.cell();
        UPar[celli] += (p.U())/(p.nParticle());
    }
    UPar.primitiveFieldRef() /= 1;//mesh_.V();
    UPar.correctBoundaryConditions();

    return tUPar;
}



 
template<class CloudType>
inline const Foam::tmp<Foam::volScalarField>
Foam::KinematicCloud<CloudType>::SumPar() const
{
    tmp<volScalarField> tSumPar
    (
        volScalarField::New
        (
            this->name() + ":SumPar",
            mesh_,
            dimensionedScalar(dimless, 0),
            extrapolatedCalculatedFvPatchScalarField::typeName
        )
    );


    volScalarField& SumPar = tSumPar.ref();
    forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
    {
        const parcelType& p = iter();
        const label celli = p.cell();
        SumPar[celli] += p.nParticle();
   
    }

    SumPar.primitiveFieldRef() /= 1;//mesh_.V();
    SumPar.correctBoundaryConditions();

    return tSumPar;
}
Afterwards 'UPar' is divided by 'SumPar' in in my pimpleLPT.C file if there is at least one particle in it.


Code:
       

    volScalarField SumParTemp=(kinematicCloud.SumPar())();
    volVectorField UParTemp=(kinematicCloud.UPar())();

    forAll(mesh.C(), i)
    {
    if(SumParTemp[i]>0)
        {
            UPar[i] = UParTemp[i]/SumParTemp[i];
        }
    else
        {
            UPar[i] = UParTemp[i];
        }   
    }
    UPar.correctBoundaryConditions(); 
    SumPar = kinematicCloud.SumPar();
    SumPar.correctBoundaryConditions();
But be carefull.... in paraview the point data is not correct. This is because of the interpolation of the velocities. You have to look at the 'cell data'.



This solution seems to be very good for my investigations. Nevertheless, if you have suggestions for improvement, I would be happy if you share them. :)

as020002 March 3, 2020 02:55
Quote:
Originally Posted by Ake (Post 749662)
Thanks for your answer. I havent tried it that way. But finally, I solved my problem myself. Here is my solution for all those who are interested.



I created a new vector field 'UPar' in the file KinematicCloudI.H similar to the given scalar field 'theta'. Additionally i created a scalar field of the particles sum 'SumPar' (=parcels sum). I need this field because it is possible that one cell contains more than one particle and then the averaged velocity is calculated.


Code:
template<class CloudType>
inline const Foam::tmp<Foam::volVectorField>
Foam::KinematicCloud<CloudType>::UPar() const
{
    tmp<volVectorField> tUPar
    (
        volVectorField::New
        (
            this->name() + ":UPar",
            mesh_,
            dimensionedVector(dimVelocity, vector(0,0,0)),
            extrapolatedCalculatedFvPatchVectorField::typeName
        )
    );


    volVectorField& UPar = tUPar.ref();
    forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
    {
        const parcelType& p = iter();
        const label celli = p.cell();
        UPar[celli] += (p.U())/(p.nParticle());
    }
    UPar.primitiveFieldRef() /= 1;//mesh_.V();
    UPar.correctBoundaryConditions();

    return tUPar;
}



 
template<class CloudType>
inline const Foam::tmp<Foam::volScalarField>
Foam::KinematicCloud<CloudType>::SumPar() const
{
    tmp<volScalarField> tSumPar
    (
        volScalarField::New
        (
            this->name() + ":SumPar",
            mesh_,
            dimensionedScalar(dimless, 0),
            extrapolatedCalculatedFvPatchScalarField::typeName
        )
    );


    volScalarField& SumPar = tSumPar.ref();
    forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
    {
        const parcelType& p = iter();
        const label celli = p.cell();
        SumPar[celli] += p.nParticle();
   
    }

    SumPar.primitiveFieldRef() /= 1;//mesh_.V();
    SumPar.correctBoundaryConditions();

    return tSumPar;
}
Afterwards 'UPar' is divided by 'SumPar' in in my pimpleLPT.C file if there is at least one particle in it.


Code:
       

    volScalarField SumParTemp=(kinematicCloud.SumPar())();
    volVectorField UParTemp=(kinematicCloud.UPar())();

    forAll(mesh.C(), i)
    {
    if(SumParTemp[i]>0)
        {
            UPar[i] = UParTemp[i]/SumParTemp[i];
        }
    else
        {
            UPar[i] = UParTemp[i];
        }   
    }
    UPar.correctBoundaryConditions(); 
    SumPar = kinematicCloud.SumPar();
    SumPar.correctBoundaryConditions();
But be carefull.... in paraview the point data is not correct. This is because of the interpolation of the velocities. You have to look at the 'cell data'.



This solution seems to be very good for my investigations. Nevertheless, if you have suggestions for improvement, I would be happy if you share them. :)
Hi Kevin,

This is very helpful, may I ask where did you put the second code in pimpleLPT.C? I added it below the "parcels.evolve" within the runTime.run loop and it didn't work.

Thanks,
RIck

Ake March 3, 2020 03:17
Hi Rick,


I put it in the same place. Within the runtime loop and outside the pimple loop.


Maybe you can try it without the temporary fields and the particle sum in a first step. Just add a velocity vectorfield "UPar2" in the createFieldsDict and then add this below the .evolve in your pimpleLPT.C file:

Code:
UPar2 = kinematicCloud.UPar();

as020002 March 3, 2020 08:41
Quote:
Originally Posted by Ake (Post 760224)
Hi Rick,


I put it in the same place. Within the runtime loop and outside the pimple loop.


Maybe you can try it without the temporary fields and the particle sum in a first step. Just add a velocity vectorfield "UPar2" in the createFieldsDict and then add this below the .evolve in your pimpleLPT.C file:

Code:
UPar2 = kinematicCloud.UPar();
Dear Kevin,

Yes I think the code in .C file is fine but I still haven't figured out where the problems are, could you post your createFields.H? It would be very appreciated.

Cheers,
Rick

as020002 March 3, 2020 09:27
1 Attachment(s)
I think I should post my error message here:

Ake March 3, 2020 10:13
It could be the name of your Cloud (.evolve)


The name of my kinematicCloud is kinematicCloud. This is my evolve in the .C file:
Code:
kinematicCloud.evolve();
And this is my createFieldsDict:
Code:
const word kinematicCloudName
(
    args.optionLookupOrDefault<word>("cloudName", "kinematicCloud")
);
Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
basicKinematicCollidingCloud kinematicCloud
(
    kinematicCloudName,
    rho,
    U,
    mu,
    g
);

as020002 March 3, 2020 11:19
Quote:
Originally Posted by Ake (Post 760286)
It could be the name of your Cloud (.evolve)


The name of my kinematicCloud is kinematicCloud. This is my evolve in the .C file:
Code:
kinematicCloud.evolve();
And this is my createFieldsDict:
Code:
const word kinematicCloudName
(
    args.optionLookupOrDefault<word>("cloudName", "kinematicCloud")
);
Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
basicKinematicCollidingCloud kinematicCloud
(
    kinematicCloudName,
    rho,
    U,
    mu,
    g
);
Dear Kevin,

Thank you it successfully compiled. However I tested it with some cases it ended up like another error, it's like there are problems in dimension:, and my createFields.H is of attachment.:(

Cheers,
Rick

as020002 March 3, 2020 11:20
1 Attachment(s)
Sorry, I forgot to attach the image, there is it.

as020002 March 4, 2020 03:14
After I correct some minor syntax mistakes and make some declaration in source code files, now it works perfectly now.

Thank you Kevin:)

arch.mr August 16, 2020 07:24
1 Attachment(s)
Dear Cheung Wing Ki
I also tried to create a scalar field of the particles sum (particle number per unit volume)
I followed Kevin's codes and tried to compile

kinematicCloudI.H

Code:
template<class CloudType>
inline const Foam::tmp<Foam::volScalarField>
Foam::KinematicCloud<CloudType>::SumPar() const
{
    tmp<volScalarField> tSumPar
    (
        volScalarField::New
        (
            this->name() + ":SumPar",
            mesh_,
            dimensionedScalar(dimless, 0),
            extrapolatedCalculatedFvPatchScalarField::typeName
        )
    );


    volScalarField& SumPar = tSumPar.ref();
    forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
    {
        const parcelType& p = iter();
        const label celli = p.cell();
        SumPar[celli] += p.nParticle();
   
    }

    SumPar.primitiveFieldRef() /= mesh_.V(); // 1;
    SumPar.correctBoundaryConditions();

    return tSumPar;
}
and add this right after rhoEff
KinematicCloud.H

Code:
inline const tmp<volScalarField> rhoEff() const;

inline const tmp<volScalarField> SumPar() const;
but I got following error message...
intL4Foam.C: In function ‘int main(int, char**)’:
intL4Foam.C:188:9: error: ‘SumPar’ was not declared in this scope
SumPar = kinematicCloud.SumPar();
^~~~~~
intL4Foam.C:188:33: error: ‘Foam::basicKinematicCollidingCloud {aka class Foam::CollidingCloud<Foam::KinematicCloud<Foam::Cl oud<Foam::CollidingParcel<Foam::KinematicParcel<Fo am::particle> > > > >}’ has no member named ‘SumPar’
SumPar = kinematicCloud.SumPar();
^~~~~~


Does anybody can help me?
Thank you

arch.mr August 16, 2020 07:54
1 Attachment(s)
I tried to get particle density output per volume
Add code in kinematicCloud.H, kinematicCloudI.H, main.C file and createFields.H
And got following error messages..

Anybody please can help me?


Thank you


All times are GMT -4. The time now is 03:21.