|
[Sponsors] |
September 28, 2010, 06:08 |
particle number concentration in OpenFOAM
|
#1 |
New Member
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15 |
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üß |
|
September 28, 2010, 11:17 |
|
#2 |
Member
Marta Lazzarin
Join Date: Jun 2009
Location: Italy
Posts: 71
Rep Power: 16 |
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 |
|
September 29, 2010, 03:43 |
|
#3 |
New Member
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15 |
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! |
|
September 29, 2010, 04:59 |
|
#4 |
Member
Marta Lazzarin
Join Date: Jun 2009
Location: Italy
Posts: 71
Rep Power: 16 |
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 |
|
September 29, 2010, 05:26 |
|
#5 |
New Member
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15 |
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 |
|
September 29, 2010, 05:55 |
|
#6 |
Member
Marta Lazzarin
Join Date: Jun 2009
Location: Italy
Posts: 71
Rep Power: 16 |
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 |
|
September 29, 2010, 06:20 |
|
#7 |
Member
Marta Lazzarin
Join Date: Jun 2009
Location: Italy
Posts: 71
Rep Power: 16 |
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 |
|
September 29, 2010, 08:08 |
|
#8 |
New Member
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15 |
Hallo, Marta
Thank you very much! With your help, everything seems clear! I will try and inform you my progress! Wentao |
|
September 29, 2010, 14:16 |
|
#9 |
New Member
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15 |
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. |
|
September 30, 2010, 13:36 |
|
#10 |
New Member
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15 |
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!!! |
|
September 30, 2010, 15:49 |
|
#11 |
Member
Marta Lazzarin
Join Date: Jun 2009
Location: Italy
Posts: 71
Rep Power: 16 |
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 |
|
September 30, 2010, 16:02 |
|
#12 |
Member
Marta Lazzarin
Join Date: Jun 2009
Location: Italy
Posts: 71
Rep Power: 16 |
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 |
|
September 30, 2010, 17:41 |
|
#13 |
New Member
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15 |
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; } |
|
September 30, 2010, 18:04 |
|
#14 |
New Member
Fabian
Join Date: Sep 2010
Posts: 8
Rep Power: 15 |
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 |
|
January 17, 2011, 18:12 |
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 |
Quote:
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 |
||
February 3, 2011, 02:04 |
how to do two way coupling
|
#16 |
Member
|
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 |
|
February 3, 2011, 10:15 |
|
#17 | |
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,225
Rep Power: 51 |
Quote:
Bernhard |
||
August 21, 2013, 23:06 |
|
#18 |
Member
India
Join Date: Oct 2012
Posts: 84
Rep Power: 13 |
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. |
|
August 22, 2013, 00:28 |
|
#19 |
Senior Member
Yogesh Bapat
Join Date: Oct 2010
Posts: 102
Rep Power: 15 |
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 |
|
August 22, 2013, 11:41 |
|
#20 |
Member
India
Join Date: Oct 2012
Posts: 84
Rep Power: 13 |
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 14:22. |
|
|
|
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 13:06 |
Problem with decomposePar tool | vinz | OpenFOAM Pre-Processing | 18 | January 26, 2011 02:17 |
Particle Tracking | Batis | CFX | 2 | October 6, 2010 15:20 |
[Commercial meshers] Trimmed cell and embedded refinement mesh conversion issues | michele | OpenFOAM Meshing & Mesh Conversion | 2 | July 15, 2005 04:15 |
particle courant number | Michiel | Siemens | 1 | August 10, 2003 05:50 |