Skip to content

Brownian motion force should use cubic distribution

Summary

The current code in BrownianMotionForce.C is calculating wrong Brownian motion force. The equation used to be correct (now it's commented out). However this reporter (ee9e7cf0,https://bugs.openfoam.org/view.php?id=2153) didn't know how to properly set up a particle tracking simulation and used an unrealistically large time step to generate wrong results. Unfortunately Mr. Henry Weller believed in him and adapted his wrong code. The ESI version is also affected.

Steps to reproduce

I'm attaching an example case to reproduce. The domain is 2um x 2um x 2um box with stationary water inside. 2000 x 1um particles are injected at (0,0,0) with BrownianMotion and sphereDrag enabled. After some time, collect all coordinates and calculate the RMS of the particles (I use paraview+excel).

Example case

This case takes 2 hours to run on my workstation. It might be OK to reduce the total time. However the time step should always be smaller than particle relaxation time (5e-8s in this case).

Tesetcase-v2212.zip

blockMesh
decomposePar
mpirun -np 8 reactingParcelFoam -parallel

What is the current bug behaviour?

The spherical distribution is generating a Brownian motion force that is sqrt(3) times smaller than it should be. The validation case provided by the original reporter made 2 mistakes: 1. Use 2D geometry to validate 3D Brownian motion. 2. Use super large timestep (20,000 times of particle relaxation time) to simulate Brownian motion.

What is the expected correct behavior?

Literature 1: Dispersion and Deposition of Spherical Particles from Point Sources in a Turbulent Channel Flow https://www.tandfonline.com/doi/ref/10.1080/02786829208959550

Literature 2: https://dl.cfdexperts.net/cfd_resources/Ansys_Documentation/Fluent/Ansys_Fluent_Theory_Guide.pdf (page 425) Fluent is using the "cubic distribution" as well.

The RMS of the particle should match the theoretical value from diffusion coefficient (https://en.wikipedia.org/wiki/Random_walk#Relation_to_Wiener_process), where RMS^2 should equal to 6Dt. Here T=293, mu=1e-3, d=1e-6, therefore using Stokes-Einstein equation D=4.3e-13.

Relevant logs and/or images

RMS

Environment information

  • OpenFOAM version : all versions since 2016
  • Operating system :
  • Hardware info :
  • Compiler :

Possible fixes

Use the following code in src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C will generate correct results.

    Random& rnd = this->owner().rndGen();
    for (direction dir = 0; dir < vector::nComponents; dir++)
    {
        value.Su()[dir] = f*rnd.GaussNormal<scalar>();
    }

This should be faster than the implementation that is commented out, as rnd.GaussNormal<scalar>() generates a pair of random numbers.

Edited by Guanyang Xue