MATLAB Answers

Colour representation of three matrices entities (SOLID, rho1 and rho2)

11 views (last 30 days)
Oluwaseyi Aliu
Oluwaseyi Aliu on 22 Jan 2021
Dear friend,
Please I am solving a mixture of two fluids in a solid domain just as presented in the code below; I want a representation of the three parameters(SOLID,rho1 and rho2) with distinct colours in the same domain. Attached are two graphs (the one with two colours is the result I got, the result with three colours is what I want).
clear clc;close all;
% define numerical parameters
%new BounceBack<T,Descriptor>(0.5);
N=256;
nx=1*N; ny=N;
tau1=1; tau2=1; % relexation time
G=-1.5;
% define weight coefficient(D2Q9)
w0=4/9;
w1=1/9; w2=1/9; w3=1/9; w4=1/9;
w5=1/36; w6=1/36; w7=1/36; w8=1/36;
% initialize variable values in the field
c=1; %lattice speed
dt=1;% delta t
% solid capturing initialisation
SOLID=rand(nx+1,ny+1)>0.7;%extremely porous random domain
%SOLID = false( 1, nx, ny );
SOLID( 1, : ) = true;
SOLID( end, : ) = true;
SOLID( :, 1 ) = true;
SOLID( :, end ) = true;
%SOLID = find( SOLID );
% initialize distribution functions for two components
delta_rho=0.001*(1-2*rand(ny+1,nx+1));
rho1=1+delta_rho;
rho2=1-delta_rho;
% distribution function for component 1
f(:,:,1)=w1*rho1;
f(:,:,2)=w2*rho1;
f(:,:,3)=w3*rho1;
f(:,:,4)=w4*rho1;
f(:,:,5)=w5*rho1;
f(:,:,6)=w6*rho1;
f(:,:,7)=w7*rho1;
f(:,:,8)=w8*rho1;
f(:,:,9)=w0*rho1;
% distribution function for component 2
g(:,:,1)=w1*rho2;
g(:,:,2)=w2*rho2;
g(:,:,3)=w3*rho2;
g(:,:,4)=w4*rho2;
g(:,:,5)=w5*rho2;
g(:,:,6)=w6*rho2;
g(:,:,7)=w7*rho2;
g(:,:,8)=w8*rho2;
g(:,:,9)=w0*rho2;
for it=1:1000
% macropic properties
% calculate interaction body forces
rho1=sum(f,3); %density of fluid 1
rho2=sum(g,3); %density of fluid 2
rho_tot=rho1+rho2;% total local density
% body forces for conhension
F221_x=-rho1.*G.*(w1*circshift(rho2,[0,1])-w3*circshift(rho2,[0,-1])+w5*circshift(rho2,[-1,1])-w6*circshift(rho2,[-1,-1])...
-w7*circshift(rho2,[1,-1])+w8*circshift(rho2,[1,1]));
F221_y=-rho1.*G.*(w2*circshift(rho2,[-1 0])-w4*circshift(rho2,[1,0])+w5*circshift(rho2,[-1,1])+w6*circshift(rho2,[-1,-1])...
-w7*circshift(rho2,[1,-1])-w8*circshift(rho2,[1,1]));
F122_x=-rho2.*G.*(w1*circshift(rho1,[0,1])-w3*circshift(rho1,[0,-1])+w5*circshift(rho1,[-1,1])-w6*circshift(rho1,[-1,-1])...
-w7*circshift(rho1,[1,-1])+w8*circshift(rho1,[1,1]));
F122_y=-rho2.*G.*(w2*circshift(rho1,[-1 0])-w4*circshift(rho1,[1,0])+w5*circshift(rho1,[-1,1])+w6*circshift(rho1,[-1,-1])...
-w7*circshift(rho1,[1,-1])-w8*circshift(rho1,[1,1]));
% velocity field
u=(sum(f(:,:,[1 5 8]),3)-sum(f(:,:,[3 6 7]),3)+sum(g(:,:,[1 5 8]),3)-sum(g(:,:,[3 6 7]),3))./rho_tot+(F221_x+F122_x)./rho_tot./2;
v=(sum(f(:,:,[2 5 6]),3)-sum(f(:,:,[4 7 8]),3)+sum(g(:,:,[2 5 6]),3)-sum(g(:,:,[4 7 8]),3))./rho_tot+(F221_y+F122_y)./rho_tot./2;
% collision step
% calculate equilibrium distribution for fluid 1
feq(:,:,1)=w1*rho1.*(1+3*u/c+9/2*u.^2/c^2-3/2*(u.^2+v.^2)/c^2);
feq(:,:,2)=w2*rho1.*(1+3*v/c+9/2*v.^2/c^2-3/2*(u.^2+v.^2)/c^2);
feq(:,:,3)=w3*rho1.*(1+3*-u/c+9/2*u.^2/c^2-3/2*(u.^2+v.^2)/c^2);
feq(:,:,4)=w4*rho1.*(1+3*-v/c+9/2*v.^2/c^2-3/2*(u.^2+v.^2)/c^2);
feq(:,:,5)=w5*rho1.*(1+3*(u+v)/c+9/2*(u+v).^2/c^2-3/2*(u.^2+v.^2)/c^2);
feq(:,:,6)=w6*rho1.*(1+3*(-u+v)/c+9/2*(-u+v).^2/c^2-3/2*(u.^2+v.^2)/c^2);
feq(:,:,7)=w7*rho1.*(1+3*(-u-v)/c+9/2*(-u-v).^2/c^2-3/2*(u.^2+v.^2)/c^2);
feq(:,:,8)=w8*rho1.*(1+3*(u-v)/c+9/2*(u-v).^2/c^2-3/2*(u.^2+v.^2)/c^2);
feq(:,:,9)=w0*rho1.*(1-3/2*(u.^2+v.^2)/c^2);
% for fluid 2
geq(:,:,1)=w1*rho2.*(1+3*u/c+9/2*u.^2/c^2-3/2*(u.^2+v.^2)/c^2);
geq(:,:,2)=w2*rho2.*(1+3*v/c+9/2*v.^2/c^2-3/2*(u.^2+v.^2)/c^2);
geq(:,:,3)=w3*rho2.*(1+3*-u/c+9/2*u.^2/c^2-3/2*(u.^2+v.^2)/c^2);
geq(:,:,4)=w4*rho2.*(1+3*-v/c+9/2*v.^2/c^2-3/2*(u.^2+v.^2)/c^2);
geq(:,:,5)=w5*rho2.*(1+3*(u+v)/c+9/2*(u+v).^2/c^2-3/2*(u.^2+v.^2)/c^2);
geq(:,:,6)=w6*rho2.*(1+3*(-u+v)/c+9/2*(-u+v).^2/c^2-3/2*(u.^2+v.^2)/c^2);
geq(:,:,7)=w7*rho2.*(1+3*(-u-v)/c+9/2*(-u-v).^2/c^2-3/2*(u.^2+v.^2)/c^2);
geq(:,:,8)=w8*rho2.*(1+3*(u-v)/c+9/2*(u-v).^2/c^2-3/2*(u.^2+v.^2)/c^2);
geq(:,:,9)=w0*rho2.*(1-3/2*(u.^2+v.^2)/c^2);
%calculate body force terms in collision step
% for fluid 1
F1(:,:,1)=w1*(1-1/2/tau1)*((6*u+3).*(F221_x)+-3*v.*(F221_y));
F1(:,:,2)=w2*(1-1/2/tau1)*(-3*u.*(F221_x)+(3+6*v).*(F221_y));
F1(:,:,3)=w3*(1-1/2/tau1)*((6*u-3).*(F221_x)+-3*v.*(F221_y));
F1(:,:,4)=w4*(1-1/2/tau1)*(-3*u.*(F221_x)+(-3+6*v).*(F221_y));
F1(:,:,5)=w5*(1-1/2/tau1)*((3+6*u+9*v).*(F221_x)+(3+9*u+6*v).*(F221_y));
F1(:,:,6)=w6*(1-1/2/tau1)*((-3+6*u-9*v).*(F221_x)+(3-9*u+6*v).*(F221_y));
F1(:,:,7)=w7*(1-1/2/tau1)*((-3+6*u+9*v).*(F221_x)+(-3+9*u+6*v).*(F221_y));
F1(:,:,8)=w8*(1-1/2/tau1)*((3+6*u-9*v).*(F221_x)+(-3-9*u+6*v).*(F221_y));
F1(:,:,9)=w0*(1-1/2/tau1)*(-3*u.*(F221_x)+-3*v.*(F221_y));
% for fluid 2
F2(:,:,1)=w1*(1-1/2/tau2)*((6*u+3).*(F122_x)+-3*v.*(F122_y));
F2(:,:,2)=w2*(1-1/2/tau2)*(-3*u.*(F122_x)+(3+6*v).*(F122_y));
F2(:,:,3)=w3*(1-1/2/tau2)*((6*u-3).*(F122_x)+-3*v.*(F122_y));
F2(:,:,4)=w4*(1-1/2/tau2)*(-3*u.*(F122_x)+(-3+6*v).*(F122_y));
F2(:,:,5)=w5*(1-1/2/tau2)*((3+6*u+9*v).*(F122_x)+(3+9*u+6*v).*(F122_y));
F2(:,:,6)=w6*(1-1/2/tau2)*((-3+6*u-9*v).*(F122_x)+(3-9*u+6*v).*(F122_y));
F2(:,:,7)=w7*(1-1/2/tau2)*((-3+6*u+9*v).*(F122_x)+(-3+9*u+6*v).*(F122_y));
F2(:,:,8)=w8*(1-1/2/tau2)*((3+6*u-9*v).*(F122_x)+(-3-9*u+6*v).*(F122_y));
F2(:,:,9)=w0*(1-1/2/tau2)*(-3*u.*(F122_x)+-3*v.*(F122_y));
% collision
f=f-1/tau1.*(f-feq)+F1;
g=g-1/tau2.*(g-geq)+F2;
%streaming
f(:,:,1)=circshift(f(:,:,1),[0,1]);
f(:,:,2)=circshift(f(:,:,2),[-1,0]);
f(:,:,3)=circshift(f(:,:,3),[0,-1]);
f(:,:,4)=circshift(f(:,:,4),[1,0]);
f(:,:,5)=circshift(f(:,:,5),[-1,1]);
f(:,:,6)=circshift(f(:,:,6),[-1,-1]);
f(:,:,7)=circshift(f(:,:,7),[1,-1]);
f(:,:,8)=circshift(f(:,:,8),[1,1]);
g(:,:,1)=circshift(g(:,:,1),[0,1]);
g(:,:,2)=circshift(g(:,:,2),[-1,0]);
g(:,:,3)=circshift(g(:,:,3),[0,-1]);
g(:,:,4)=circshift(g(:,:,4),[1,0]);
g(:,:,5)=circshift(g(:,:,5),[-1,1]);
g(:,:,6)=circshift(g(:,:,6),[-1,-1]);
g(:,:,7)=circshift(g(:,:,7),[1,-1]);
g(:,:,8)=circshift(g(:,:,8),[1,1]);
% add BC
%L
f(:,1,1)=f(:,nx+1,1);
f(:,1,5)=f(:,nx+1,5);
f(:,1,8)=f(:,nx+1,8);
g(:,1,1)=g(:,nx+1,1);
g(:,1,5)=g(:,nx+1,5);
g(:,1,8)=g(:,nx+1,8);
%R
f(:,nx+1,3)=f(:,1,3);
f(:,nx+1,6)=f(:,1,6);
f(:,nx+1,7)=f(:,1,7);
g(:,nx+1,3)=g(:,1,3);
g(:,nx+1,6)=g(:,1,6);
g(:,nx+1,7)=g(:,1,7);
%T
f(1,:,4)=f(ny+1,:,4);
f(1,:,7)=f(ny+1,:,7);
f(1,:,8)=f(ny+1,:,8);
g(1,:,4)=g(ny+1,:,4);
g(1,:,7)=g(ny+1,:,7);
g(1,:,8)=g(ny+1,:,8);
%B
f(ny+1,:,2)=f(1,:,2);
f(ny+1,:,5)=f(1,:,5);
f(ny+1,:,6)=f(1,:,6);
g(ny+1,:,2)=g(1,:,2);
g(ny+1,:,5)=g(1,:,5);
g(ny+1,:,6)=g(1,:,6);
if rem(it,5)==0
imagesc(rho1);
%imagesc(rho2);
%imagesc(solid);
fff=getframe;
axis equal off;
end
end

  0 Comments

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!