# SaffmanMeiLiftForce

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

 October 13, 2013, 11:16 SaffmanMeiLiftForce #1 Member   Evangelos Join Date: Sep 2011 Posts: 67 Rep Power: 5 Hi Can anyone provide an example of SaffmanMeiLiftForce dictionary ? thanks a lot !

 October 13, 2013, 14:51 #2 Super Moderator   Bruno Santos Join Date: Mar 2009 Location: Lisbon, Portugal Posts: 8,253 Blog Entries: 34 Rep Power: 84 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: The "banana" way: http://openfoamwiki.net/index.php/Op...de/Use_bananas The source code way, namely by looking at the source code and or the code documentation: http://www.openfoam.org/docs/cpp/ Best regards, Bruno __________________ OpenFOAM: Frequently Asked Questions | Useful links for building and using Forum: How to ask for help | Posting code and output with [CODE] When will I answer questions? Check this page for dates: http://wyldckat.github.io And please: Read this before sending private messages to me

 October 13, 2013, 16:00 #3 Member   Evangelos Join Date: Sep 2011 Posts: 67 Rep Power: 5 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

 October 14, 2013, 04:06 #4 Member   Evangelos Join Date: Sep 2011 Posts: 67 Rep Power: 5 I need something like that pressureGradient dictionary but for saffmanmeiliftforce!!!

 October 14, 2013, 16:47 #5 Super Moderator   Bruno Santos Join Date: Mar 2009 Location: Lisbon, Portugal Posts: 8,253 Blog Entries: 34 Rep Power: 84 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 __________________ OpenFOAM: Frequently Asked Questions | Useful links for building and using Forum: How to ask for help | Posting code and output with [CODE] When will I answer questions? Check this page for dates: http://wyldckat.github.io And please: Read this before sending private messages to me

 October 4, 2014, 23:34 #6 Senior Member     Join Date: Jan 2010 Posts: 347 Blog Entries: 2 Rep Power: 8 Hi all, Do you find any dictionary for saffmanMeiLiftForce?

October 5, 2014, 04:20
#7
Super Moderator

Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,253
Blog Entries: 34
Rep Power: 84
Greetings Maysam,

Quote:
 Originally Posted by maysmech 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 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

 October 5, 2014, 04:33 #8 Senior Member     Join Date: Jan 2010 Posts: 347 Blog Entries: 2 Rep Power: 8 For example for DPMFoam

October 5, 2014, 06:58
#9
Super Moderator

Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,253
Blog Entries: 34
Rep Power: 84
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
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```
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 Foam::scalar Foam::SaffmanMeiLiftForce::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

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

 October 5, 2014, 20:25 #10 Senior Member     Join Date: Jan 2010 Posts: 347 Blog Entries: 2 Rep Power: 8 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.

 October 18, 2014, 13:38 #11 Super Moderator   Bruno Santos Join Date: Mar 2009 Location: Lisbon, Portugal Posts: 8,253 Blog Entries: 34 Rep Power: 84 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.

November 21, 2014, 08:54
#12
New Member

Michael
Join Date: Dec 2011
Location: Geneva
Posts: 28
Rep Power: 5
Quote:
 Originally Posted by wyldckat 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
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.

Last edited by lithos; November 21, 2014 at 10:11.

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

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