Sphere settling at Re = 1.5

Hello,

I'm trying to reconstruct cases described in article: "Simulation of Sedimentation of a Sphere in a Viscous Fluid Using the Lattice Boltzmann Method Combined with the Smoothed Profile Method" by Suresh Alapati, Woo Seong Che, and Yong Kweon Suh, and here A. ten Cate, C. H. Nieuwstad, J. J. Derksen, and H. E. A. van den Akker, “Particle imaging velocimetry experiments and lattice-Botlzmann simulations on a single sphere settling under gravity,” Physics of Fluids, vol. 14, no. 11, pp. 4012–4025, 2002.

Paper introduces the LBM method capable of handling solid particles within the fluid. The validation case is a particle settling due to gravitational force.

I would like to use cfdemSolverIB for that.
I have problem because the velocity of the particle is not constant after some time and I'm searching for the root cause.

The domain size is (0.1, 0.1, 0.16). On all boundaries I've set up the no-slip BC (U = (0, 0, 0) and p = zeroGradient). Domain is divided into (50, 50, 75) cells. The picture of voidFraction filed is attached to this message. Base on that I think that the resolution is acceptable for the first tests.

The Reynolds number is achieved by adjusting kinematic viscosity. Using simple formula for Re with particle diameter = 0.015, u_{\infty} = 0.038 m/s and Re = 1.5 we get \nu = 3.8e-5. As it is written in the papers the u_{\infty} is calculated using Stokes theorem and it is particle terminal velocity in case of infinite channel. In my simulations at some point I expect that particle velocity will be somewhere around this value, but unfortunately it is much larger and I can't see why.

The time step for CFD is dtc = 1.32E-5, and for DEM side dtd = 1.32e-6; It is calculated taking 10% of Rayleigh time for nylon particle.

I'm not using ArhimedesIB forcing model therefore I assumed that resultant acceleration acting on a particle will be reduced m_p * a = F_b- F_g where m_p is particle mass (Volume * \rho_p) F_b = V \rho_f g - buoyancy force and F_g = V * rho_p * g is the gravity force where g is gravitational acceleration. This gives me a = g * ( (\rho_f / \rho_s) - 1) which for \rho_f = 970 and \rho_s = 1120 a = 1.517 m / s^2

Additional questions:

• Is the cfdemSolverIB capable for Re ~= 1 to 11? In this document http://www.cfd.com.au/cfd_conf15/PDFs/140GON.pdf It is written that this solver fails for low Re but is mean around 1 or much less than 1 like stokes flow?
• Will the cfdemSolverForceIB be released in PUBLIC repository?
• About ArhimedesIB. If I enable ArhimedesIB forces I should use normal gravitational acceleration (9.81) in my simulation setup both LIGGHTS config file and in constant / g file is it right?

I would be verry happy if someone could give me some hints what is wrong with my simulation.
The simulation case is also attached to this message.

Thank you in advance for your help.
Regards,
Jakub.