View Issue Details
ID  Project  Category  View Status  Date Submitted  Last Update 

0002036  OpenFOAM  Bug  public  20160329 20:37  20160425 08:07 
Reporter  bijan darbari  Assigned To  henry  
Priority  normal  Severity  minor  Reproducibility  have not tried 
Status  resolved  Resolution  fixed  
Platform  GNU/Linux  OS  Ubuntu  OS Version  15.04 
Fixed in Version  dev  
Summary  0002036: 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 194197): https://github.com/OpenFOAM/OpenFOAM3.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??  
Tags  No tags attached.  

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

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 OpenFOAMhistory repository, where we can find that this class was first introduced in this commit: https://github.com/OpenCFD/OpenFOAMhistory/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:209226, 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. 

BrownianMotionForce.H (4,850 bytes)
/**\ =========  \\ / F ield  OpenFOAM: The Open Source CFD Toolbox \\ / O peration  \\ / A nd  Copyright (C) 20112016 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 = {209226}, 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 // ************************************************************************* // 

@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. 

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). Oneequation subgrid scale (SGS) modelling for Euler–Euler large eddy simulation (EELES) of dispersed bubbly flow. Chemical Engineering Science, 63(15), 39233931. \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. 

@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. 

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. 

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? 

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 openfoamdev 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. 
Date Modified  Username  Field  Change 

20160329 20:37  bijan darbari  New Issue  
20160331 18:13  bijan darbari  Note Added: 0006067  
20160417 00:41  wyldckat  Note Added: 0006134  
20160417 00:55  wyldckat  File Added: BrownianMotionForce.H  
20160417 00:58  wyldckat  Note Added: 0006135  
20160417 15:20  henry  Note Added: 0006136  
20160417 15:42  henry  Note Added: 0006137  
20160424 19:44  bijan darbari  Note Added: 0006166  
20160424 19:45  bijan darbari  Note Edited: 0006166  
20160424 19:47  bijan darbari  Note Edited: 0006166  
20160424 19:47  henry  Note Added: 0006167  
20160424 19:47  bijan darbari  Note Edited: 0006166  
20160424 19:47  henry  Note Edited: 0006167  
20160425 07:49  bijan darbari  Note Added: 0006170  
20160425 07:50  bijan darbari  Note Edited: 0006170  
20160425 07:57  bijan darbari  Note Edited: 0006170  
20160425 08:00  bijan darbari  Note Edited: 0006170  
20160425 08:02  bijan darbari  Note Edited: 0006170  
20160425 08:06  henry  Status  new => resolved 
20160425 08:06  henry  Fixed in Version  => dev 
20160425 08:06  henry  Resolution  open => fixed 
20160425 08:06  henry  Assigned To  => henry 