# problem to implement multiphase fluid flow

We have try to implement a multiphase fluid flow with Matlab, but this code involves look like a Single Phase fluid flow.
We think that the problem is on the initialization of the two fluid:

clear

% GENERAL FLOW CONSTANTS
lx = 200;
ly = 42;
uMax1 = 0.02;
uMax2 = 0.015;
nu1 = 0.4;
nu2 = 0.06;
omega1 = 1. / (3nu1+1./2.);
omega2 = 1. / (3
nu2+1./2.);
maxT = 100000;
tPlot = 5;
G = -1.9;
% D2Q9 LATTICE CONSTANTS
t = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cx = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
cy = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
opp = [ 1, 4, 5, 2, 3, 8, 9, 6, 7];
col = [22:ly-1];
lid=[21:40];
lid1=[1:20];

obst=zeros(lx,ly);
obst([1:20,41:lx],[1:ly-21,ly]) = 1;
obst([21:40],ly)=1;
bbRegion = find(obst);

Gomega1 = G/omega1;
Gomega2 = G/omega2;

fIn = reshape( t’ * ones(1,lxly), 9, lx, ly);
gIn = reshape(t’ * zeros(1,lx
ly), 9, lx, ly);

% MAIN LOOP (TIME CYCLES)
for cycle = 1:maxT
% MACROSCOPIC VARIABLES
rho1 = sum(fIn);
rho2 = sum(gIn);
jx1 = reshape ( (cx * reshape(fIn,9,lxly)), 1,lx,ly);
jy1 = reshape ( (cy * reshape(fIn,9,lx
ly)), 1,lx,ly);
jx2 = reshape ( (cx * reshape(gIn,9,lxly)), 1,lx,ly);
jy2 = reshape ( (cy * reshape(gIn,9,lx
ly)), 1,lx,ly);
rhoTot_OMEGA = rho1omega1 + rho2omega2;
uTotX = (jx1omega1+jx2omega2) ./ rhoTot_OMEGA;
uTotY = (jy1omega1+jy2omega2) ./ rhoTot_OMEGA;
rhoContrib1x = 0.0;
rhoContrib2x = 0.0;
rhoContrib1y = 0.0;
rhoContrib2y = 0.0;
for i=2:9
rhoContrib1x = rhoContrib1x + circshift(rho1*t(i), [0,cx(i),cy(i)])cx(i);
rhoContrib1y = rhoContrib1y + circshift(rho1
t(i), [0,cx(i),cy(i)])*cy(i);

``````    rhoContrib2x = rhoContrib2x + circshift(rho2*t(i), [0,cx(i),cy(i)])*cx(i);
rhoContrib2y = rhoContrib2y + circshift(rho2*t(i), [0,cx(i),cy(i)])*cy(i);
end
uTotX1 = uTotX - Gomega1.*rhoContrib2x;
uTotY1 = uTotY - Gomega1.*rhoContrib2y;

uTotX2 = uTotX - Gomega2.*rhoContrib1x;
uTotY2 = uTotY - Gomega2.*rhoContrib1y;

%Inlet
L=18;y=lid1-1.5;
uTotX1(:,1,col) = 20*uMax1 / (L*L) * (y.*L-y.*y);
uTotY1(:,1,col) = 0;
rho1(:,1,col) = 1 ./ (1-uTotX1(:,1,col)) .* ( ...
sum(fIn([1,3,5],1,col)) + ...
2*sum(fIn([4,7,8],1,col)) );

uTotX2(:,lid,1)= 0;
uTotY2(:,lid,1)= 25*uMax2;
rho2(:,lid,1) = 1 ./ (1-uTotY2(:,lid,1)) .* ( ...
sum(gIn([1,2,4],lid,1)) + 2*sum(gIn([5,9,8],lid,1)) );
%Outlet:
rho1(:,lx,col) = rho1(:,lx-1,col);
uTotY1(:,lx,col)  = 0;
uTotX1(:,lx,col)  = uTotX1(:,lx-1,col);
rho2(:,lx,col) = rho2(:,lx-1,col);
uTotY2(:,lx,col)  = 0;
uTotX2(:,lx,col)  = uTotX2(:,lx-1,col);
``````

% COLLISION STEP
for i=1:9
cu1 = 3*(cx(i)uTotX1+cy(i)uTotY1);
cu2 = 3
(cx(i)uTotX2+cy(i)uTotY2);
fEq(i,:, = rho1 .
t(i) .

( 1 + cu1 + 1/2
(cu1.cu1) …
- 3/2
(uTotX1.^2+uTotY1.^2) );
fOut(i,:, = fIn(i,:, - …
omega1 .* (fIn(i,:,:)-fEq(i,:,:));

``````   gEq(i,:,:)  = rho2 .* t(i) .* ...
( 1 + cu2 + 1/2*(cu2.*cu2) ...
- 3/2*(uTotX2.^2+uTotY2.^2) );
gOut(i,:,:) = gIn(i,:,:) - ...
omega2 .* (gIn(i,:,:)-gEq(i,:,:));
end

% MICROSCOPIC BOUNDARY CONDITIONS
for i=1:9

fOut(i,1,col) = fEq(i,1,col) + ...
18*t(i)*cx(i)*cy(i)* ( fIn(8,1,col) - ...
fIn(7,1,col)-fEq(8,1,col)+fEq(7,1,col) );
``````

%
fOut(i,lid,1) = fEq(i,lid,1) + …
18*t(i)cx(i)cy(i) ( fIn(8,lid,1) - …
fIn(9,lid,1)-fEq(8,lid,1)+fEq(9,lid,1) );
gOut(i,lid,1) = gEq(i,lid,1) + …
18
t(i)*cx(i)cy(i) ( gIn(8,lid,1) - …
gIn(9,lid,1)-gEq(8,lid,1)+gEq(9,lid,1) );

``````       fOut(i,lx,col) = fEq(i,lx,col) + ...
18*t(i)*cx(i)*cy(i)* ( fIn(6,lx,col) - ...
fIn(9,lx,col)-fEq(6,lx,col)+fEq(9,lx,col) );
gOut(i,lx,col) = gEq(i,lx,col) + ...
18*t(i)*cx(i)*cy(i)* ( gIn(6,lx,col) - ...
gIn(9,lx,col)-gEq(6,lx,col)+gEq(9,lx,col) );

fOut(i,bbRegion) = fIn(opp(i),bbRegion);
gOut(i,bbRegion) = gIn(opp(i),bbRegion);
end
% STREAMING STEP
for i=1:9
fIn(i,:,:) = ...
circshift(fOut(i,:,:), [0,cx(i),cy(i)]);
gIn(i,:,:) = ...
circshift(gOut(i,:,:), [0,cx(i),cy(i)]);
end

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

end
end

I see!

it’s shan-chen model.

Hi Gabriele,

I’m M.S.c student in Iran and work in on LBM,

I have some experience in flow past a rotating cylinder with LBM.

Now I want to start 2-phase flows simulation with LBM.