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. / (3nu2+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,lxly), 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,lxly)), 1,lx,ly);
jx2 = reshape ( (cx * reshape(gIn,9,lxly)), 1,lx,ly);
jy2 = reshape ( (cy * reshape(gIn,9,lxly)), 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(rho1t(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) + …
18t(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