Dear kk,

I am doing multi component simulation. I was doing it using the simple SC EOS and so i am trying to convert it into a better EOS. I thought it would be best to try it on a simple SCMP problem first and then move on to MCMP. the code i shared from palabos had the same problem i had, just was in SCMP instead of MCMP.

I tried what you asked. Changing G to -1 helped make it run a bit longer but i still cannot observe any phase separation for any value of T. I will give my SCMP code below. i am not able to pinpoint the error.

I adapted the code from another one which i had written some time ago. So ignore the code written for the walls and the wall interaction as i did not place any walls.

clear all

close all

nx=400; %matrix x dimension

ny=400; %matrix y dimension

density=20; %density

t=0.0403;

tr=t/0.0869;

a=2/49;

b=2/21;

alpha=1.7405;

tend=2000;

tau=1;

G=-1;

Gads=-327;

% D2Q9 LATTICE CONSTANTS

to = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];

ex = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];

ey = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];

oppf = [ 1, 4, 5, 2, 3, 8, 9, 6, 7];

f1=3.0;

f2=9/2;

f3=3/2;

f=ones(9,ny,nx); %density distribution function for air

feq=f;

forcing=2.2222e-6;

rho=ones(ny,nx); %density of air

r=zeros(8,ny,nx);

fx=zeros(ny,nx); %intra fluid x force

fy=zeros(ny,nx); %intra fluid y force

fwx=zeros(ny,nx); %intra fluid x force

fwy=zeros(ny,nx); %intra fluid y force

psi=zeros(ny,nx);

solid=zeros(ny,nx);

ux=rho;

uy=rho;

ueqx=rho;

ueqy=rho;

ran=rand([ny nx]);

ran=ran+density;

for i=1:9

f(i,:,:)=to(i).*ran;

end

for i=1:ny

for j=1:nx

if solid(i,j)==1

f(:,i,j)=to(:);

end

end

end

for z=1:tend

disp(z);

%calculating macrovariables

rho(:,:)=sum(f,1);

ux(:,:)=(f(2,:,:)+f(6,:,:)+f(9,:,:)-f(4,:,:)-f(7,:,:)-f(8,:,:));

uy(:,:)=(f(3,:,:)+f(6,:,:)+f(7,:,:)-f(5,:,:)-f(8,:,:)-f(9,:,:));

ux=ux./rho;

uy=uy./rho;

p=(rho.*t)./(1-b.*rho)-(a.*alpha.*rho.^2)./(1+b.**rho);*

psi=sqrt(abs(2.(p-rho./3)./(6G)));

```
fx=to(2)*circshift(psi,[0 -1])+to(6)*circshift(psi,[-1 -1])+to(9)*circshift(psi,[1 -1])-to(4)*circshift(psi,[0 1])-to(8)*circshift(psi,[1 1])-to(7)*circshift(psi,[-1 1]);
fy=to(3)*circshift(psi,[-1 0])+to(6)*circshift(psi,[-1 -1])+to(7)*circshift(psi,[-1 1])-to(5)*circshift(psi,[1 0])-to(8)*circshift(psi,[1 1])-to(9)*circshift(psi,[1 -1]);
```

% fx=to(2).*2.*circshift(psi,[0 -1])+to(6)./2.*circshift(psi,[-1 -1])+to(9)./2.*circshift(psi,[1 -1])-to(4).*2.*circshift(psi,[0 1])-to(8)./2.*circshift(psi,[1 1])-to(7)./2.*circshift(psi,[-1 1]);

% fy=to(3).*2.*circshift(psi,[-1 0])+to(6)./2.*circshift(psi,[-1 -1])+to(7)./2.*circshift(psi,[-1 1])-to(5).*2.*circshift(psi,[1 0])-to(8)./2.*circshift(psi,[1 1])-to(9)./2.*circshift(psi,[1 -1]);

% not sure if it is this one or the other (commented above)

```
fwx=to(2)*circshift(solid,[0 -1])+to(6)*circshift(solid,[-1 -1])+to(9)*circshift(solid,[1 -1])-to(4)*circshift(solid,[0 1])-to(8)*circshift(solid,[1 1])-to(7)*circshift(solid,[-1 1]);
fwy=to(3)*circshift(solid,[-1 0])+to(6)*circshift(solid,[-1 -1])+to(7)*circshift(solid,[-1 1])-to(5)*circshift(solid,[1 0])-to(8)*circshift(solid,[1 1])-to(9)*circshift(solid,[1 -1]);
```

%above wall interaction will be zero as solid=0

```
fx=-6*G.*psi.*fx;
fy=-6*G.*psi.*fy;
fwx=-1*Gads.*psi.*fwx;
fwy=-1*Gads.*psi.*fwy;
ux=ux+(fx+fwx).*tau./rho;
uy=uy+(fy+fwy+forcing).*tau./rho;
% calculating equilibrium conditions and collision
for i=1:9
cuNS = 3*(ex(i)*ux+ey(i)*uy);
feq(i,:,:) = rho .* to(i) .* ...
( 1 + cuNS + 1/2*(cuNS.*cuNS) - 3/2*(ux.^2+uy.^2) );
for j=1:ny
for k=1:nx
if solid(j,k)==1
feq(:,j,k)=f(:,j,k);
end
end
end
f(i,:,:) = f(i,:,:) - (f(i,:,:)-feq(i,:,:))./tau;
end
%solid collision
for i=1:ny
for j=1:nx
if solid(i,j)==1
f(:,i,j)=f(oppf,i,j);
end
end
end
%streaming
for i=1:9
f(i,:,:)=circshift(f(i,:,:),[0,ey(i) ex(i)]);
end
%visualising
u=sqrt(ux.^2+uy.^2);
figure(1);
imagesc(rho);
axis equal off;
drawnow;
```

end

I greatly appreciate your help.

Regards,

Githin