particle number concentration in OpenFOAM
 

Hallo, everybody

I am using OpenFOAM-1.6.
I combined a two-way coupling Lagrangian class kinematicCloud with sonicFoam, getting a compressible Eulerian-Lagrangian solver.
I want to calculate the particle number concentration in each control volume (N/mesh.V) and the particle flux through the face of the control volume (N/(deltaT*face.area)).

Somebody can give some help on C++ code? how to input the two terms into the main code?:confused: The main code is as follows,
-----------------------------------------------------------------------------------------------------------------------------------------------------
#include "fvCFD.H"
#include "basicPsiThermo.H"
#include "turbulenceModel.H"

#include "basicThermoCloud.H"
#include "basicKinematicCloud.H"

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

int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "createFields.H"
#include "createClouds.H"
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "readTimeControls.H"
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"

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

Info<< "\nStarting time loop\n" << endl;

while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;

#include "readPISOControls.H"
#include "compressibleCourantNo.H"

thermoCloud1.evolve();
thermoCloud1.info();

kinematicCloud1.evolve();
kinematicCloud1.info();

#include "rhoEqn.H"

#include "UEqn.H"

#include "hsEqn.H"


// --- PISO loop

for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
}

turbulence->correct();

rho = thermo.rho();

runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

Info<< "End\n" << endl;

return 0;
}
-----------------------------------------------------------------------------------------------------------------------------------------------------------------

Thank U very much!
Grüß

Hi!
You can create functionObjects and compile them as libraries, the way also other libraries for post processing are created and used in OFoam.
For example you can pick forces.C and forces.H files, that are used to calculate aerodynamic forces, and substitute the code written there with the calculations you need to perform, acting on private member functions and modifying the code accordingly.

After you have created your own libraries, you also have to include them into you controlDict file, this way:

functions
(
forces
{
type forces;
functionObjectLibs ("libforces.so"); //Lib to load
patches
(
intra
extra
); // change to your patch names
CofR (0 0 0);
outputControl timeStep;
outputInterval 50;
rhoInf ...;
pName p;
UName U;
rhoName rhoInf;

}

);

I've taken this example from the way you can calculate the aerodynamic forces.

This procedure will allow you to have the fields you need calculated in your domain, and you will find them also in paraview, so that it will be possible to plot their value.

Marta

Hallo, Marta

Thank you very much. That's a good idea but i have poor C++ knowledge.

In fact when i post-process the results in Paraview, I have found a variable named 'nParticle'.
'nParticle' means particle number in a parcel. Therefore, i want to collect all the 'nParticle' in a control volume, using expression sum(nParticle)/mesh.V(). The error is that nParticle is not declared.

I ever used the similar expression, sum(mesh.V()*alpha1).value(), to collect the liquid volume fraction for interFoam. It worked well.

I know nParticle really exists in my results. The key is I am not familiar with the c++ syntax.

May you give me some more advice?
And other experts can give some help?
Thank you very much!

Dear Wentao, can you post the entire file where you did your code modification here so that I can have a look at it?
Sometimes the only problem is that the functions or parameters you want to use for your operations are not directly accessible from a specific class, so that you need to use the lookup method in a specific way.

See also this thread here:
Accessing fields and dictionaries from within a BC implementation

http://www.cfd-online.com/Forums/ope...mentation.html

Marta

Hallo, Marta
-----------------------------------------------------------------------------------------------------------------------------------------------------
#include "fvCFD.H"
#include "basicPsiThermo.H"
#include "turbulenceModel.H"

#include "basicThermoCloud.H"
#include "basicKinematicCloud.H"

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

int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "createFields.H"
#include "createClouds.H"
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "readTimeControls.H"
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"

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

Info<< "\nStarting time loop\n" << endl;

while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;

#include "readPISOControls.H"
#include "compressibleCourantNo.H"

thermoCloud1.evolve();
thermoCloud1.info();

kinematicCloud1.evolve();
kinematicCloud1.info();

#include "rhoEqn.H"

#include "UEqn.H"

#include "hsEqn.H"


// --- PISO loop

for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
}

turbulence->correct();

rho = thermo.rho();

Info << " pnc: " << sum(nParticle)/mesh.V() << endl;

runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

Info<< "End\n" << endl;

return 0;
}
-----------------------------------------------------------------------------------------------------------------------------------------------------------------

Thank you very much!
Wentao

Ok, the error the compiler gives you derives from the fact that you are trying to access the parameter nParticle, but this parameter is a local variable, so your main file does not recognise it.

This is also the reason why it says that you have not declared it: because it does not know this parameter nor its type.

You should see where inside your classes this parameter and mesh are used and accessible, and for example add in that class a method performing the calculation you need.

Marta

I had a better look at the classes: kinematicCloudI.H and kinematicCloud.H.
You have the template methods in kinematicCloudI.H. I think you can try modifying something there (adding you own piece of code), because you have the possibility to access the particle data and therefore nParticle. You also can have access to mesh data, so it should be fine. This is the piece of code I'm referring to.

template<class ParcelType>
inline const Foam::tmp<Foam::volScalarField>
Foam::KinematicCloud<ParcelType>::rhoEff() const
{
tmp<volScalarField> trhoEff
(
new volScalarField
(
IOobject
(
this->name() + "RhoEff",
this->db().time().timeName(),
this->db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh_,
dimensionedScalar("zero", dimDensity, 0.0)
)
);

scalarField& rhoEff = trhoEff().internalField();
forAllConstIter(typename KinematicCloud<ParcelType>, *this, iter)
{
const ParcelType& p = iter();
const label cellI = p.cell();

rhoEff[cellI] += p.nParticle()*p.mass();
}

rhoEff /= mesh().V();

return trhoEff;
}


Marta

Hallo, Marta

Thank you very much! With your help, everything seems clear!
I will try and inform you my progress!

Wentao

Hallo, Marta

With your help, everything seems clear. But I still want to confirm whether what I have done is right or not!

Take that piece of code you referred to as an example. Function rhoEff() is declared in KinematicCloud.H but defined in KinematicCloudl.H. At the same time KinematicCloud.H is inserted into basicKinematicCloud.H which is applied in my main code as follows.
-----------------------------------------------------------------------------------------------------------------------------------------------------
#include "fvCFD.H"
#include "basicPsiThermo.H"
#include "turbulenceModel.H"
#include "basicKinematicCloud.H"
-----------------------------------------------------------------------------------------------------------------------------------------------------
Therefore, rhoEff() can be regarded as a member function in basicKinematicCloud.H. And then I create an object for basicKinematicCloud class, named kinematicCloud.

In the creatField.H, I also build a volScalarField variable named pncon.
---------------------------------------------------------------------------------------------
Info << "Reading field pncon\n" << endl;
volScalarField pncon
(
IOobject
(
"pncon",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh
);
----------------------------------------------------------------------------------------

In the main code, I assigned kinematicCloud.rhoEff() to pncon as follows.
--------------------------------------------------------------------------------------------------
pncon = kinematicCloud.rhoEff();

runTime.write();
---------------------------------------------------------------------------------------------------
It is compiled successfully.

Hallo, Marta

I am sorry to disturb you again.
When I build a volScalarField, it is compiled very well. But for surfaceScalarField, I have some problem.
I want to calculate the particle number flux (particle number through unit cell face in unit time). It should be a surfaceScalarField.
Therefore, in KinematicCloud_H, i wrote
---------------------------------------------------
inline const tmp<surfaceScalarField> nflux() const;
----------------------------------------------------------
In KinematicCloudl_H, I wrote
----------------------------------------------------------
template<class ParcelType>
inline const Foam::tmp<Foam::surfaceScalarField>
Foam::KinematicCloud<ParcelType>::nflux() const
{
tmp<surfaceScalarField> tnflux
(
new surfaceScalarField
(
IOobject
(
this->name() + "Nflux",
this->db().time().timeName(),
this->db(),
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
mesh_,
dimensionedScalar("zero", dimDensity*dimLength/dimMass/dimTime, 0.0)
)
);

scalarField& nflux = tnflux().internalField();
forAllConstIter(typename KinematicCloud<ParcelType>, *this, iter)
{
const ParcelType& p = iter();
const label faceI = p.face();

nflux[faceI] += p.nParticle()*(p.U()&mesh().Sf())/(mag(p.U())*this->db().time().deltaT());
}

nflux /= (mesh().magSf()*mesh().magSf());

return tnflux;
}
----------------------------------------------------------------------
And compile it, the error is
----------------------------------------------------------------------
/OpenFOAM/OpenFOAM-1.6.x/src/lagrangian/intermediate/lnInclude/KinematicCloudI.H:393: error: no match for ‘operator+=’ in ‘((Foam::scalarField*)nflux)->Foam::Field<double>::<anonymous>.Foam::List<doubl e>::<anonymous>.Foam::UList<T>::operator[] [with T = double](faceI) += Foam::operator/(const Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh> >&, const Foam::dimensioned<double>&) [with Type = double, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh](((const Foam::dimensioned<double>&)((const Foam::dimensioned<double>*)(& Foam::operator*(const Foam::dimensioned<double>&, const Foam::dimensioned<Type>&) [with Type = double](((const Foam::dimensioned<double>&)((const Foam::dimensioned<double>*)(& Foam::TimeState::deltaT() const()))))))))’
make: *** [Make/linux64GccDPOpt/intermediateFoam.o] Error 1
--------------------------------------------------------------------------------------
I can't find an example on how to introduce a surfaceScalarField, therfore just follow the volScalarField syntax......
May you help me find the problem? And other experts?
Thank you very much!!!

Dear Wentao, first of all i will go back to the rhoEff and kinematicCloud topic. You said that your code compiles, but apart from that, are you sure that it works correctly?
Because i think that, if you create an object and then call function rhoEff from your main file, what you get out of it is what is returned by the function, which is not what you are trying to calculate.

The method returns in fact: trhoEff.

Have you verified this aspect?

I will have a look at the rest of the code for you!

Regards

Marta

Dear Wentao, I think that the error you get from the surfaceScalarField operation, is due to the fact that you are actually summing a parameter defined as a surfaceScalarField (nFlux) to another one of a different type. That's why it says: 'no match for operator +'.

To perform your calculation correctly, i think you should cast the parameter which does not correspond to a surfaceScalarField as a type.


Hope this helps!


Marta

Hallo, dear Marta

Thank you very much for your reply!
You are right! trhoEff is not what i am trying to calculate. I just made an excecise using it.
In fact, i built another function for particle concentration (particle.number/cell.volume)
in KinematicCloud_H, i wrote
---------------------------------------------------
inline const tmp<volScalarField> pnc() const;
----------------------------------------------------------
In KinematicCloudl_H, I wrote
----------------------------------------------------------
forAllConstIter(typename KinematicCloud<ParcelType>, *this, iter)
{
const ParcelType& p = iter();
const label cellI = p.cell();

pnc[cellI] += p.nParticle());
}

pcn /= mesh().V();

return tpnc;
}

As for number flux, you advised me to ''cast the parameter which does not correspond to a surfaceScalarField as a type''.
You mean I should change the type of other parameters such as p.number() and deltaT()? But the types of these parameters have fixed in other classes. I don't understand what you said. May you give me an example?
I have checked other threads on how to treat a surfaceScalarField, for example, phi()= U(interpolate)&mesh().Sf(). here U should be the velocity of continuous phase.
But for particle number flux, may I use an interpolated particle velocity?
Best wishes!
Wentao

scalar-flux and nParticle Gradient normal to face?
 

Wentao,

What happened, did you figure out your answer? I am interested in getting a surfaceScalarField of scalar-flux of particles through a face along with the gradient of the particle concentration normal to the face of each cell. I'm looking at KinematicCloudI.H for some inspiration and I've only spent a few hours looking at the problem. Thoughts...suggestions?

Dan

how to do two way coupling
 

Dear all,
I have done some small lagrangian particle tracking using solidParticleFoam and icoLagrangianFoam. These simulations are one way coupled. Now I would like to do two way coupling. This I would like to do for bubble column case where gas bubbles can be injected as lagrangian particle. I am wondering if some one can tell me how to do two way coupling.
Thanks in advance.
Raghu

What do you mean with two-way coupling? The particles influence the fluid? As far as I remember the different incarnations of icoLagrangianFoam (and most lagrangian solvers in the distribution) support that. Have a look at UEqn.H of whatever solver you're using. If there is a source term in the equation then it is already implemented. Maybe it is not switched on

Bernhard

Hello Fabian and Marta,

I see you guys were making an attempt to count number of particle coming out of a domain. I am doing a particle simulation using icoUncoupledKinematicParcelFoam and I want to calculate the number of particles coming out of each outlet. I am stuck on that problem although simulation is running fine. Can you please guide me how to calculate particles moving across boundary patches?

Regards,
Mayank.

Hello,

You can use local interaction for patchInteraction of particles. Then you need to specify particle interaction type for each BC. Specify outlet(Use BC name in your case) as
outlet
{
type escape;
}
This will report number of particles escaping from outlet BC in log file.

Regards,
-Yogesh

Hi Yogesh,

Thanks a lot for your reply. I was looking for a hint for many days.

Is this the complete piece of code or anything else has to be added?

Could you be a little more specific where to place this piece of code, I mean inside solver or somewhere in kinematicCloudProperties or associated file?

I apologize for asking for specific information as I am not very good at C++ or OpenFOAM but I am trying my level best to gain proficiency.

Regards,
Mayank