CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

particle number concentration in OpenFOAM

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   September 28, 2010, 07:08
Default particle number concentration in OpenFOAM
  #1
New Member
 
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15
wentao is on a distinguished road
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? 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üß
wentao is offline   Reply With Quote

Old   September 28, 2010, 12:17
Default
  #2
Member
 
Marta's Avatar
 
Marta Lazzarin
Join Date: Jun 2009
Location: Italy
Posts: 71
Rep Power: 16
Marta is on a distinguished road
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
Marta is offline   Reply With Quote

Old   September 29, 2010, 04:43
Default
  #3
New Member
 
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15
wentao is on a distinguished road
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!
wentao is offline   Reply With Quote

Old   September 29, 2010, 05:59
Default
  #4
Member
 
Marta's Avatar
 
Marta Lazzarin
Join Date: Jun 2009
Location: Italy
Posts: 71
Rep Power: 16
Marta is on a distinguished road
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
Marta is offline   Reply With Quote

Old   September 29, 2010, 06:26
Default
  #5
New Member
 
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15
wentao is on a distinguished road
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
wentao is offline   Reply With Quote

Old   September 29, 2010, 06:55
Default
  #6
Member
 
Marta's Avatar
 
Marta Lazzarin
Join Date: Jun 2009
Location: Italy
Posts: 71
Rep Power: 16
Marta is on a distinguished road
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
Marta is offline   Reply With Quote

Old   September 29, 2010, 07:20
Default
  #7
Member
 
Marta's Avatar
 
Marta Lazzarin
Join Date: Jun 2009
Location: Italy
Posts: 71
Rep Power: 16
Marta is on a distinguished road
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
Alaska1964 and chpjz0391 like this.
Marta is offline   Reply With Quote

Old   September 29, 2010, 09:08
Default
  #8
New Member
 
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15
wentao is on a distinguished road
Hallo, Marta

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

Wentao
wentao is offline   Reply With Quote

Old   September 29, 2010, 15:16
Default
  #9
New Member
 
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15
wentao is on a distinguished road
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.
wentao is offline   Reply With Quote

Old   September 30, 2010, 14:36
Default
  #10
New Member
 
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15
wentao is on a distinguished road
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>:perator[] [with T = double](faceI) += Foam:perator/(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:perator*(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!!!
wentao is offline   Reply With Quote

Old   September 30, 2010, 16:49
Default
  #11
Member
 
Marta's Avatar
 
Marta Lazzarin
Join Date: Jun 2009
Location: Italy
Posts: 71
Rep Power: 16
Marta is on a distinguished road
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
Marta is offline   Reply With Quote

Old   September 30, 2010, 17:02
Default
  #12
Member
 
Marta's Avatar
 
Marta Lazzarin
Join Date: Jun 2009
Location: Italy
Posts: 71
Rep Power: 16
Marta is on a distinguished road
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
Marta is offline   Reply With Quote

Old   September 30, 2010, 18:41
Default
  #13
New Member
 
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15
wentao is on a distinguished road
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;
}
wentao is offline   Reply With Quote

Old   September 30, 2010, 19:04
Default
  #14
New Member
 
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15
wentao is on a distinguished road
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
wentao is offline   Reply With Quote

Old   January 17, 2011, 19:12
Default scalar-flux and nParticle Gradient normal to face?
  #15
Senior Member
 
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 621
Rep Power: 0
chegdan will become famous soon enoughchegdan will become famous soon enough
Quote:
Originally Posted by wentao View Post
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
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
chegdan is offline   Reply With Quote

Old   February 3, 2011, 03:04
Default how to do two way coupling
  #16
Member
 
Raghavendra
Join Date: Mar 2009
Location: Goteborg, Sweden
Posts: 95
Rep Power: 17
raagh77 is on a distinguished road
Send a message via Yahoo to raagh77 Send a message via Skype™ to raagh77
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
raagh77 is offline   Reply With Quote

Old   February 3, 2011, 11:15
Default
  #17
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,225
Rep Power: 51
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by raagh77 View Post
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
gschaider is offline   Reply With Quote

Old   August 22, 2013, 00:06
Default
  #18
Member
 
India
Join Date: Oct 2012
Posts: 84
Rep Power: 13
mayank.dce2k7 is on a distinguished road
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.
mayank.dce2k7 is offline   Reply With Quote

Old   August 22, 2013, 01:28
Post
  #19
Senior Member
 
Yogesh Bapat
Join Date: Oct 2010
Posts: 102
Rep Power: 15
ybapat is on a distinguished road
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
ybapat is offline   Reply With Quote

Old   August 22, 2013, 12:41
Default
  #20
Member
 
India
Join Date: Oct 2012
Posts: 84
Rep Power: 13
mayank.dce2k7 is on a distinguished road
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

Last edited by mayank.dce2k7; August 22, 2013 at 15:22.
mayank.dce2k7 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
mesh file for flow over a circular cylinder Ardalan Main CFD Forum 7 December 15, 2020 14:06
Problem with decomposePar tool vinz OpenFOAM Pre-Processing 18 January 26, 2011 03:17
Particle Tracking Batis CFX 2 October 6, 2010 16:20
[Commercial meshers] Trimmed cell and embedded refinement mesh conversion issues michele OpenFOAM Meshing & Mesh Conversion 2 July 15, 2005 05:15
particle courant number Michiel Siemens 1 August 10, 2003 06:50


All times are GMT -4. The time now is 19:39.