Shan Chen Model for High Density Ratio Multi Component Flow

Palabos General
#1

Hello everyone,

I’m currently programming a multicomponent Lattice Boltzmann model using Matlab. The problem I keep encountering is that at high density ratios, the simulation becomes unstable and yields weird results. I am using the Shan Chen model for the intermolecular forces.

Can any one help me in this regard and explain why the problem is occuring, especially at high density ratios? I’ve heard that to model meta stable states you need to take two forces into account - short range and mid range. Is that the problem here?

I’m relatively new to the Lattice Boltzmann method, and any help would be greatly appreciated.

#2

Here is the program I’m using. It’s a simple program, and it is a modification of the 2D Segregation program available in the LBM wiki.

It initializes two fluids, with the Fluid 1 concentrated at the centre in a circle (of 25 lu radius) and the Fluid 2 surrounding it in a 200X200 box.

Fluid 1 is subjected to gravitational force in the downward direction.


clear
clf

% GENERAL FLOW CONSTANTS

ly           = 201;
lx           = 201;

G = -4;  % Amplitude of the molecular interaction force

omega1 = 1.;  % Relaxation parameter for fluid 1
omega2 = 1.;  % Relaxation parameter for fluid 2

maxT   = 80000;    % total number of iterations
tPlot  = 10;       % iterations between successive graphical outputs

% D2Q9 LATTICE CONSTANTS
tNS   = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cxNS  = [  0,   1,  0, -1,  0,    1,  -1,  -1,   1];
cyNS  = [  0,   0,  1,  0, -1,    1,   1,  -1,  -1];
oppNS = [  1,   4,  5,  2,  3,    8,   9,   6,   7];

[y,x] = meshgrid(1:ly,1:lx);

rho1 = zeros(lx,ly);
rho2 = ones(lx,ly);
for i=1:lx
    for j=1:ly
        if ((i-100)^2 + (j-100)^2 < 625)
            rho1(i,j) = 2;
            rho2(i,j) = 0;
        end
    end
end

% INITIAL CONDITION FOR BOTH DISTRIBUTION FUNCTIONS: (T=0) ==> TIn(i) = t(i)
for i=1:9
    fIn(i,1:lx,1:ly) = tNS(i).*rho1;
    gIn(i,1:lx,1:ly) = tNS(i).*rho2;
end

rho1  = reshape(sum(fIn),lx,ly);
imagesc(rho1');
colorbar
title('Fluid 1 density');
axis equal off; drawnow

% MAIN LOOP (TIME CYCLES)
Gomega1 = G/omega1;
Gomega2 = G/omega2;
for cycle = 1:maxT
    display(cycle);
    % MACROSCOPIC VARIABLES
    rho1 = sum(fIn);
    rho2 = sum(gIn);
    jx1  = reshape ( (cxNS * reshape(fIn,9,lx*ly)), 1,lx,ly);
    jy1  = reshape ( (cyNS * reshape(fIn,9,lx*ly)), 1,lx,ly);
    jx2  = reshape ( (cxNS * reshape(gIn,9,lx*ly)), 1,lx,ly);
    jy2  = reshape ( (cyNS * reshape(gIn,9,lx*ly)), 1,lx,ly);
   
    rhoTot_OMEGA = rho1*omega1 + rho2*omega2;
    uTotX = (jx1*omega1+jx2*omega2) ./ rhoTot_OMEGA;
    uTotY = (jy1*omega1+jy2*omega2) ./ rhoTot_OMEGA;
	
    rhoContrib1x = 0.0;
    rhoContrib2x = 0.0;
    
    rhoContrib1y = 0.0;
    rhoContrib2y = 0.0;
    
    rho1change = (1 - exp(-rho1));
    rho2change = (1 - exp(-rho2));
    for i=2:9
        rhoContrib1x = rhoContrib1x + circshift(rho1change*tNS(i), [0,cxNS(i),cyNS(i)])*cxNS(i);
        rhoContrib1y = rhoContrib1y + circshift(rho1change*tNS(i), [0,cxNS(i),cyNS(i)])*cyNS(i);
        
        rhoContrib2x = rhoContrib2x + circshift(rho2change*tNS(i), [0,cxNS(i),cyNS(i)])*cxNS(i);
        rhoContrib2y = rhoContrib2y + circshift(rho2change*tNS(i), [0,cxNS(i),cyNS(i)])*cyNS(i);
    end
    
    uTotX1 = uTotX - Gomega1.*rhoContrib2x; %POTENTIAL CONTRIBUTION OF FLUID 2 ON 1
    uTotY1 = uTotY - Gomega1.*rhoContrib2y +0.001/omega1;
    
    uTotX2 = uTotX - Gomega2.*rhoContrib1x; %POTENTIAL CONTRIBUTION OF FLUID 2 ON 1
    uTotY2 = uTotY - Gomega2.*rhoContrib1y;

   % COLLISION STEP FLUID 1 AND 2
   for i=1:9
      cuNS1        = 3*(cxNS(i)*uTotX1+cyNS(i)*uTotY1);
      cuNS2        = 3*(cxNS(i)*uTotX2+cyNS(i)*uTotY2);
      
      fEq(i,:,:)   = rho1 .* tNS(i) .* ...
                       ( 1 + cuNS1 + 0.5*(cuNS1.*cuNS1) - 1.5*(uTotX1.^2+uTotY1.^2) );
                       
      gEq(i,:,:)   = rho2 .* tNS(i) .* ...
                       ( 1 + cuNS2 + 0.5*(cuNS2.*cuNS2) - 1.5*(uTotX2.^2+uTotY2.^2) );
                       
      fOut(i,:,:)  = fIn(i,:,:) - omega1 .* (fIn(i,:,:)-fEq(i,:,:));
      gOut(i,:,:)  = gIn(i,:,:) - omega2 .* (gIn(i,:,:)-gEq(i,:,:));
   end
         
   % STREAMING STEP FLUID 1 AND 2
   for i=1:9
      fIn(i,:,:) = circshift(fOut(i,:,:), [0,cxNS(i),cyNS(i)]);
      gIn(i,:,:) = circshift(gOut(i,:,:), [0,cxNS(i),cyNS(i)]);
   end

   % VISUALIZATION
   if(mod(cycle,tPlot)==0)
       rho1 = reshape(rho1,lx,ly);
       img = imagesc(rho1'); 
       colorbar
       title('Fluid 1 density');
       axis equal off; drawnow
   end
end

I am trying to understand multi component flows and why this code won’t work for higher density ratios (rho1/rho2). Been trying everything I can all day, and still no results.

This is probably the only forum I know which deals in Lattice Boltzmann problems. Any help would be greatly appreciated.

#3

hi my friend
it is not your code problem becuase shan chen model in very sensitive to high density.
As you know when density chenges from 1 to 1000 there is large gradient is the domain so that"s why your code becomes unstabe.
there are several scheme you can use to simulate. for example use EOS in your model.
best regards

#4

Hi to all!
I want to simulate mixing of two fluids (like water and ethanol) in a T-type micromixer using LBM. I have some question and i would appreciate anyone who can help me.
1-Am i supposed to solve a multiphase flow?
2-Can i use shan-chen model?
3-Can i use Palabos to solve my problem?