[ansubru]

Hey everybody!!

I am new to openfoam and I have been having a lot of difficulty in trying to understand and implement the drag formulation. I am trying to simulate a system of particles (dense packing, vol fraction 0.85) in a box. This simulation is an extension of the hopper tutorial (icoUncoupledKinematicParcelFoam).

For a given system, particle gravity is constant but gas drag changes dramatically (in a strong non-linear way) with the local particle number density. The presence of other particles restricts fluid space, creates a sharp velocity gradient in the surrounding gas phase and, as a result, yields increased shear stress acting on the particle surface. This is the effect i want to quantify.

There are several relations in literature to account for expressing drag in the presence of other particles. i would like to implement one/several of these. I am having a hard time trying to figure out in which file/s are the expression for drag force being estimated (I have examined the file Spheredragforce.C in the location ~./src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag). But this file only provides the expression for the drag coe-eff Cd, what i want is the expression for the drag force Fdrag, which I am unable to locate.

Could anyone help me in identifying this source file??

Regards

ansubru

[wyldckat]

Hi ansubru,

As promised via private message, I've taken a better look into this.

I've taken a quick look at the source code and I think that the first step is to determine at what point do you need to perform this calculation.
I ask this, because from your description, what you seem to want to do is a post-processing feature, not something that is done for the solution of the simulation itself.

Best regards,
Bruno

[ansubru]

Hi!!

Thank you for looking into this!!

I am looking at this change more at the "source".. as in.. I would like to change the way drag is calculated by including a volume fraction function (f(epsilon)). This would mean that the drag calculation at every time step must include this function.

Does this help???

[wyldckat]

Hi ansubru,

I've finally managed to give a better look into this.

If you have a look into the classes defined in the folder "src/lagrangian/intermediate/submodels/Kinematic/ParticleForces", there is an interesting pattern:

• The sub-folder "ParticleForce" has the base class for all possible forces. It has 3 important methods:
  

               //- Calculate the coupled force
              virtual forceSuSp calcCoupled
              (
                  const typename CloudType::parcelType& p,
                  const scalar dt,
                  const scalar mass,
                  const scalar Re,
                  const scalar muc
              ) const;
  
              //- Calculate the non-coupled force
              virtual forceSuSp calcNonCoupled
              (
                  const typename CloudType::parcelType& p,
                  const scalar dt,
                  const scalar mass,
                  const scalar Re,
                  const scalar muc
              ) const;
  
              //- Return the added mass
              virtual scalar massAdd
              (
                  const typename CloudType::parcelType& p,
                  const scalar mass
              ) const;

• All rely in the class "forceSuSp", located in the sub-folder of the same name, which has the following description:
  

      Helper container for force Su and Sp terms.
  
          F = Sp(U - Up) + Su
  
      Explicit contribution, Su specified as a force
      Implicit coefficient, Sp specified as force/velocity

• For example, the sub-folder "Lift/LiftForce" has the base class structure for lift, which relies in the method "calcCoupled" and defines only the "Su" part and "Sp" is set to "0.0".
• The sub-folder "Drag" has the classes for drag, from which you've pointed out "SphereDrag". These classes also rely in the method "calcCoupled" and defines only the "Sp" part.
• The sub-folder "Gravity" has a single class that relies in the method "calcNonCoupled" and affects only "Su".
• The class in the sub-folder "PressureGradient" also relies in the method "calcCoupled" and also affects "Su" only.

OK, so this gets us the individual accountabilities of each force component. The remaining question is: how are all of them added up?
It's in the class "KinematicParcel", method "calcVelocity", located in the folder "src/lagrangian/intermediate/parcels/Templates/KinematicParcel/". It relies on the class "TrackData::cloudType::forceType", where "TrackData" is a template name... OK, I got lost here, because this section is deeply based in templates, which is a considerable pain to figure out who-is-what-and-where-and-how
But the point is that this line:

const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);

is probably the addition of all forces that is listed inside "forces".

And the rest... is up to you!

Best regards,
Bruno

[Sören Sander]

Hi ansu, Bruno,

I am also investigating this solver. I still have many remaining questions, but maybe the follwing gives you a starting point:

A. For information on template for class Trackdata look into /intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C.

B. If you want to include a new force to the solver, you may try the following (which I haven't tried yet):

• Copy the intermediate library, recompile it to something like "mylib" in your user folder and link the new user library to your own copy of the solver (myIco... )
• in submodels/Kinematic/ParticleForces/.. copy a drag model and rename it (do not forget to rename the classes inside the file)
• In intermediate/parcels/include/makeParcelForces.H include the file "myNewForce.H" and add the force to the definition (register it to the library):
  

  ...
  #include "myNewForce.H"
  
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  
  #define makeParcelForces(CloudType)                                           \
                                                                                \
      makeParticleForceModel(CloudType);                                        \
      makeParticleForceModelType(myNewForce, CloudType);                   \
  ...

• In the case file /constant/kinematicCloudProperties add
  
  

  ...
  subModels
  {
      particleForces
      {
          sphereDrag;
          gravity;
          myNewForce;
      }
  ...

I might find some time at the weekend to implement a new force and check if the works. I will keep you updated.

Best regards
Sören

[ansubru]

Hi Soren!!

Thank you so much for that update. Your views have indeed helped me. I have compiled my own user libraries and am modifying the source code there. That idea was really a very profitable one.

In addition to that, I have encountered another huge error during my computation. I have started a thread about the same, could you give me your impressions about these?? The thread in question is --

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

Best regards

ansubru

[Sören Sander]

I implemented a new force during the weekend and it seems to work without bigger problems. The only thing, which I could not figure out, is how additional fields of the carrier phase are implemented to the library (e.g. the volVectorField U). Any hints on this?

All I found was

template <class ParcelType>
inline const Foam::vector& Foam::KinematicParcel<ParcelType>::Uc() const
{
    return Uc_;
}

[sfigato]

Hi Soren and Ananda,

I am trying to implement a new drag model as well. Can you explain better which steps must I do after I compiled the new library??

I try to explain myself better. After I copied the Myintermediate lib in my user library folder, how can I link this one with the solver myIcoUncoupledKinematicParcelFoam. What must I change in the file Make/option of the solver?

Thanks in advance for the answer
Regards
Marco

[wyldckat]

Greetings to all!

@Marco: It really depends on what exactly you're trying to achieve, and knowing which OpenFOAM version you're using would also help .

Using OpenFOAM 2.3 as a reference: having a new library "lagrangianIntermediate" means that you'll need to replace in "Make/options" the following lines:

• Replace this:
  

  -llagrangianIntermediate \

  with this:
  

  -L$(FOAM_USER_LIBBIN) -lmylagrangianIntermediate \

• And this:
  

  -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \

  with this:
  

  -I$(WM_PROJECT_USER_DIR)/myintermediate/lnInclude \

  To know where the variable "WM_PROJECT_USER_DIR" points to, run this command:
  

  echo $WM_PROJECT_USER_DIR

But you might trigger a problem, where another library tried to load up the original "lagragianIntermediate", leading to a collision of objects of the same name, when running the solver.


So, my suggestion is that:

1. Try to adapt the instructions I wrote above.
2. If the solver crashes due to an error message about something like "link ... duplicate object", then try to diagnose what exact changes you really need to do. Because it might be possible to simply create a new drag model, without having to have an entire copy of the original library source code.
   For example, it might be as simple as creating a new boundary condition, as exemplified here: http://www.cfd-online.com/Forums/ope...tml#post446451 - post #10

Best regards,
Bruno

[sfigato]

Hi Bruno,

I thank you for the answer. The OpenFoam version is 2.3.0. I try to follow your suggestions.
Anyway, I would you like to recompile the whole library because I want to compute as well the volume fraction of the lagrangian phase.

Regarding this, do you have some suggestion on in which part of the code is available the number of particles, their diameter and the mesh volume? So which part of the code I need to modify or at least have look?

Thanks in advance
Best Regards
Marco

[ansubru]

Hi Guys!!

I have compiled my own dynamic linked library from scratch.. Most of the steps needed are available in previous posts, just as Bruno suggested ... Isnt the volume fraction of the carrier phase already computed (look at alpha.air)?? Maybe that can help??

Anyways, I had some questions for all you guys.. I want to modify the DPMFoam solver so that it can handle dynamic mesh computations (I need to introduce a moving wall at one of the boundaries and from most of the tutorials pimpleDyMFoam test tube mixer vessel) .I have deduced that the best way to go about this is by using a 'DyM' approach.. My questions are as follows -

1) Can the DPMFoam solver handle 'DyM* cases.. I guess this is possible as the DPMFoam solver is nothing but a combination of pimpleFoam and icoUncoupledKinematicFoam ( both of these solvers have DyM capability)

2) How could I go about with this.. As of now I am trying to adapt the DPMFoam code pimple velocity loop to the corresponding loop from pimpleDyMFoam (i.e. replaced the loop in the original solver with this one.. ) Is this the correct way forward..?? Can somebody help me with this??

I will paste the DPMDymFoam solver script here and highlight the bits I have altered,

Application
    DPMDyMFoam

Description
    Transient solver (incorporating moving meshes) for the coupled transport of a single kinematic particle cloud including the effect of the volume fraction of particles on the continuous phase.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "singlePhaseTransportModel.H"
#include "PhaseIncompressibleTurbulenceModel.H"
#include "pimpleControl.H"
#include "fixedFluxPressureFvPatchScalarField.H"
#include "fvIOoptionList.H"


#ifdef MPPIC
    #include "basicKinematicMPPICCloud.H"
    #define basicKinematicTypeCloud basicKinematicMPPICCloud
#else
    #include "basicKinematicCollidingCloud.H"
    #define basicKinematicTypeCloud basicKinematicCollidingCloud
#endif

int main(int argc, char *argv[])
{
    argList::addOption
    (
        "cloudName",
        "name",
        "specify alternative cloud name. default is 'kinematicCloud'"
    );

    #include "setRootCase.H"
    #include "createTime.H"
    #include "createDynamicFvMesh.H"
    #include "readGravitationalAcceleration.H"
    #include "createFields.H"
    #include "createUf.H"
    #include "createFvOptions.H"
    #include "initContinuityErrs.H"

    pimpleControl pimple(mesh);

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

    while (runTime.run())
    {
        
	#include "readControls.H"
        #include "CourantNo.H"
        #include "setDeltaT.H"
        #include "createPcorrTypes.H"
        runTime++;

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

        mesh.update();

        // Calculate absolute flux from the mapped surface velocity
        phi = mesh.Sf() & Uf;

        if (mesh.changing() && correctPhi)
        {
            #include "correctPhi.H"
        }

        // Make the flux relative to the mesh motion
        fvc::makeRelative(phi, U);

        if (mesh.changing() && checkMeshCourantNo)
        {
            #include "meshCourantNo.H"
        }

        U.correctBoundaryConditions();

        continuousPhaseTransport.correct();
        muc = rhoc*continuousPhaseTransport.nu();

        Info<< "Evolving " << kinematicCloud.name() << endl;
        kinematicCloud.evolve();

        // Update continuous phase volume fraction field
        alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
        alphac.correctBoundaryConditions();
        alphacf = fvc::interpolate(alphac);
        alphaPhic = alphacf*phic;

        fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));
        volVectorField cloudVolSUSu
        (
            IOobject
            (
                "cloudVolSUSu",
                runTime.timeName(),
                mesh
            ),
            mesh,
            dimensionedVector
            (
                "0",
                cloudSU.dimensions()/dimVolume,
                vector::zero
            ),
            zeroGradientFvPatchVectorField::typeName
        );

        cloudVolSUSu.internalField() = -cloudSU.source()/mesh.V();
        cloudVolSUSu.correctBoundaryConditions();
        cloudSU.source() = vector::zero;

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UcEqn.H"

            // --- PISO loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

            if (pimple.turbCorr())
            {
                continuousPhaseTurbulence->correct();
            }
        }

        runTime.write();

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

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

    return 0;
}


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

I have included the correctPhi.H file from the PimpleDym Solver??

Could anyone go through and tell me if this methodology is correct?? Moreover I get the following errors during compilation of this solver -

In file included from /home/ask/OpenFOAM/OpenFOAM-2.3.0/src/fvOptions/lnInclude/fvOptionList.H:41:0,
                 from /home/ask/OpenFOAM/OpenFOAM-2.3.0/src/fvOptions/lnInclude/fvIOoptionList.H:38,
                 from DPMDyMFoam_mod.C:40:
/home/ask/OpenFOAM/OpenFOAM-2.3.0/src/fvOptions/lnInclude/fvOption.H:53:24: fatal error: meshToMesh.H: No such file or directory
compilation terminated.
make: *** [Make/linux64GccDPDebug/DPMDyMFoam_mod.o] Error 1

I guess the fvOptions.H is not defined within this scope (although I have included it as a header just as in the pimpleDymfoam solver)?? Could someone tell me how this can be fixed??

Thanks guys

Regards

ansubru

[ansubru]

Hi Again Guys!!

I tried to remove some of of the redundant headers in the above code and now i get this error while compilation -

In file included from DPMDyMFoam_mod.C:65:0:
/home/ask/OpenFOAM/OpenFOAM-2.3.0/src/finiteVolume/lnInclude/createUf.H: In function ‘int main(int, char**)’:
/home/ask/OpenFOAM/OpenFOAM-2.3.0/src/finiteVolume/lnInclude/createUf.H:49:22: error: ‘U’ was not declared in this scope
/home/ask/OpenFOAM/OpenFOAM-2.3.0/src/finiteVolume/lnInclude/createUf.H:54:14: error: ‘phi’ was not declared in this scope
In file included from DPMDyMFoam_mod.C:66:0:
/home/ask/OpenFOAM/OpenFOAM-2.3.0/src/fvOptions/lnInclude/createFvOptions.H:1:1: error: ‘IOoptionList’ is not a member of ‘Foam::fv’
/home/ask/OpenFOAM/OpenFOAM-2.3.0/src/fvOptions/lnInclude/createFvOptions.H:1:18: error: expected ‘;’ before ‘fvOptions’
DPMDyMFoam_mod.C:87:9: error: ‘phi’ was not declared in this scope
make: *** [Make/linux64GccDPDebug/DPMDyMFoam_mod.o] Error 1

Clearly this is not the way forward as both U and phi are not defined in the scope of this solver.. could someone then tell me how can I introduce this dynamic capability??

Thanks again

Regards

ansubru

[wyldckat]

Greetings to all!
And sorry for the reaaaaally late reply, but I'm still catching up on threads I wasn't able to answer since July.

@Marco:

Originally Posted by sfigato

Regarding this, do you have some suggestion on in which part of the code is available the number of particles, their diameter and the mesh volume? So which part of the code I need to modify or at least have look?

I don't know if you've found the answer to this or not?

------------------

@ansubru: I've searched for other posts you've made since the post above and I believe you've gotten an answer in this post: http://www.cfd-online.com/Forums/ope...tml#post501356 post #3 ?

Best regards,
Bruno

[sfigato]

Hi Bruno,

@Bruno

I don't know if you've found the answer to this or not?

thanks for the answer. No, I did not find any suggestions on it.

I am using uncoupledKinematicParcelFoam. There is just a deterministic collision model and I tryed to implement (copy) there a sthocastic collision model, which accounts as well for the coalescence of the particles. I copied the O'Rourke collision model and when i tryed to compile it I got some error messages due the lack of the reaction and thermo partin uncoupledKinematicParcelFoam. Dou you know how I can implement and compile a sthocastic collision model with coalescence for uncoupledKinematicParcelFoam?

Best Regards
Marco Longhitano

[wyldckat]

Hi Marco,

OK, I did a quick search and based myself on some stuff I already know from my studies for this thread: http://www.cfd-online.com/Forums/ope...mppicfoam.html
In other words, I'm not an expert on this topic, I only know some bits and pieces, which hopefully can guide you in the right direction.

Originally Posted by sfigato

which part of the code is available the number of particles, their diameter and the mesh volume?

• Number of particles is managed by the Cloud class: "src/lagrangian/basic/Cloud/"
• Diameter: It's not contemplated in the model directly as an exact diameter. Closest there is are densities and maximum distances. It's a bit of a confusing code, so I strongly advise you to study C++ templates in more detail, before looking directly at OpenFOAM's source code. Either way, files of importance on this topic of dimensions of the particle:
  

  src/lagrangian/basic/particle/
  src/lagrangian/intermediate/parcels/Templates/KinematicParcel/
  src/lagrangian/intermediate/clouds/derived/basicKinematicCloud/

  From the last file on this folder list you can see that it uses the other two classes.
  • Essentially, this line from the last file:
    

    typedef KinematicCloud<Cloud<basicKinematicParcel> > basicKinematicCloud;

    defines:
    • KinematicCloud - the main cloud interaction model
      • Cloud - the cloud manager
        • basicKinematicParcel - the parcel/particle, which is in turn divided into "KinematicParcel<particle>", where "KinematicParcel" manages the particles' own existence, as an entity with density and mass.
• Mesh volume: not directly addressed, as far as I can understand. It might be possible to infer the proportionally occupied volume, based on the total mass, perhaps?
  Either way, the closest is to look at the source code for the "cloudInfo" function object, at:
  

  src/postProcessing/functionObjects/cloud/cloudInfo/

Originally Posted by sfigato

There is just a deterministic collision model and I tryed to implement (copy) there a sthocastic collision model, which accounts as well for the coalescence of the particles. I copied the O'Rourke collision model and when i tryed to compile it I got some error messages due the lack of the reaction and thermo partin uncoupledKinematicParcelFoam. Dou you know how I can implement and compile a sthocastic collision model with coalescence for uncoupledKinematicParcelFoam?

I feel like I can't process all of those expensive words Sorry, it's just since I'm not very familiar with this topic, there are a lot of terms that feel like they go straight past me

Nonetheless, if I can understand you correctly, you want to add features to a collision model.
First a quick overview:

1. You need to keep in mind that OpenFOAM might not fully support the kind of modelling you want to do, or at least maybe not directly.
2. I believe that the cloud management system comes some a long time ago, mostly used for DSMC, which was introduced in OpenFOAM 1.6: http://www.openfoam.org/archive/1.6/...ease-notes.php - search for the description for dsmcFoam shown there.
3. Read-through the thread I mentioned before, for more ideas on the topic of how OpenFOAM handles clouds and particles: http://www.cfd-online.com/Forums/ope...mppicfoam.html

Now, if I'm not mistaken, you will need to create a new cloud model. You'll need to study the following folders for more ideas:

src/lagrangian/intermediate/clouds/baseClasses/
src/lagrangian/intermediate/clouds/derived/
src/lagrangian/intermediate/clouds/Templates/

On the other hand, a quick look at the file "basicReactingCloud.H" located at "src/lagrangian/intermediate/clouds/derived/basicReactingCloud/", makes me believe that all you need to do to change this block of code:

    basicKinematicCloud kinematicCloud
    (
        kinematicCloudName,
        rho,
        U,
        thermo.mu(),
        g
    );

by changing "basicKinematicCloud" to "basicReactingCloud". For ideas on how to do this, look at the source code for the solver reactingParcelFilmFoam.


Beyond this, I'm all out of ideas, at least for now.

Best regards,
Bruno

[sfigato]

Hi Bruno,

I thank you for the answer. Your suggestions will be really helpful. I will go ahead with my implementation and let the thread update.

I would like also to ask to you an opinion about another post that I recently opened. Here is the link:

http://www.cfd-online.com/Forums/ope...tml#post510287

When/If you have time please have a look.

Thanks
Regards

Marco