How Can I plot Streamlines for the following code. You can reduce the time steps to save your time. Thanks!

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

[Cris LaPierre]

I would use the streamline function.

[Mahendra Yadav]

That's all I'm asking. While using Streamline function I'm getting error. If possible could you give me correct command.

[Cris LaPierre]

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.

Answer 1

[Cris LaPierre]

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])

[Mahendra Yadav]

    >> startx = 0:0.1:1;
    >> starty = 0:0.1:1;
    >> streamline(x,y,v,u,startx,starty)
    Error using stream2 (line 46)
    The size of X must match the size of U or the number of columns of U.
    
    Error in streamline (line 62)
          verts = stream2(x,y,u,v,sx,sy,options);
     

[Mahendra Yadav]

I want to say that you defined another mesh grid to use streamline function.

But could you tell me how can I plot streamlines with my orginal dimensions i.e. with non - uniform meshgrids. Please use timesteps = 500 in the code to save your time.

[Cris LaPierre]

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

• https://www.mathworks.com/matlabcentral/answers/762791-streamlines-do-not-plot?s_tid=srchtitle
• https://www.mathworks.com/matlabcentral/answers/790449-how-can-i-plot-streamline-from-a-stream-function