HSD model for multiphase flow

Hi everybody,

	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.

Thanks again.


Dear Andrea:
I worked with the HSD model for the 2D 2phase poiseuille flow simulation. In my experience, the intializtions of phi and rho are very important, it determines whether your result is good, and the gradient of psi plays the most important role in this model. In the places where are far from the interfaces, the gradient of rho(phi) being approximately zero makes sure that the density profile in these places is flat; and the gradient of psi at the interfaces makes sure that you can keep a sharp interface. According to your questions, I have some suggestions, which I hope might help you. First, you could examine your equation of state, and Carnahan and Starling equation of state could not be the only choices. Second, you should try your initializations of rho and phi for more times according to your equation of state, and make sure that the low and high values of them locate on different zones in the equation of state. Third, given the importance of the gradient of psi, I suggest you implement some accurate algorithms to simulate it.
It cost me lots of time to simulate this model, and I know it was very hard. I hope my suggestions might be helpful and wish you good luck~~

Thank you very much cubicmatrixist,

      By the way .... Could you please tell me if the manner I use to  initialize the pressure field looks fine to you?  I mean I first define the rho for the two fluids and I define the initial rho field. From rho field I find the index field phi ... and from phi, using a EOS (in my case Carnahan and Starling), I find the pressure field p. Now I use the macroscopic fields to initialize the particle distribution function at the equilibrium.

Is my way of implementing the initial condition close to the one that you usually implement?

thank you very much…