View Issue Details

IDProjectCategoryView StatusLast Update
0002036OpenFOAMBugpublic2016-04-25 08:07
Reporterbijan darbari Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityhave not tried
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSUbuntuOS Version15.04
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)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Class
    Foam::BrownianMotionForce

Description
    Calculates particle Brownian motion force.

    From:
    \verbatim
        \article{ALi_GAhmadi,
        author = {Amy Li and Goodarz Ahmadi},
        title = {Dispersion and Deposition of Spherical Particles from Point
                 Sources in a Turbulent Channel Flow},
        journal = {Aerosol Science and Technology},
        volume = {16},
        number = {4},
        pages = {209-226},
        year = {1992},
        url = {http://dx.doi.org/10.1080/02786829208959550},
        doi = {10.1080/02786829208959550}
        }
    \endverbatim


SourceFiles
    BrownianMotionForceI.H
    BrownianMotionForce.C

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

#ifndef BrownianMotionForce_H
#define BrownianMotionForce_H

#include "ParticleForce.H"
#include "cachedRandom.H"

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

namespace Foam
{

/*---------------------------------------------------------------------------*\
                    Class BrownianMotionForce Declaration
\*---------------------------------------------------------------------------*/

template<class CloudType>
class BrownianMotionForce
:
    public ParticleForce<CloudType>
{
    // Private data

        //- Reference to the cloud random number generator
        cachedRandom& rndGen_;

        //- Molecular free path length [m]
        const scalar lambda_;

        //- Turbulence flag
        bool turbulence_;

        //- Pointer to the turbulence kinetic energy field
        const volScalarField* kPtr_;

        //- Flag that indicates ownership of turbulence k field
        bool ownK_;


    // Private Member Functions

        //- Inverse erf for Gaussian distribution
        scalar erfInv(const scalar y) const;

        //- Return the k field from the turbulence model
        tmp<volScalarField> kModel() const;


public:

    //- Runtime type information
    TypeName("BrownianMotion");


    // Constructors

        //- Construct from mesh
        BrownianMotionForce
        (
            CloudType& owner,
            const fvMesh& mesh,
            const dictionary& dict
        );

        //- Construct copy
        BrownianMotionForce(const BrownianMotionForce& bmf);

        //- Construct and return a clone
        virtual autoPtr<ParticleForce<CloudType>> clone() const
        {
            return autoPtr<ParticleForce<CloudType>>
            (
                new BrownianMotionForce<CloudType>(*this)
            );
        }


    //- Destructor
    virtual ~BrownianMotionForce();


    // Member Functions

        // Access

            //- Return const access to the molecular free path length [m]
            inline scalar lambda() const;

            //- Return const access to the turbulence flag
            inline bool turbulence() const;


        // Evaluation

            //- Cache fields
            virtual void cacheFields(const bool store);

            //- 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;
};


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

} // End namespace Foam

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

#include "BrownianMotionForceI.H"

#ifdef NoRepository
    #include "BrownianMotionForce.C"
#endif

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

#endif

// ************************************************************************* //
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

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

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

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
2016-04-24 19:47 bijan darbari Note Edited: 0006166
2016-04-24 19:47 henry Note Added: 0006167
2016-04-24 19:47 bijan darbari Note Edited: 0006166
2016-04-24 19:47 henry Note Edited: 0006167
2016-04-25 07:49 bijan darbari Note Added: 0006170
2016-04-25 07:50 bijan darbari Note Edited: 0006170
2016-04-25 07:57 bijan darbari Note Edited: 0006170
2016-04-25 08:00 bijan darbari Note Edited: 0006170
2016-04-25 08:02 bijan darbari Note Edited: 0006170
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