View Issue Details
|ID||Project||Category||View Status||Date Submitted||Last Update|
|0002153||OpenFOAM||Bug||public||2016-07-16 21:37||2016-07-23 11:59|
|Reporter||bijan darbari||Assigned To||henry|
|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.67e-8 ev/K as a Boltzmann constant, while the value of 1.38e-23 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:
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):
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 Box-Moller Transform is used to generate the zero mean unit variance Gaussian random number. These numbers are presented by the following expressions:
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 0-1. I utilized this transform for Brownian force submodel of the Openfoam and also replaced the Boltzmann constant by value of 1.38e-23.
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.38e-23*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:
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.|
Testcase.zip (284,131 bytes)
BrownianMotion.zip (4,619 bytes)
1000nm.png (31,417 bytes)
1000nm.png (31,417 bytes)
Brownian.pdf (230,134 bytes)
Your changes are for an older release of OpenFOAM and it is hard to integrate your changes into OpenFOAM-dev. Could you upgrade to OpenFOAM-dev?
It is not clear why you have changed the method to calculate the random vector, what was wrong with the method Bruno implemented?
value.Su() = mass*f*G1;
value.Su() = mass*f*G2;
value.Su() = mass*f*G3;
value.Su() = mass*f*G4;
is clearly incorrect because it set the elements 1-4 of a vector with elements 0-2, could you check?
I have corrected the Boltzmann constant:
OpenFOAM-dev commit ac67f6a84bb07ed27897e4cefa724d198965e3e2
and I am running the test-case 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 openfoam-dev but I need a bit more time to compile openfoam-dev 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 Box-muller 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 0-2 but the results were not validated. Finally, I choose random vector of 1-4 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 Box-Muller 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:
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.38e-23
Pi = 3.1415
dynamic viscousity = 0.001 Pa.s
diameter = 1e-6 m
then: D = 4e-13
When you run the test case, you should compute particles Y-RMS displacement in a specific time from the start of solution by using the following correlation:
Y-RMS = ((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 Y-RMS.
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.88e-7”
4- Use the theory equation and calculate Y RMS = (2*4e-13*1)^0.5 = 8.94e-7
5- You can see that the RMS displacement is validated for t=1 sec.
6- You can repeat this step for other times.
Y-RMS-fortran.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
Y-RMS-fortran.f90 (507 bytes)
> 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 1-4 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 full-debug 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 full-debug 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 y-displacement? Why not calculate the RMS of the vector displacement?|
|Because in the file Brownian.PDF ,Y-RMS were calculated for validation process.|
lagrangian::BrownianMotionForce: Changed from a cubic to a spherical distribution
See also: StochasticDispersionRAS
Resolved in OpenFOAM-dev 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.
|2016-07-16 21:37||bijan darbari||New Issue|
|2016-07-16 21:37||bijan darbari||File Added: Testcase.zip|
|2016-07-16 21:37||bijan darbari||File Added: BrownianMotion.zip|
|2016-07-16 21:38||bijan darbari||File Added: 10nm.png|
|2016-07-16 21:38||bijan darbari||File Added: 100nm.png|
|2016-07-16 21:38||bijan darbari||File Added: 1000nm.png|
|2016-07-16 21:40||bijan darbari||File Added: Brownian.pdf|
|2016-07-22 15:43||henry||Note Added: 0006558|
|2016-07-22 16:19||henry||Note Added: 0006559|
|2016-07-22 20:43||bijan darbari||Note Added: 0006561|
|2016-07-22 20:45||bijan darbari||File Added: Y-RMS-fortran.f90|
|2016-07-22 20:52||henry||Note Added: 0006562|
|2016-07-22 21:13||henry||Note Added: 0006563|
|2016-07-22 21:31||bijan darbari||Note Added: 0006564|
|2016-07-22 21:35||henry||Note Added: 0006565|
|2016-07-22 21:50||bijan darbari||Note Added: 0006566|
|2016-07-22 21:52||henry||Note Added: 0006567|
|2016-07-22 21:57||bijan darbari||Note Added: 0006568|
|2016-07-23 11:58||henry||Note Added: 0006571|
|2016-07-23 11:58||henry||Status||new => resolved|
|2016-07-23 11:58||henry||Fixed in Version||=> dev|
|2016-07-23 11:58||henry||Resolution||open => fixed|
|2016-07-23 11:58||henry||Assigned To||=> henry|