I'm trying to implement the HSD scheme (He, Shan, Doolen) for multiphase fluid flow and I'm having some tiny little trouble during the runs (catastrophic numerical divergence).
I’m playing with the case of a single 2D static bubble placed in the middle of my numerical domain… and, assuming that the algorithm I wrote is correct (I would be surprised if it would… but) I’m wondering if I’m messing something up during the variable initialization step.
I’ll try to briefly summarize what I’m doing. The HSD scheme requires to initialize the following macroscopic quantities:
. scalar density field rho;
. scalar phase index phi;
. scalar pressure field p;
. vector velocity (momentum) field u;
Now I define the densities for the heavy (rho_h=0.5) and light phase (rho_l=0.1), the extreme indexes for phi between which we can have phase segregation (when intermolecular forces are strong enough to keep the two phases separated) and I chose for the heavy phase phi_h=0.2589 and light phase phi_l= 0.0403 ( these values are taken from Quingming Chang PhD thesis… but in the He 1999 paper is explained how to find them ).
Now, given that I know the relation between the index phase phi(x) and the density rho(x):
rho(x)= rho_l + (phi(x)-phi_l)*(rho_h-rho_l)/(phi_h-phi_l), (1)
I first define the density field for the 2D bubble (concentric region at the middle of the numerical domain with rho=rho_l and external region with rho=rho_h ) and then find phi(x) using eq. 1
For the pressure field p I use the Carnahan and Starling equation of state equation of state:
p= cs2rho (1+ Brho/4+ (Brho/4)^2 - (Brho/4)^3 )/(1-Brho/4)^3 + A*rho^2
where cs2 is the square speed of sound of my lattice, B and A are constants. I choose their values as B=4 and A=12cs2 (like in Chang thesis or He1999 paper).
Moreover I set the initial velocity equal to zero.
From the definitions of rho,p,u,phi I initialize the 2 particle distribution function fin and gin at the equilibrium state.
Then I let run the code … but given that after a few iteration my bubble blow up, I thin I’m making something wrong …
One more thing. Let’s suppose that I initialize the problem in a correct way… well what are the tricky steps of the HSD scheme ? For example, in the evolution eqs. for fin and gin, we are required the calculate the gradient of the potential psi. However for the particle distribution function f_i it seems we need the gradient of psi where psi is a function of phi and in the second distribution function g_i we need the gradient of psi as a function of phi. I thought that the two were the same … but maybe they are not. There’s somebody there who has any clue about that?
thank you very much for your help.
Any kind hints would be great, even if you are implementing others problems than bubbles or drops (maybe some of you play with Rayleigh-Taylor…) it would be great to have some hints.