View Issue Details

IDProjectCategoryView StatusLast Update
0002036OpenFOAM[All Projects] Bugpublic2016-04-25 08:07
Reporterbijan darbariAssigned Tohenry 
PrioritynormalSeverityminorReproducibilityhave not tried
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSUbuntuOS Version15.04
Product Version 
Fixed in Versiondev 
Summary0002036: BrownianMotionForce of particles in lamiar flow, incorrect correlation
Description
In DPM model. BrownianMotionForce is one of the main forces that act on particles.

look at this file (line 194-197):
https://github.com/OpenFOAM/OpenFOAM-3.0.x/blob/master/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C

the BrownianMotionForce for laminar flow have been defined as:
""const scalar rhoRatio = p.rho()/p.rhoc();
  const scalar s0 =
            216*muc*sigma*Tc/(sqr(mathematical::pi)*pow5(dp)*(rhoRatio)*cc);
  f = eta*sqrt(mathematical::pi*s0/dt);""


in lamiar flow, the correlation for this force is:
http://uploadgoogle.ir/uploads/145915541121.jpg

note: in above picture, rho and rho_f is for continuous phase and rho_p is for particles.

If you look to above picture, in the denominator of S0, the square of rho_p (particle density) exist. but when I look into BrownianMotionForce.C, the power one of rho_p is in the denomiantor of S0!!! ( hope I don't mistake )

Is this a bug in code??

TagsNo tags attached.

Activities

bijan darbari

2016-03-31 18:13

reporter   ~0006067

Excuse me for poor English. If possible please reply this post.

wyldckat

2016-04-17 00:41

updater   ~0006134

OK, this is a really tough question. Specially because the code documentation in OpenFOAM in this class is missing where this data actually came from...
Although we at least have the OpenFOAM-history repository, where we can find that this class was first introduced in this commit: https://github.com/OpenCFD/OpenFOAM-history/commit/36441393ee98317fb9dc0aee04c9f592d460a7bc


So I've done some online searching for the actual expression and I ended up finding in Fluent's documentation a reference to the following paper for the same expression:

  A. Li and G. Ahmadi.
  Dispersion and Deposition of Spherical Particles from Point Sources in a Turbulent Channel Flow.
  Aerosol Science and Technology, 16:209-226, 1992.

Luckily it can be found available online for free here: http://www.tandfonline.com/doi/abs/10.1080/02786829208959550 - and it's not the only place where this is available.


Now I see what you mean, because according to that paper, the expression in the source code should actually be this:

  const scalar rhoRatio = p.rho()/p.rhoc();
  const scalar s0 =
            216*muc*sigma*Tc/(sqr(mathematical::pi)*p.rhoc()*pow5(dp)*sqr(rhoRatio)*cc);

Problem is that it's very possible that some simplifications were put into practice in the code and figuring them out can be a bit of a pain... although... hold on, let's see, the paper does state that it should be "nu" and not "mu", therefore, it should in fact be this:

  const scalar rhoRatio = p.rho()/p.rhoc();
  const scalar s0 =
            216*nuc*sigma*Tc/(sqr(mathematical::pi)*p.rhoc()*pow5(dp)*sqr(rhoRatio)*cc);

Notice that instead of "muc" it has "nuc", as in the paper. But since:

  nuc = muc/p.rhoc()

means that it then becomes:

  const scalar rhoRatio = p.rho()/p.rhoc();
  const scalar s0 =
            216*muc*sigma*Tc/(sqr(mathematical::pi)*sqr(p.rhoc())*pow5(dp)*sqr(rhoRatio)*cc);

which by dropping "rhoRatio", it would become this:

  const scalar s0 =
            216*muc*sigma*Tc/(sqr(mathematical::pi)*sqr(p.rhoc())*pow5(dp)*sqr(p.rho()/p.rhoc())*cc);

And finally, by simplifying terms:

  const scalar s0 =
            216*muc*sigma*Tc/(sqr(mathematical::pi)*pow5(dp)*sqr(p.rho())*cc);



I need a second and maybe even a third opinion on this one. It's very possible that when the implementation was done in OpenFOAM, a mistake was accidentally done when simplifying the equation.

wyldckat

2016-04-17 00:55

updater  

BrownianMotionForce.H (4,850 bytes)

wyldckat

2016-04-17 00:58

updater   ~0006135

@Henry: In the meantime, until someone can take a second look into this, attached is the file "BrownianMotionForce.H", meant to update the file "src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H":

  - I added the information about the document that I managed to find and used the convention used in the file "azizChen.H", namely to use the Bibtex convention of citing a paper/article.

  - Also fixed a typo in the comment for the "calcCoupled" method.

henry

2016-04-17 15:20

manager   ~0006136

Thanks for looking into this.

"convention used in the file "azizChen.H", namely to use the Bibtex convention of citing a paper/article"

It appears the author of this file chose to use a reference "convention" different to the rest of OpenFOAM. After some consideration of the alternatives we have settled on the APA citation style but not yet completed the update of all references in OpenFOAM. For examples of the APA style take a look at the TurbulenceModel library, e.g.

    Reference:
    \verbatim
        Niceno, B., Dhotre, M. T., & Deen, N. G. (2008).
        One-equation sub-grid scale (SGS) modelling for
        Euler–Euler large eddy simulation (EELES) of dispersed bubbly flow.
        Chemical Engineering Science, 63(15), 3923-3931.
    \endverbatim

I have updated BrownianMotionForce.H and azizChen.H (even the naming of this file is inconsistent with the OpenFOAM conventions) and will push shortly.

henry

2016-04-17 15:42

manager   ~0006137

@Bruno, I have studied the paper and agree with your assessment. There is absolutely no reason for the particle density to appear in the expression at all.

I have applied the change:
commit 04fe88592012870be5cacb94cb5c9ce2aa84339e

@bijan darbari: Could you test this change and if all is well I will close this report.

bijan darbari

2016-04-24 19:44

reporter   ~0006166

Last edited: 2016-04-24 19:47

View 4 revisions

Dear Dr @henry
Thanks your very much for your attention.
I tested and it worked well. so you can close this report. although I don't know how the class generate "gaussian random number with zero mean and unit variance"!!!

Dear @wyldckat
thank you very much for your understanding. I was disappointed about this bug. because nobody reply my post for a long time. and your reply make me so happy.

henry

2016-04-24 19:47

manager   ~0006167

Last edited: 2016-04-24 19:47

View 2 revisions

It is not clear what you mean by

> I tested and it worked well. although there is some ambiguity about that how
> the class generate "gaussian random number with zero mean and unit variance"!!!

could you clarify and let me know what changes are needed?

bijan darbari

2016-04-25 07:49

reporter   ~0006170

Last edited: 2016-04-25 08:02

View 5 revisions

Dear Dr @henry

Now I saying these class is true for these reasons:

In openfoam 3.0.1, particles escaped back from inlet, when I actived the Brownian force in my dpm simulation. After this commit. I replaced BrownianMotionForce.C with edited one in openfoam-dev version and compiled that library.I tested again with the same conditions and no any particles escaped back from the inlet, so I think Brownian force work well in this class and this is the true one. also the correlation for S0 is now true.


About the ambiguity:
As I mentioned earlier, this is the expression for the Brownian force
 
link: http://uploadgoogle.ir/uploads/145915541121.jpg

As you see at above, the square multiply by Gi, which are three Gaussian random numbers with zero mean and unit variance (from reference paper by “Li and Ahmadi”). These are special sequence of random numbers approximately between -1 and 1.
but I don’t know how the code generate this sequence of random numbers. Mybe these numbers generated by below lines of BrownianMotionForce.C that convert series of uniform random numbers between 0 and 1 to series of Gaussian random numbers with zero mean and unit variance.

const scalar sqrt2 = sqrt(2.0);

for (label i = 0; i < 3; i++)
{
const scalar x = rndGen_.sample01<scalar>();
const scalar eta = sqrt2*erfInv(2*x - 1.0);
value.Su()[i] = mass*f*eta;
}


I didn’t mentioned that there is still something wrong, but I guess you know how these series of Gaussian random numbers with zero mean and unit variance generated by this class.
In my opinion this class is now true and there is no need to do other changes on it.

Excuse me. mybe my poor english cause for you to don't understand well what I sayed.

Issue History

Date Modified Username Field Change
2016-03-29 20:37 bijan darbari New Issue
2016-03-31 18:13 bijan darbari Note Added: 0006067
2016-04-17 00:41 wyldckat Note Added: 0006134
2016-04-17 00:55 wyldckat File Added: BrownianMotionForce.H
2016-04-17 00:58 wyldckat Note Added: 0006135
2016-04-17 15:20 henry Note Added: 0006136
2016-04-17 15:42 henry Note Added: 0006137
2016-04-24 19:44 bijan darbari Note Added: 0006166
2016-04-24 19:45 bijan darbari Note Edited: 0006166 View Revisions
2016-04-24 19:47 bijan darbari Note Edited: 0006166 View Revisions
2016-04-24 19:47 henry Note Added: 0006167
2016-04-24 19:47 bijan darbari Note Edited: 0006166 View Revisions
2016-04-24 19:47 henry Note Edited: 0006167 View Revisions
2016-04-25 07:49 bijan darbari Note Added: 0006170
2016-04-25 07:50 bijan darbari Note Edited: 0006170 View Revisions
2016-04-25 07:57 bijan darbari Note Edited: 0006170 View Revisions
2016-04-25 08:00 bijan darbari Note Edited: 0006170 View Revisions
2016-04-25 08:02 bijan darbari Note Edited: 0006170 View Revisions
2016-04-25 08:06 henry Status new => resolved
2016-04-25 08:06 henry Fixed in Version => dev
2016-04-25 08:06 henry Resolution open => fixed
2016-04-25 08:06 henry Assigned To => henry