CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Pre-Processing

SaffmanMeiLiftForce

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

Like Tree6Likes
  • 2 Post By wyldckat
  • 1 Post By maysmech
  • 2 Post By wyldckat
  • 1 Post By lithos

Reply
 
LinkBack Thread Tools Display Modes
Old   October 13, 2013, 11:16
Default SaffmanMeiLiftForce
  #1
Member
 
Evangelos
Join Date: Sep 2011
Posts: 67
Rep Power: 5
Danath is on a distinguished road
Hi

Can anyone provide an example of SaffmanMeiLiftForce dictionary ?

thanks a lot !
Danath is offline   Reply With Quote

Old   October 13, 2013, 14:51
Default
  #2
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,253
Blog Entries: 34
Rep Power: 84
wyldckat is just really nicewyldckat is just really nicewyldckat is just really nicewyldckat is just really nice
Hi Evangelos,

I had no idea there was such a thing as "SaffmanMeiLiftForce" in OpenFOAM...
But if you can provide an example case, something simple based on one of OpenFOAM's tutorials, in which this dictionary can be used, then it should be pretty simple to reverse engineer it!
Actually, OpenFOAM helps with this in two ways:
Best regards,
Bruno
wyldckat is offline   Reply With Quote

Old   October 13, 2013, 16:00
Default
  #3
Member
 
Evangelos
Join Date: Sep 2011
Posts: 67
Rep Power: 5
Danath is on a distinguished road
Hello Bruno !

I try to calculate the orbit of a particle when injected into a domain without fluid motion.I use icoUncoupledKinematicFoam but i am a little bit confused!
All the solvers take into account dragforce and gravity without calculate the lift from the fluid !

In openfoam 2.2.0 particle lift is available

So in the kinematicTransportproperties/particleforces under the drag force and gravity i add particle force and the program show me the valid particle forces

There are two " lift methods " TomiyamaLift and SaffmanMeiLiftForce

About the dimension of Saffman-Mei particle lift force model
Danath is offline   Reply With Quote

Old   October 14, 2013, 04:06
Default
  #4
Member
 
Evangelos
Join Date: Sep 2011
Posts: 67
Rep Power: 5
Danath is on a distinguished road
I need something like that

pressureGradient dictionary

but for saffmanmeiliftforce!!!
Danath is offline   Reply With Quote

Old   October 14, 2013, 16:47
Default
  #5
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,253
Blog Entries: 34
Rep Power: 84
wyldckat is just really nicewyldckat is just really nicewyldckat is just really nicewyldckat is just really nice
Hi Evangelos,

Do you mean "icoUncoupledKinematicParcelFoam"?

Do you think that the tutorial "lagrangian/icoUncoupledKinematicParcelFoam/hopper" can be used to test this dictionary you are asking about? If so, I can give it a try next weekend.

Best regards,
Bruno
wyldckat is offline   Reply With Quote

Old   October 4, 2014, 23:34
Question
  #6
Senior Member
 
maysmech's Avatar
 
Join Date: Jan 2010
Posts: 347
Blog Entries: 2
Rep Power: 8
maysmech is on a distinguished road
Hi all,
Do you find any dictionary for saffmanMeiLiftForce?
maysmech is offline   Reply With Quote

Old   October 5, 2014, 04:20
Default
  #7
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,253
Blog Entries: 34
Rep Power: 84
wyldckat is just really nicewyldckat is just really nicewyldckat is just really nicewyldckat is just really nice
Greetings Maysam,

Quote:
Originally Posted by maysmech View Post
Do you find any dictionary for saffmanMeiLiftForce?
At least I haven't, because Evangelos didn't answer back

Therefore my questions and offer still stand:
Quote:
Originally Posted by wyldckat View Post
Do you mean "icoUncoupledKinematicParcelFoam"?

Do you think that the tutorial "lagrangian/icoUncoupledKinematicParcelFoam/hopper" can be used to test this dictionary you are asking about? If so, I can give it a try next weekend.
Best regards,
Bruno
wyldckat is offline   Reply With Quote

Old   October 5, 2014, 04:33
Question
  #8
Senior Member
 
maysmech's Avatar
 
Join Date: Jan 2010
Posts: 347
Blog Entries: 2
Rep Power: 8
maysmech is on a distinguished road
For example for DPMFoam
maysmech is offline   Reply With Quote

Old   October 5, 2014, 06:58
Default
  #9
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,253
Blog Entries: 34
Rep Power: 84
wyldckat is just really nicewyldckat is just really nicewyldckat is just really nicewyldckat is just really nice
OK, so as I mentioned in one of the posts above, having an example case as a basis and using the "banana" trick, here's what I've done (used OpenFOAM 2.3.x for this example):
  1. Made a copy of the tutorial "lagrangian/DPMFoam/Goldschmidt" and ran the utility blockMesh.
  2. Then I had to look at which file needed to modified. All pointed to the file "constant/kinematicCloudProperties".
  3. Now I had to figure out where the Lift comes in. After going around looking at the other tutorials and a few threads, the conclusion was that this block of dictionary is where we can stack up the various forces at work:
    Code:
        particleForces
        {
            ErgunWenYuDrag
            {
                alphac alpha.air;
            }
            gravity;
        }
  4. And here's where the "banana" comes in:
    Code:
        particleForces
        {
            ErgunWenYuDrag
            {
                alphac alpha.air;
            }
            gravity;
            banana;
        }
  5. As soon as I run DPMFoam, it complains about "banana" not being valid:
    Code:
    --> FOAM FATAL ERROR: 
    Unknown particle force type banana, constructor not in hash table
    
        Valid particle force types are:
    
    13
    (
    ErgunWenYuDrag
    PlessisMasliyahDrag
    SRF
    SaffmanMeiLiftForce
    TomiyamaLift
    WenYuDrag
    gravity
    nonInertialFrame
    nonSphereDrag
    paramagnetic
    pressureGradient
    sphereDrag
    virtualMass
    )
    There's our pretty lift model "SaffmanMeiLiftForce".
  6. OK, rename "banana" to "SaffmanMeiLiftForce" and added brackets with a "banana" inside :
    Code:
        particleForces
        {
            ErgunWenYuDrag
            {
                alphac alpha.air;
            }
            gravity;
            SaffmanMeiLiftForce
            {
                banana;
            }
        }
    mmm... the solver crashed. So much for the "banana".
  7. OK, if we look at the source code for this model:
    Code:
    template<class CloudType>
    Foam::SaffmanMeiLiftForce<CloudType>::SaffmanMeiLiftForce
    (
        CloudType& owner,
        const fvMesh& mesh,
        const dictionary& dict,
        const word& forceType
    )
    :
        LiftForce<CloudType>(owner, mesh, dict, forceType)
    {}
    there is no coefficient necessary.
  8. Therefore, this should be enough:
    Code:
        particleForces
        {
            ErgunWenYuDrag
            {
                alphac alpha.air;
            }
            gravity;
            SaffmanMeiLiftForce;
        }
  9. Er, nope... it complains it should be a dictionary itself. Then back to this:
    Code:
        particleForces
        {
            ErgunWenYuDrag
            {
                alphac alpha.air;
            }
            gravity;
            SaffmanMeiLiftForce
            {
            }
        }
  10. So we go back to the crash. The error given is this:
    Code:
    --> FOAM FATAL ERROR: 
    
        request for volVectorField U from objectRegistry region0 failed
        available objects of type volVectorField are
    1(U.air)
    OK, it's looking for the "U" field and it only found "U.air". So, either the model is just wrong, or something else is missing.
  11. Looking at the class from which this model derives, the constructor has a vital clue to what's missing:
    Code:
    template<class CloudType>
    Foam::LiftForce<CloudType>::LiftForce
    (
        CloudType& owner,
        const fvMesh& mesh,
        const dictionary& dict,
        const word& forceType
    )
    :
        ParticleForce<CloudType>(owner, mesh, dict, forceType, true),
        UName_(this->coeffs().template lookupOrDefault<word>("U", "U")),
        curlUcInterpPtr_(NULL)
    {}
    It's looking for a new name for "U", something like:
    Code:
    U U.air;
  12. So let's try this:
    Code:
        particleForces
        {
            ErgunWenYuDrag
            {
                alphac alpha.air;
            }
            gravity;
            SaffmanMeiLiftForce
            {
                U U.air;
            }
        }
    Side note: curiously enough, Maysam got a similar answer a few minutes ago in post #4 at pressureGradient dictionary
  13. Running DPMFoam once again gives this error message:
    Code:
    --> FOAM FATAL IO ERROR: 
    keyword curlUcDt is undefined in dictionary "/home/ofuser/OpenFOAM/ofuser-2.3.x/run/tutorials/lagrangian/DPMFoam/Goldschmidt/constant/kinematicCloudProperties.solution.interpolationSchemes"
    
    file: /home/ofuser/OpenFOAM/ofuser-2.3.x/run/tutorials/lagrangian/DPMFoam/Goldschmidt/constant/kinematicCloudProperties.solution.interpolationSchemes from line 27 to line 29.
    OK, we need to look at this this block of dictionary:
    Code:
    solution
    {
        active          true;
        coupled         true;
        transient       yes;
        cellValueSourceCorrection off;
    
        interpolationSchemes
        {
            rho.air         cell;
            U.air           cellPoint;
            mu.air          cell;
        }
    
        integrationSchemes
        {
            U               Euler;
        }
    
        sourceTerms
        {
            schemes
            {
                U semiImplicit 1;
            }
        }
    }
  14. Let's try:
    Code:
        interpolationSchemes
        {
            rho.air         cell;
            U.air           cellPoint;
            mu.air          cell;
            curlUcDt        banana;
        }
    The error message stated by DPMFoam:
    Code:
    --> FOAM FATAL ERROR: 
    Unknown interpolation type banana for field curlUcDt
    
    Valid interpolation types : 
    
    6
    (
    cell
    cellPatchConstrained
    cellPoint
    cellPointFace
    cellPointWallModified
    pointMVC
    )
  15. I have no idea. Let's try the most conservative/complex one "cellPointFace"... nope, it crashes. Let's try all others... nope, all of them crash.
  16. The error message of the crash is something like this:
    Code:
    Solving 3-D cloud kinematicCloud
    #0  Foam::error::printStack(Foam::Ostream&) in "/home/ofuser/OpenFOAM/OpenFOAM-2.3.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
    #1  Foam::sigFpe::sigHandler(int) in "/home/ofuser/OpenFOAM/OpenFOAM-2.3.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
    #2   in "/lib/x86_64-linux-gnu/libc.so.6"
    #3  Foam::SaffmanMeiLiftForce<Foam::KinematicCloud<Foam::Cloud<Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > > > >::Cl(Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > const&, Foam::Vector<double> const&, double, double) const in "/home/ofuser/OpenFOAM/OpenFOAM-2.3.x/platforms/linux64GccDPOpt/lib/liblagrangianIntermediate.so"
    #4  Foam::LiftForce<Foam::KinematicCloud<Foam::Cloud<Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > > > >::calcCoupled(Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > const&, double, double, double, double) const in "/home/ofuser/OpenFOAM/OpenFOAM-2.3.x/platforms/linux64GccDPOpt/lib/liblagrangianIntermediate.so"
    #5  
     in "/home/ofuser/OpenFOAM/OpenFOAM-2.3.x/platforms/linux64GccDPOpt/bin/DPMFoam"
    #6
    and it continues with more information, but doesn't help much more...
    So there was a "sigFpe", hence division by zero or infinite or something like that, in the method "::Cl" (it's in line #3).
  17. OK, looking at the source code for this method:
    Quote:
    https://github.com/OpenFOAM/OpenFOAM...MeiLiftForce.C
    Code:
    template<class CloudType>
    Foam::scalar Foam::SaffmanMeiLiftForce<CloudType>::SaffmanMeiLiftForce::Cl
    (
        const typename CloudType::parcelType& p,
        const vector& curlUc,
        const scalar Re,
        const scalar muc
    ) const
    {
        scalar Rew = p.rhoc()*mag(curlUc)*sqr(p.d())/(muc + ROOTVSMALL);
        scalar beta = 0.5*(Rew/(Re + ROOTVSMALL));
        scalar alpha = 0.3314*sqrt(beta);
        scalar f = (1.0 - alpha)*exp(-0.1*Re) + alpha;
    
        scalar Cld = 0.0;
        if (Re < 40)
        {
            Cld = 6.46*f;
        }
        else
        {
            Cld = 6.46*0.0524*sqrt(beta*Re);
        }
    
        return 3.0/(mathematical::twoPi*sqrt(Rew))*Cld;
    }
    there are a few possible locations where it might have broken.
  18. OK, let's go old-school debugging:
    Code:
    template<class CloudType>
    Foam::scalar Foam::SaffmanMeiLiftForce<CloudType>::SaffmanMeiLiftForce::Cl
    (
        const typename CloudType::parcelType& p,
        const vector& curlUc,
        const scalar Re,
        const scalar muc
    ) const
    {
        scalar Rew = p.rhoc()*mag(curlUc)*sqr(p.d())/(muc + ROOTVSMALL);
        Info << "Rew: " << Rew << endl;
        
        scalar beta = 0.5*(Rew/(Re + ROOTVSMALL));
        Info << "Re: " << Re << endl;
        Info << "beta: " << beta << endl;
    
        scalar alpha = 0.3314*sqrt(beta);
        Info << "alpha: " << alpha << endl;
    
        scalar f = (1.0 - alpha)*exp(-0.1*Re) + alpha;
        Info << "f: " << f << endl;
    
        scalar Cld = 0.0;
        if (Re < 40)
        {
            Cld = 6.46*f;
        }
        else
        {
            Cld = 6.46*0.0524*sqrt(beta*Re);
        }
    
        Info << "Cld: " << Cld << endl;
    
        return 3.0/(mathematical::twoPi*sqrt(Rew))*Cld;
    }
    we print out everything so that we can have a look at each value. It's good to know what's happening under-the-hood, instead of having to guess: http://openfoamwiki.net/index.php/HowTo_debugging
  19. And we have to build the modified library "src/lagrangian/intermediate":
    Code:
    wmake libso $FOAM_SRC/lagrangian/intermediate
  20. Let's have another go with DPMFoam... the output we now get:
    Code:
    Solving 3-D cloud kinematicCloud
    Rew: 0
    Re: 71.5702
    beta: 0
    alpha: 0
    f: 0.00077937
    Cld: 0
    that's waaaay too many zeros... this leads to this expression to go haywire:
    Code:
    return 3.0/(mathematical::twoPi*sqrt(Rew))*Cld;
    Because it becomes:
    Code:
    3.0/(2*Pi*sqrt(0.0))*0.0
    where it divides by zero.
  21. OK, we can try and insert a division protection like it's done for the others:
    Code:
    return 3.0/(mathematical::twoPi*sqrt(Rew+SMALL))*Cld;
  22. Build the library again and gave DPMFoam another run... didn't crash, but gave a lot of junk output. So I Ctrl+C to cancel the run and went back to the code and removed the "Info" lines I had added. And built the library again.
  23. Another run with DPMFoam and... and it's finally working! As to whether it's giving good or bad values, I have no idea.
In the meantime I'll report this bug, although I'm not sure if the fix I presented is correct or not .

-------------------
edit: Bug reported here: http://www.openfoam.org/mantisbt/view.php?id=1408

Additional note: Letting the case run first without any particles should initialize the "U.air" field, should no longer require you to modify the source code, since the flow field would all not be zero.
-------------------

Best regards,
Bruno
Mojtaba.a and matin like this.

Last edited by wyldckat; October 5, 2014 at 07:16. Reason: see "edit:" and "Additional note:"
wyldckat is offline   Reply With Quote

Old   October 5, 2014, 20:25
Default
  #10
Senior Member
 
maysmech's Avatar
 
Join Date: Jan 2010
Posts: 347
Blog Entries: 2
Rep Power: 8
maysmech is on a distinguished road
Thank you very much Bruno. I used a lot from your step by step debugging.
Especially Banana and Old-School debugging were fantastic.
This problem is for other forces too such as pressure gradient force.
However thanks for submit of the bug.
wyldckat likes this.
maysmech is offline   Reply With Quote

Old   October 18, 2014, 13:38
Default
  #11
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,253
Blog Entries: 34
Rep Power: 84
wyldckat is just really nicewyldckat is just really nicewyldckat is just really nicewyldckat is just really nice
FYI: Bug has been fixed: http://www.openfoam.org/mantisbt/view.php?id=1408
The official fix can be seen here: https://github.com/OpenFOAM/OpenFOAM...1d371bd9be3000
Mojtaba.a and lithos like this.
wyldckat is offline   Reply With Quote

Old   November 21, 2014, 08:54
Thumbs up
  #12
New Member
 
Michael
Join Date: Dec 2011
Location: Geneva
Posts: 28
Rep Power: 5
lithos is on a distinguished road
Quote:
Originally Posted by wyldckat View Post
Dear Bruno, you made my day. Thank you very much for this little detail, explained very understandable. Cheers, Michael

P.S.: I used this with 2.2.x.
wyldckat likes this.

Last edited by lithos; November 21, 2014 at 10:11.
lithos is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 20:56.