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
