Clear Filters
Clear Filters

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

5 views (last 30 days)
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

Answers (1)

Anurag Ojha
Anurag Ojha on 3 Oct 2024 at 8:34
Hey
To represent the three parameters (SOLID, rho1, and rho2) with distinct colors in the same domain, you can combine them into a single array where each parameter corresponds to a different channel in the RGB color model. For instance:
  • Use SOLID as the red channel.
  • Use rho1 as the green channel.
  • Use rho2 as the blue channel.
This way, areas of the plot where all three variables are significant will appear as a combination of their respective colors. Below is a simplified code which I have written taking certain assumption. Kindly modify it as per your use case:
% Clear workspace and close all figures
clear;
clc;
close all;
% Define numerical parameters
N = 100; % Grid size
nx = N;
ny = N;
% Initialize solid capturing: Randomly assign solid areas
SOLID = rand(nx, ny) > 0.7; % Extremely porous random domain
% Add boundaries as solid
SOLID(1, :) = true;
SOLID(end, :) = true;
SOLID(:, 1) = true;
SOLID(:, end) = true;
% Initialize densities for fluid 1 (rho1) and fluid 2 (rho2)
delta_rho = 0.01 * (1 - 2 * rand(ny, nx)); % Small random perturbations
rho1 = 1 + delta_rho; % Fluid 1 density
rho2 = 1 - delta_rho; % Fluid 2 density
% Normalize the values for visualization
rho1_norm = (rho1 - min(rho1(:))) / (max(rho1(:)) - min(rho1(:)));
rho2_norm = (rho2 - min(rho2(:))) / (max(rho2(:)) - min(rho2(:)));
SOLID_norm = double(SOLID); % Ensure SOLID is numeric
% Combine the parameters into an RGB image
% Red for SOLID, Green for rho1, Blue for rho2
RGB_image = cat(3, SOLID_norm, rho1_norm, rho2_norm);
% Display the RGB image
figure;
imagesc(RGB_image);
axis equal off;
title('Representation of SOLID, rho1, and rho2 in RGB');

Categories

Find more on Fluid Dynamics in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!