[Atzori]

Dear all,
I’d like to implement a particles stochastic model which requires to add random generators in the integration both of U_particle and of X_particle.
I though that the best option could be to find where U_p and X_p are updated to insert there my equations, but I didn't managed to find where the position of the single parcel is updated...
I shortly summed up what I considered the most important passages in the "flow" of the solver:

1. kinematicCloud.evolve() (called in the solver)
2. TrackingData created (called in evolve() in KinematicCloud.C)
3. solve(td) (called in evolve() in KinematicCloud.C)
4. preEvolve (called in solve() in KinematicCloud.C)
5. evolveCloud(td) (called in solve() in KinematicCloud.C)
6. td.cloud().motion(td); (called in evolveCloud(td) in KinematicCloud.C)
7. CloudType::move(td, solution_.trackTime()); (called in motion(td) in KinematicCloud.C)
8. bool keepParticle = p.move(td, trackTime); (called in move() in Cloud.C)
9. p.calc(td, dt, celli); (called in move() in KinematicParcel.C)
10. this->U_ = calcVelocity(td, dt, celli, Re, muc_, mass0, Su, dUTrans, Spu); (called in move() in KinematicParcel.C) (here is called the integrator for U, where I plan to set the new equation for U_p)
11. ...
12. ...
    

In the following I hoped to find something equivalent to:

new_X = current_X+ Delta_t * new_U

Where I can find it?

It seems to me that "something happens" in the TrackingData, but I did not managed to divide the "physical part" (i.e.: the new position) and "service part" (interaction with boundaries, change of processor and so on...)

On addition: do you have any suggestions for a better approach in coding?
(e.g.: I think that, if only U_p had to be changed, the best choice would have been to create a new force model, but it seems to me that it cannot help with X_p)
Please, consider that the "strategic task" is to obtain somethig as (in a very very very idealized picture):

X_p(n+1)=X_p(n) + \Delta t * U(n) + Random1
U_p(n+1)=(Drag term) + Random2


Thanks in advance!!

[A_Pete]

Hi Atzori,

go to "particleTemplates.C". Here, you will find the track() and trackToFace() functions. The global variable "position_" within the particle class will be set based on former evaluations. You may also have a look at the hitWallFaces() function within the same file. If a tri is intersected, the "position_" value will possibly be set here.

Does that help you?

Best regards,
Andreas

[Atzori]

First, thank you so much for your reply

I have allready check, but, due to the name "endPosition", it seems to me that "somewhere else" I should found the evaluation of the the physical position, while the track function (and similar) have only a "secondary" aim - to match the pysical position with the mesh and the boundary constrain and so on...

Foam::label Foam::particle::track(const vector& endPosition, TrackData& td)

Also "when" this function is called is a bit misterious for me...

I attached an incomplete "flowchart" (it is a work in progess: "TO DO" means only "to write", I check), but where you can read "move" I do not manage to see something directly related with the position...

Thanks again for your help!

[A_Pete]

I don't really get what you mean with "endPosition". I am currently working in the same class, so I can try to elaborate on my thoughts:

position_ is definately the current position of the particle as can be seen in particle.H. endPosition is the position that the particle is supposed to be moved to within the given time step. This does now depend on the particle crossing a face, which is evaluated in hitWallFaces(). If one of these conditions is met, position_ is changed accordingly. If none of these conditions is met, position_ will be set to "endPosition" e.g.:

// Did not hit any tet tri faces, and no wall face has been
// found to hit.
if (tris.empty() && faceI_ < 0)
{
position_ = endPosition;
}

In KinematicParcel.C::move() the velocity U_ of a particle is set and the stepFraction_ accordingly, which IMO is based on some possibility of a particle crossing a face during this time step at the given velocity. So this does just set the time step fraction as well as the velocity, but the position_ IMO is really set in track() and trackToFace().

[Atzori]

I perfectly agree with you, but my question is, when you said:

"endPosition is the position that the particle is supposed to be moved"

endPosition where is evaluated?

Thanks again!

[A_Pete]

Ah, now I got you.

There seem to be two spots. The one you are looking for is probably this one in KinematicParcel.H (line 330):

dt1 *= dCorr/d*p.trackToFace(p.position() + dCorr*U_/magU, td);

This is also kind of confusing for me, since dt1 is set here. Though, by calling p.tracktoFace() the position_ variable should automatically be operated on, since that is what happens in trackToFace(). So the first argument "p.position() + dCorr*U_/magU" should be the local "endPosition" variable in trackToFace(). But there is this other spot in the track() function in particleTemplates.C:

stepFraction_ += trackToFace(endPosition, td)*(1.0 - stepFraction_);

You should probably be looking at the first spot, but that is just a guess.

Good luck!

[Atzori]

Dear Pete,
Forgive my lack of clarity: due to the desire to find the solution, I was answering during a meeting

I think (and it seems to me that you agree) that:
1) OF evaluates the new velocity (this part is clear for me);
2) OF evaluates the “endPosition”, (let call it the “theoretical position” or the “physical” one);
3) using several functions lead by “track”, OF decides the “practical” position, that depends on several things, as you said (hit a face? If yes, is it a boundary? If no, is it a subdivision between parallel domains? and so on).

Now, lines such as:

if (tris.empty() && faceI_ < 0) 
{ 
position_ = endPosition; 
}

Represent one of the possibile result of the process (3).

What I’d like to do is to identify where the “physical” new position is decided, so I’m looking for something as:

  endPosition = position_ + U_ * Delta_t;

That in my opinion should be “the only possibile end” of the (phantom ) process (2).

In other words, where is the code corrisponding to this description?
"endPosition is the position that the particle is supposed to be moved to within the given time step."

Cheers!!
Marco Atzori

[Atzori]

Sorry XDD

Thanks for the "Good luck"

[sisetrun]

Hello Marco, hello Andreas,

I found your thread while I was searching for kinematikcloud.move()

For my current task which I am working on, I want to do particle tracking where molecular diffusion is taken into account (solver: icoUncoupledKinematicParcelFoam). So the equation for the particle movement should look similar to this equation.
x(x+dt) = x(t) + u(x(t))*dt+f(D_mol)
with f(D_mol) = rand(vector)*sqrt(2*Dmol*/dt), which looks quiet similar to Marco's equation.
Have you made any progress in implementing the random"stuff"?

Since I have not developed any OF code yet (trying and learning at the moment), it is quit hard to follow the the steps in the code...
Thanks a lot for your help and info
Best regards,
Sebastian

[Atzori]

Hi sebastian!

I managed but in a "barbarian way":

The problem is that OpenFOAM has to check several things between when a new position is setted and when it is finally written, so a "final position" does not exist.

In few words: I use my equations (which indeed are very similar to yours ) to calc x(t+dt), then I calc a sort of "fake velocity":

U_ = x(t+dt)-x(t) / dt

So that U_ can be used as a "classical" velocity for the solver (where "classical" means "not stochastic").

(in this case it could be better to add another variable to the particle class to save the physical velocity)

I hope that it helps, do you know how to implement what I described or you need help?

Cheers!

Marco

(I'm very busy in this precise moment because of my thesis, If needed I'll post a more complete answer as soon as possible, likely in the evening)

[sisetrun]

Thanks Marco,

well, your "barbarian way" this is still a very clever idea to solve the problem !

Since I just started programming OF (basic C++, good python knowledge) any hint/help is thankfully appreciated !

At the moment, I am looking for a starting point. I am still a bit confused with all the classes and .H files in OpenFOAM and getting them linked together for a complete application!

Thanks a lot for your time

[Atzori]

Hi!

Have a look at the attached pdf and try to use the links in the figures (at the end, I thought that the easiest way would be to share one of my internal report: forgive the style pedantic a bit )

If something is not clear in the pdf ask me (anyway, do not consider all the details and be aware of possible errors)

Cheers!!

Marco

[sisetrun]

Thank you so much...
I will try to work myself through the code!

Have a nice weekend!

Cheers back

[openfoammaofnepo]

Hi guys,

In the srouce files for the lagrangian solver in OpenFOAM, I found that there is an important class called "TrackData", such as in the the following source files, it is used:

https://github.com/OpenFOAM/OpenFOAM...actingParcel.C

However, I did not find the corresponding source file for TrackData, such as TrackData.[H, C] in OpenFOAM. This is important since it is used by lots of the functions in the OF lagrangian solvers.

Could you please give some hints about where "TrackData" is defined? Thank you.

best regards, OFFO

Originally Posted by Atzori

Dear all,
I’d like to implement a particles stochastic model which requires to add random generators in the integration both of U_particle and of X_particle.
I though that the best option could be to find where U_p and X_p are updated to insert there my equations, but I didn't managed to find where the position of the single parcel is updated...
I shortly summed up what I considered the most important passages in the "flow" of the solver:

1. kinematicCloud.evolve() (called in the solver)
2. TrackingData created (called in evolve() in KinematicCloud.C)
3. solve(td) (called in evolve() in KinematicCloud.C)
4. preEvolve (called in solve() in KinematicCloud.C)
5. evolveCloud(td) (called in solve() in KinematicCloud.C)
6. td.cloud().motion(td); (called in evolveCloud(td) in KinematicCloud.C)
7. CloudType::move(td, solution_.trackTime()); (called in motion(td) in KinematicCloud.C)
8. bool keepParticle = p.move(td, trackTime); (called in move() in Cloud.C)
9. p.calc(td, dt, celli); (called in move() in KinematicParcel.C)
10. this->U_ = calcVelocity(td, dt, celli, Re, muc_, mass0, Su, dUTrans, Spu); (called in move() in KinematicParcel.C) (here is called the integrator for U, where I plan to set the new equation for U_p)
11. ...
12. ...
    

In the following I hoped to find something equivalent to:

new_X = current_X+ Delta_t * new_U

Where I can find it?

It seems to me that "something happens" in the TrackingData, but I did not managed to divide the "physical part" (i.e.: the new position) and "service part" (interaction with boundaries, change of processor and so on...)

On addition: do you have any suggestions for a better approach in coding?
(e.g.: I think that, if only U_p had to be changed, the best choice would have been to create a new force model, but it seems to me that it cannot help with X_p)
Please, consider that the "strategic task" is to obtain somethig as (in a very very very idealized picture):

X_p(n+1)=X_p(n) + \Delta t * U(n) + Random1
U_p(n+1)=(Drag term) + Random2


Thanks in advance!!

[Atzori]

Hi!

Unluckily, I cannot give a complete answer... If you are "just" trying to understand how it works, maybe it's enough to have a look at the KinematicCloud class:

https://github.com/OpenFOAM/OpenFOAM...KinematicCloud

I guess a couple of the methods of the TrackData class have been specified here -- at least, as it concerns the particles classes. (e.g.: see line 246 of KinematicCloud.H and at line 756 of KinematicCloud.C is explained the meaning of "trackFraction").

Bests!!

Marco

(P.S.: I'm really interested in understanding the Lagrangian tracking in OF: we can discuss about it in the limits of my experience.)

[openfoammaofnepo]

Dear Marco,

Thank you so much for your reply. I check the lines you pointed out, but I did not figure them out about where "TrackData" is defined. I will think about it and will let you know.

OffO

[DustExplosion]

Hello All,

Thank you for a great thread so far. Attached is my spread sheet of the class hierarchy in the Lagrange solver. Hopefully, you find it useful.

I am trying to verify the submodels independently and for this I am simulating a single parcel (Np = 1) falling in a tube. The velocity and displacement results are shown here:

displacement.png
velocity.png

As you can see gravity is increasing the "velocity" that is output (in time/lagrange/U); however, the position is not being updated. It also looks like the drag model is not getting an updated velocity either as terminal velocity does not occur. Note that the lines in the plots are the analytical answer.

Does anyone see what might be going wrong from an input standpoint? I am using the coalChermistySolver with OF 3.0.1. Just diving into the code now.

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      limestoneCloud1Properties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solution
{
    active          true;
    coupled         true;
    transient       yes;
    cellValueSourceCorrection on;
    maxCo           0.3;

    sourceTerms
    {
        schemes
        {
            U               explicit 1;
            h               explicit 1;
            radiation       explicit 1;
        }
    }

    interpolationSchemes
    {
        rho             cell;
        thermo:mu       cell;
        U               cellPoint;
        Cp              cell;
        kappa           cell;
        T               cell;
        G               cell;
    }

    integrationSchemes
    {
        U               Euler;
        T               analytical;
    }
}

constantProperties
{
    parcelTypeId    2;

    rho0            2500;
    T0              300;
    Cp0             900;

    epsilon0        1;
    f0              0.5;
}

subModels
{
    particleForces
    {
        sphereDrag;
    gravity;
    }

    injectionModels
    {
        model1
        {
            type            manualInjection;
            massTotal       0.0;
            //parcelBasisType mass;
        parcelBasisType fixed;
        nParticle 1;
            SOI             0;
            positionsFile   "limestonePositions";
            U0              (0 0 0);
            sizeDistribution
            {
                type        fixedValue;
                fixedValueDistribution
                {
                    value        10e-06;
                }
            }
        }
    }

    dispersionModel stochasticDispersionRAS;

    patchInteractionModel standardWallInteraction;

    heatTransferModel RanzMarshall;

    stochasticCollisionModel none;

    surfaceFilmModel none;

    radiation       off; // on;

    standardWallInteractionCoeffs
    {
        type            rebound;
        e               1;
        mu              0;
    }

    RanzMarshallCoeffs
    {
        BirdCorrection  false;
    }
}


cloudFunctions
{}


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

[DustExplosion]

Attached is a list of lagrange sub models and which class they "belong to" FYI.

[Atzori]

Hi Chris!

Thank you so much for sharing your work!

Unluckily, I'm not able to see the problem in your case: actually, I'm finding it awkward a bit and I'll try to reproduce it a similar case as soon as possible.

Have you used a uniform 0 velocity in all the "tube", right?

Let stay in touch to solve the problem!

See you soon!

Marco

(The only "strange" thing that I found in your input file is the fact that you turn on the "dispersion model", which in my personal opinion is not appropriate considering your aim... but I do not see how it can be related with the "infinity" terminal velocity)

[DustExplosion]

Hi Atzori,

Good eye on the dispersion model, but I cannot seem to remove it (the solver complains). I do not have turbulence on so I do no think it is doing anything.

In the case I posted, the velocity is initially zero but increases during the simulation. This is why the drag "goes away" However, I am not sure why the fluid velocity increases as the particle mass is only 1.308996939e-12 (my domain is 3x3x3 cm which I guess is small but the volume fraction of the particle is still much smaller!).

Regardless, I changed the problem bit to make sure that boundary conditions were not the issue. Now I have a horizontal tube with a constant inflow velocity of 1 m/s and no gravity. The particle still seems to effect the flow a bit which I do not understand, again because the mass is so small. The particle acceleration does not look right either, compared to the analytical solution.

Displacement_DragOnly.png
Velocity_DragOnly.png

Attached is the case if you want to take a look. If you have python installed you should be able to change the solver in batchRun and run using:

./batchRun > console.out
python velocity.py

These commands will give you the plots above.

Looking forward to figuring this out together!