Units for Shan-Chen interaction strength parameter

Lattice Boltzmann Method LBM: Theory
#1

Hi All,

I am in the process of converting an LBM code from simulation to physical units. I know the basics of the procedure, but I run into some difficulties when working with the Shan-Chen models for multiphase/multicomponent flows. In order to transform the pressure to physical units I use the non-ideal EOS as shown in Sukop and Thorne.
This reads: p = c_2^2rho + c_s^2G/2psi(rho)^2. However, when I set G to achieve dimensional consistency, its units differ from those needed for dimensional consistency of the interaction force density. For simplicity, following Martys and Chen 1996, I set psi(rho) = rho. This gives units of G as L^5M^{-1}T^{-2}. Now, if I use the same effective mass (psi(rho)=rho) in the interaction strength equation, which reads F = -Gpsisum_a(w_ipsie_i), G has units L^3M^{-1}*T^{-1} ! In the force equation w_i is dimensionless and e_i has units of speed.

Has anyone had any similar problems? Any advice will be greatly appreciated.

Thank you,

Peter

#2

Dear PeterH,

this is a really good question. The answer is not so trivial and longer than you think. I am assuming that T=time, M=mass, L=length. With this in mind, I end up with

[G] = L^3 * M^-1

which is the inverse of the units of the density. And this actually fits both equations: the interaction strength equation and the pressure equation.
First important thing is that the shan chen force is actually a volume force, thus, its units are Newton/meter^3 or in your notation M^1 * L^-2 * T^-2. Now, the force equation reads

F = -G \psi sum_i w_i psi(x + e_i) e_i

where the sum is a finite difference approximation for grad(psi). The problem now is that you can find two sets of weighting coefficients w_i in the literature. Take the following discretisation as an example:

e_i = (1,0) (0,1) (-1,0) (0,-1) (1,1) (-1,1) (-1,-1) (1,-1)

then usually the corresponding weights are

1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36

However, these actually contain the squared lattice speed of sound! In the derivation of the finite difference stencils, the weights are usually normalized by the squared speed of sound, i.e.,

1/3, 1/3, 1/3, 1/3, 1/12, 1/12, 1/12, 1/12

(For the derivation, see e.g. Shan 2006 Phys Rev E 73, 047701 (2006) or Sbragaglia et al. Phys Rev E 75, 026702 (2007).
Depending on the weights you use, the weighting coefficients actually would have units.
Normalized -> w_i is dimensionless
Not normalized -> w_i has dimensions L^2*T^-2

Whats the consequence of this? If you use NORMALIZED weights, your finite difference formula will approximate

F = -G psi grad(psi)

On the other hand, if you use the non-normalized weights (1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36), you will approximate

F = -G psi c_s^2 grad(psi)

This formula is actually used to derive the pressure equation in your post (p = c_s^2 rho + c_s^2G0.5*psi^2). If I now use the force equation with the squared speed of sound in it to obtain the units of G, I will end up with my result. Using the pressure equation, I will end up at the same result. The problem is that (depending on the weights) there will be a c_s^2 appearing in front of G in the pressure or not. This leads to the inconsistency.

I hope this is understandable.

Regards,

Knut