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

0002153  OpenFOAM  Bug  public  20160716 21:37  20160723 11:59 
Reporter  bijan darbari  Assigned To  henry  
Priority  normal  Severity  minor  Reproducibility  always 
Status  resolved  Resolution  fixed  
Platform  GNU/Linux  OS  Ubuntu  OS Version  15.10 
Fixed in Version  dev  
Summary  0002153: Activating Brownian force in all lagrangian solvers cause particles to highly speedup  
Description  The simulation of Brownian force between Lagrangian particles for laminar flow suffers from two problems in Openfoam. These problems are introduced as follows: 1 Openfoam applies an inappropriate algorithm in the code to calculate the unit variance zero mean Gaussian random number. 2 In the code, it considers the value of 5.67e8 ev/K as a Boltzmann constant, while the value of 1.38e23 m2Kg/s2K is suggested in reference (e.g. Li, A., & Ahmadi, G. (1992)  Dispersion and deposition of spherical particles from point sources in a turbulent channel flow). The source file can be found at the following adress: src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C) I considered the necessary changes for this force and corrected it accordingly. The corrected Brownian force submodel is attached here. It should be noted that the corrected submodel has been tested and it works to be good very well. Moreover, the corrected submodel is verified with theory equation of below and a good agreement between two results is achieved (See additional information). Theory equation: RMS displacement = (2*D*t)^0.5  
Steps To Reproduce  An absurd value (extremely high) for particle velocity (10^2 m/s to 10^3 m/s ) is observed by activating this force when using any lagrangian solver (e.g. “reactingparcelfoam”, “reactingparcelfilmfoam”, etc.. ). Note that this phenomenon occurs for particle diameters in the range of micrometer or smaller than micrometer. For observation, I attached a simple test case here. These absurd values for velocities can be observed by running the “reactingparcelfoam” solver and activating the Brownian force. Note: I applied the last commit of Brownian force before this bug report.  
Additional Information  Brownian force as a random force between the particles causes the random motion of particles in laminar flow with Lagrangian particles. This force should be considered in calculations for the particles with diameters less than 100 μm as the particles with this size are mainly affected by this force. The simulation of the particle Brownian force in computer programing is introduced systematically by Prof. Goodarz Ahmadi in the following file (Please see page 4 and 5): web2.clarkson.edu/projects/crcd/me537/downloads/2_Brownian.pdf For more information, the authors refer to the paper of ( Li, A., & Ahmadi, G. (1992)  Dispersion and deposition of spherical particles from point sources in a turbulent channel flow). As mentioned in the instruction of this file, the BoxMoller Transform is used to generate the zero mean unit variance Gaussian random number. These numbers are presented by the following expressions: G1=((2*ln(eta1))^0.5)*cos(2*pi*eta2) G2=((2*ln(eta1))^0.5)*sin (2*pi*eta2) where, G1 and G2 are two pair of zero mean unit variance Gaussian random number. Moreover, eta1 and eta2 indicate two uniform random numbers in the range of 01. I utilized this transform for Brownian force submodel of the Openfoam and also replaced the Boltzmann constant by value of 1.38e23. As mentioned in the above file, there is a simple procedure to validate the particles Brownian force. This procedure is briefly introduced as follows: 1 Injection of some particles from a point with a simple flow direction and uniform velocity and temperature 2 Calculation of RMS displacement for a specific time from their injection (Note that this calculation should be performed for displacement of particles normal to direction of the flow) 3 It should be noted that the RMS displacement should be equal to the following expression: RSM displacement = (2*D*t)^0.5 where, t is the time elapsed from particles injection. 4 D = (1.38e23*absolute temperature of fluid)/(3*pi*fluid dynamic viscousity*particle diameter) Moreover, in the validation procedure, the following values have been considered: absolute temperature of fluid = 273 .K fluid dynamic viscousity = 0.001 Pa.s particle diameter = 10, 100, and 1000 nm I performed above procedures for validation in Openfoam and used corrected Brownian force submodel. The variations of RMS displacement with time for three diameters are plotted in the following figure: http://pasteboard.co/chPR5tlaX.png http://pasteboard.co/1ebuMLNZg.png http://pasteboard.co/chRz032EE.png To validate the obtained results by Openfoam, the calculated results by theory formula are also presented in this figure. It can be seen that there is a very good agreement between the results of Openfoam simulation and theory formula. As a result, the new Brownian force submodel in Openfoam works well by considering above corrections.  
Tags  No tags attached.  













Your changes are for an older release of OpenFOAM and it is hard to integrate your changes into OpenFOAMdev. Could you upgrade to OpenFOAMdev? It is not clear why you have changed the method to calculate the random vector, what was wrong with the method Bruno implemented? Your implementation value.Su()[1] = mass*f*G1; value.Su()[2] = mass*f*G2; value.Su()[3] = mass*f*G3; value.Su()[4] = mass*f*G4; is clearly incorrect because it set the elements 14 of a vector with elements 02, could you check? 

I have corrected the Boltzmann constant: OpenFOAMdev commit ac67f6a84bb07ed27897e4cefa724d198965e3e2 and I am running the testcase you supplied but it is not clear how the results should be processed. 

Dear Prof Weller Thank you very much for your kind reply.  Yes, I can upgrade to openfoamdev but I need a bit more time to compile openfoamdev on my system. Please kindly realize that this process takes at least about 24 hr.  The random vector that have been implemented by Bruno didn’t generate zero mean unit variance Gaussian random number efficiently. I’m spent much time to test Brownian force of openfoam and I’m sure that it’s incorrect. For example you can see that scalar “f” twice multiply by a random number “eta” ( in line 194 and 202 ) that is incorrect, because in the reference paper “f” multiply by a random number once. If you test it, you can see that particles highly speed up that is physically unreal. If you are in doubt, I can plot the random numbers generated by Bruno implemented method and Boxmuller method to compare them and see their difference. I have just used the method suggested by Prof Goodarz Ahmadi in the attached PDF file.  First, I used the random vector of 02 but the results were not validated. Finally, I choose random vector of 14 by using trial and error method. It was compiled without any error and the results were validated. Moreover, it is better to use an even random vector, because the BoxMuller transform generates a pair of random number and it is more efficient to use the pair. Process of result: As noted earlier, the RMS displacement of particles in Y direction should be equal to the following theory equation: Y_RMS=(2*D*time)^0.5 D=(Kb*temperature)/(3*pi*dynamic viscosity*diameter) In attached test case, 2000 particles are injected from position of (0.5 0.5) at initial time (t=0). The drag and Brownian forces act on these particles. Accordingly, the particles move only due to Brownian random force and their Y RMS displacement can be used to validate the process. Temperature = 273 K Kb = 1.38e23 Pi = 3.1415 dynamic viscousity = 0.001 Pa.s diameter = 1e6 m then: D = 4e13 When you run the test case, you should compute particles YRMS displacement in a specific time from the start of solution by using the following correlation: YRMS = ((1/number of particles )*(sum (Y(i)Y0)^2))^0.5 for i=1,number of particles In the test case n=2000 and Y0=0.5 . I used a simple fortran code. This code is attached here. This code reads the particles positions file and computes YRMS. For example at t=1 sec: 1 Run the test case 2 Copy the file “testcase/1/lagrangian/reactingcloud1/positions” and paste it in the fortran project directory. 3 Run the fortran code. This code reads the positions file and shows value “ Y RMS = 8.88e7” 4 Use the theory equation and calculate Y RMS = (2*4e13*1)^0.5 = 8.94e7 5 You can see that the RMS displacement is validated for t=1 sec. 6 You can repeat this step for other times. 

YRMSfortran.f90 (507 bytes)
program RMSposition implicit none real,dimension(2000)::y character(len=20)::pr1 real::sumy integer::i open(UNIT = 1, FILE = "positions", ACTION = "read", position = "ASIS") do i=1,19 read(1,*) end do do i=1,2000 read(1,*),pr1,y(i) end do sumy=0 do i=1,2000 sumy=sumy+((y(i)0.5)**2) end do sumy=sqrt(sumy/2000) write(*,*),' RMS Displacement = ',sumy read(*,*) end 

> For example you can see that scalar “f” twice multiply by a random number “eta” ( in line 194 and 202 ) that is incorrect Agreed, and if you take the first eta out are the results then correct? > Finally, I choose random vector of 14 by using trial and error method. It was compiled without any error and the results were validated. This is clearly incorrect; if you compile with fulldebug and range checking you will find the code stops with an error accessing the vector elements out of range and the first element is not set. 

I have removed the spurious additional 'eta': commit 96645225ed6e088d5ab77fc3dd5411cd28994cff 

Dear Prof Weller Thank you for your fast reply > Agreed, and if you take the first eta out are the results then correct? No. I tried that but the result were incorrect. > This is clearly incorrect; if you compile with fulldebug and range checking you will find the code stops with an error accessing the vector elements out of range and the first element is not set. Sorry. I'm not expert in C++. only I know that brownian random forces must act on particles sequential ( each pair for one direction ) and my result validated for Y direction. 

The force vector has 3 components, you cannot set 4 irrespective of the programming language you use. 

Sorry. You are right. but I say the truth. the result were validated for Y direction. you can test that. may be just set component 1 and 2 of the force vector cause the result to validated for Y direction. 

Why are you calculating the RMS of the ydisplacement? Why not calculate the RMS of the vector displacement? 

Because in the file Brownian.PDF ,YRMS were calculated for validation process. 

lagrangian::BrownianMotionForce: Changed from a cubic to a spherical distribution See also: StochasticDispersionRAS Resolved in OpenFOAMdev by commit 46d69e1deaed6b035b96524a12670ee0602da542 Note that to compare with the RMS displacement in the 1D analysis cited you need to calculate the RMS of the magnitude of the displacement in the spherical distribution. 
Date Modified  Username  Field  Change 

20160716 21:37  bijan darbari  New Issue  
20160716 21:37  bijan darbari  File Added: Testcase.zip  
20160716 21:37  bijan darbari  File Added: BrownianMotion.zip  
20160716 21:38  bijan darbari  File Added: 10nm.png  
20160716 21:38  bijan darbari  File Added: 100nm.png  
20160716 21:38  bijan darbari  File Added: 1000nm.png  
20160716 21:40  bijan darbari  File Added: Brownian.pdf  
20160722 15:43  henry  Note Added: 0006558  
20160722 16:19  henry  Note Added: 0006559  
20160722 20:43  bijan darbari  Note Added: 0006561  
20160722 20:45  bijan darbari  File Added: YRMSfortran.f90  
20160722 20:52  henry  Note Added: 0006562  
20160722 21:13  henry  Note Added: 0006563  
20160722 21:31  bijan darbari  Note Added: 0006564  
20160722 21:35  henry  Note Added: 0006565  
20160722 21:50  bijan darbari  Note Added: 0006566  
20160722 21:52  henry  Note Added: 0006567  
20160722 21:57  bijan darbari  Note Added: 0006568  
20160723 11:58  henry  Note Added: 0006571  
20160723 11:58  henry  Status  new => resolved 
20160723 11:58  henry  Fixed in Version  => dev 
20160723 11:58  henry  Resolution  open => fixed 
20160723 11:58  henry  Assigned To  => henry 