View Issue Details

IDProjectCategoryView StatusLast Update
0002153OpenFOAMBugpublic2016-07-23 11:59
Reporterbijan darbari Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSUbuntuOS Version15.10
Fixed in Versiondev 
Summary0002153: Activating Brownian force in all lagrangian solvers cause particles to highly speedup
DescriptionThe 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:
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 ReproduceAn 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 InformationBrownian 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 Box-Moller 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 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:

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.
TagsNo tags attached.

Activities

bijan darbari

2016-07-16 21:37

reporter  

Testcase.zip (284,131 bytes)

bijan darbari

2016-07-16 21:37

reporter  

BrownianMotion.zip (4,619 bytes)

bijan darbari

2016-07-16 21:38

reporter  

10nm.png (30,720 bytes)   
10nm.png (30,720 bytes)   

bijan darbari

2016-07-16 21:38

reporter  

100nm.png (37,963 bytes)   
100nm.png (37,963 bytes)   

bijan darbari

2016-07-16 21:38

reporter  

1000nm.png (31,417 bytes)   
1000nm.png (31,417 bytes)   

bijan darbari

2016-07-16 21:40

reporter  

Brownian.pdf (230,134 bytes)

henry

2016-07-22 15:43

manager   ~0006558

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?

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 1-4 of a vector with elements 0-2, could you check?

henry

2016-07-22 16:19

manager   ~0006559

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.

bijan darbari

2016-07-22 20:43

reporter   ~0006561

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:

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

bijan darbari

2016-07-22 20:45

reporter  

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)   

henry

2016-07-22 20:52

manager   ~0006562

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

henry

2016-07-22 21:13

manager   ~0006563

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

bijan darbari

2016-07-22 21:31

reporter   ~0006564

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.

henry

2016-07-22 21:35

manager   ~0006565

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

bijan darbari

2016-07-22 21:50

reporter   ~0006566

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.

henry

2016-07-22 21:52

manager   ~0006567

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

bijan darbari

2016-07-22 21:57

reporter   ~0006568

Because in the file Brownian.PDF ,Y-RMS were calculated for validation process.

henry

2016-07-23 11:58

manager   ~0006571

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.

Issue History

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