CFD Online Discussion Forums - About the dimension of Saffman-Mei particle lift force model

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - About the dimension of Saffman-Mei particle lift force model (https://www.cfd-online.com/Forums/openfoam-programming-development/116510-about-dimension-saffman-mei-particle-lift-force-model.html)

About the dimension of Saffman-Mei particle lift force model

I'm working with flows with discrete, Lagrangian, particles in OpenFOAM. In the latest 2.2-release of OF, there is a particle lift force model, i.e. cross-stream forces on particles from velocity gradients.

However, I am slightly unsure if the implemented force model can be correct. I have not yet filed a bug report on it, in case my interpretation of the code is wrong. The piece of code I am looking at is found in src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C, line 144:

Code:

value.Su() = mass/p.rho()*p.d()/2.0*p.rhoc()*Cl*((p.Uc() - p.U())^curlUc);

If I use "normal" units for the terms mass, rho, d, Uc and curlUc, and assume that the lift coefficient is dimensionless, i get that the lift force gets the unit kg*m^2/s^2. How can this be correct? Shouldn't it be kg*m/s^2, i.e. mass times acceleration?

The "correct" behaviour can for example be seen in src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Gravity/GravityForce.C:

Code:

value.Su() = mass*g_*(1.0 - p.rhoc()/p.rho());

As g_ is an acceleration (with unit m/s^2), the dimension of this expression is clearly mass times acceleration.

Does any of you have any reasonable explanation on how the expression in LiftForce.C can be correct?

I have now come to the conclusion that it defineatly is a bug in the above mentioned piece of code. The bug is reported here: http://www.openfoam.org/mantisbt/view.php?id=822

For the reference, removing the particle diameter from the expression gives you the correct force. See the bug report above for details. The corrected line of code then become:

Code:

value.Su() = mass/p.rho()*p.rhoc()/2*Cl*((p.Uc() - p.U())^curlUc);