How Can I plot Streamlines for the following code. You can reduce the time steps to save your time. Thanks!
8 views (last 30 days)
Show older comments
Mahendra Yadav
on 28 May 2021
Commented: Cris LaPierre
on 30 May 2021
clc
clear all
close all
M = 1000; N = 40;
xE = M/N; yE = N/N;
dx = xE/M; dy = yE/N;
x = 0:dx:xE;
y = 0:dy:yE;
tfinal = 4000;
twall = 1.0;
u0 = 0.2;
nu = 0.02;
Cs = 1/sqrt(3);
omega = 1/(nu/Cs^2 + 0.5);
Re = u0*N/nu;
w=[1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36 4/9];
cx = [1 0 -1 0 1 -1 -1 1 0];
cy = [0 1 0 -1 1 1 -1 -1 0];
f = zeros(M+1,N+1,9);
feq = f;
rho = ones(M+1,N+1);
u = zeros(M+1,N+1);
v = u;
% Velocity Boundary Conditions
u_east = 0.0; v_east = 0.0;
u_west = u0; v_west = 0.0;
u_north = 0.0; v_north = 0.0;
u_south = 0.0; v_south = 0.0;
u(1,:) = u_west; v(1,:) = v_west;
u(M+1,:) = u_east; v(M+1,:) = v_east;
u(:,1) = u_south; v(:,1) = v_south;
u(:,N+1) = u_north; v(:,N+1) = v_north;
% LBM Simulation
tic
for t = 1:tfinal
[f] = collision(M,N,cx,cy,w,omega,rho,f,u,v);
[f] = stream(f);
[f] = BConditions(M,N,f,u_west);
[rho,u,v] = res(M,N,f);
end
toc
% Plotting
for i = 1:M+1
for j = 1:N+1
U(i,j) = sqrt(u(i,j) + v(i,j));
end
end
for j = 1:N+1
u1(j) = u(M/2,j)/u0;
end
for i = 1:M+1
v1(i) = v(i,N/2)/u0;
end
figure(1)
contourf(u')
colorbar;
title('Velocity Distribution')
figure(2)
plot(y,u1,'b');
grid on;
title('Horizontal Velocity Distribution along the x = 0.5')
xlabel('y')
ylabel('u')
% figure(3)
% plot(x,v1,'r');
% grid on;
% title('Vertical Velocity Distribution along the y = 0.5')
% xlabel('x')
% ylabel('v')
% ---------------------------------------------------------------------
function [f] = collision(M,N,cx,cy,w,omega,rho,f,u,v)
for i = 1:M+1
for j = 1:N+1
t1 = u(i,j)*u(i,j) + v(i,j)*v(i,j);
for k = 1:9
t2 = cx(k)*u(i,j) + cy(k)*v(i,j);
feq(i,j,k) = w(k)*(rho(i,j))*(1.0 + 3.0*t2 + (9/2)*t2^2 - (3/2)*t1);
f(i,j,k) = (1-omega)*f(i,j,k) + omega*feq(i,j,k);
end
end
end
end
%----------------------------------------------------------------------
function [f] = stream(f)
f(:,:,1)=circshift( squeeze(f(:,:,1)), [+1,+0] );
f(:,:,2)=circshift( squeeze(f(:,:,2)), [+0,+1] );
f(:,:,3)=circshift( squeeze(f(:,:,3)), [-1,+0] );
f(:,:,4)=circshift( squeeze(f(:,:,4)), [+0,-1] );
f(:,:,5)=circshift( squeeze(f(:,:,5)), [+1,+1] );
f(:,:,6)=circshift( squeeze(f(:,:,6)), [-1,+1] );
f(:,:,7)=circshift( squeeze(f(:,:,7)), [-1,-1] );
f(:,:,8)=circshift( squeeze(f(:,:,8)), [+1,-1] );
end
%----------------------------------------------------------------------
function [f] = BConditions(M,N,f,u_west)
% Top Wall and Bottom Wall
for i = 1:M+1
f(i,N+1,4) = f(i,N+1,2);
f(i,N+1,7) = f(i,N+1,5);
f(i,N+1,8) = f(i,N+1,6);
f(i,1,2) = f(i,1,4);
f(i,1,5) = f(i,1,7);
f(i,1,6) = f(i,1,8);
end
% Right Wall
for j = 2:N
f(N+1,j,3) = f(N,j,3);
f(N+1,j,6) = f(N,j,6);
f(N+1,j,7) = f(N,j,7);
end
% Left Wall
for j = 2:N
rhow(1,j) = (f(1,j,9) + f(1,j,2) + f(1,j,4) + 2*(f(1,j,3) + f(1,j,6) + f(1,j,7)))/(1-u_west);
f(1,j,1) = f(1,j,3) + (2/3)*rhow(1,j)*u_west;
f(1,j,5) = f(1,j,7) + rhow(1,j)*u_west/6;
f(1,j,8) = f(1,j,6) + rhow(1,j)*u_west/6;
end
end
%-------------------------------------------------------------------------------
function [rho,u,v] = res(M,N,f)
rho = sum(f,3);
% for j = 1:N+1
% rho(1,j) = (f(1,j,9) + f(1,j,2) + f(1,j,4) + 2*(f(1,j,3) + f(1,j,6) + f(1,j,7)))/(1-0.2);
% end
u = (sum(f(:,:,[1 5 8]),3) - sum(f(:,:,[3 6 7]),3))./rho;
v = (sum(f(:,:,[2 5 6]),3) - sum(f(:,:,[4 7 8]),3))./rho;
end
3 Comments
Cris LaPierre
on 29 May 2021
The code you shared does not use the streamline function. Share the code that is giving you an error, as well as the full text of that error.
Accepted Answer
Cris LaPierre
on 29 May 2021
Edited: Cris LaPierre
on 30 May 2021
To use streamlines, you need a grid of x and y vertices, a grid of u and v vector components, and 'seed' x and y values for where the streamline(s) should start. Your code creates x,y,u and v variables, so here's what streamlines might look like. Note that since u and v are small, the streamlines don't show much deviation from horizontal.
Since it takes almost an hour for your code to run, I'm just loading the final variables for this example.
load MYvars
% create grid
[X,Y] = meshgrid(y,x);
% Quiver plot to show vector field
quiver(X,Y,u,v)
% Now show streamlines
figure
streamline(X,Y,u,v,ones(1,6)*0.1,[0 5 10 15 20 25])
3 Comments
Cris LaPierre
on 30 May 2021
You can use non-uniform grids to create X and Y. However, X and Y need to be matrices of the same size, and U and V also need to be the same size and X and Y. Your startx and starty must be the same size, but this can be a single point, a vector, or a matrix. This is just the starting point(s) of a streamline. The actual line is drawn automatically following the vectors created by U and V.
Perhaps these answers may be insightful
More Answers (0)
See Also
Categories
Find more on Axis Labels 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!