How to convert plot a contour plot or streamline plot with this code?
6 views (last 30 days)
Show older comments
% ***** PROGRAM 2.7 *****
% CHANNEL FLOW PAST A RECTANGULAR CYLINDER, AN EXAMPLE ON CURVED SURFACE AND
% THE DERIVATIVE BOUNDARY CONDITION. PLOTS OBTAINED USING THE SEARCH PROGRAM
clc;
clf;
clear all;
format short e;
global M N H JFIRST;
M = 13;
N = 5;
SYMBOL = zeros(11,1);
CONST = zeros(11,1);
JFIRST = zeros(M,1);
ID = zeros(M,N);
PSI = zeros(M,N);
PSIF = zeros(M,N);
A = zeros(M,N);
B = zeros(M,N);
C = zeros(M,N);
D = zeros(M,N);
X = zeros(M,1);
Y = zeros(N,1);
% ***** ASSIGN INPUT DATA *****
H = 0.25;
SYMBOL = '+osdv^x*ph.';
CONST = [0.0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1.0];
% ***** DEFINE THE VALUE JAT THE LOWEST INTERIOR GRID POINT ALONG EACH
% VERTICAL GRID LINE *****
for I=1:M
JFIRST(I) = 2;
if (I >= 5 && I <= 7)
JFIRST(I) = 4;
end
end
% ***** ASSIGN PSI AT GRID POINTS ON OR INSIDE SOLID BOUNDARIES *****
for I=1:M
J1 = JFIRST(I) - 1;
for J=1:J1
PSI(I,J) = 0.0;
end
end
for I=1:M
PSI(I,N) = 1.0;
end
for J=2:N-1
PSI(1,J) = 1.0;
end
% ***** ASSIGN ID NUMBERS TO ALL INTERIOR POINTS *****
for I=2:M
JA = JFIRST(I);
for J=JA:N-1
ID(I,J) = 0;
if (I == M)
ID(I,J) = -1;
end
end
end
for I=5:7
J = JFIRST(I);
ID(I,J) = 1;
end
ID(1,2) = 2;
ID(1,3) = 2;
ID(1,4) = 2;
ID(1,5) = 2;
% ***** COMPUTE A, B, C AND D FOR THE POINTS AT WHICH (2.11.8) WILL BE
% USED *****
for I=5:7
for J=2:5
if (ID(I,J) == 1)
A(I,J) = H;
B(I,J) = H;
C(I,J) = H;
D(I,J) = H;
if (I == 5)
B(I,J) = B(5,4);
end
if (I == 7)
A(I,J) = A(5,4);
end
if (I >= 5 && I <= 4)
C(I,J) = C(5,7);
end
end
end
end
% ***** ASSIGN GUESSED VALUES TO PSI *****
for I=1:M
JA = JFIRST(I);
for J=JA:N-1
PSI(I,J) = 1;
end
end
% ***** ITERATE UNTIL SUM OF ERROS IS LESS THAN OR EQUAL TO 0.0005 *****
ITER = 0;
ERROR = 1.0;
while (ERROR > 0.0005)
ITER = ITER + 1;
ERROR = 0.0;
for I=1:M
JA = JFIRST(I);
for J=JA:N-1
PSIPRV = PSI(I,J);
if (ID(I,J) == -1)
%(2.11.13) IS USED TO INCORPORATE DERIVATIVE BOUNDARY CONDITION AT THE EXIT:
PSI(I,J) = 0.25 * ( 2.*PSI(I-1,J) + PSI(I,J-1) + PSI(I,J+1) );
elseif (ID(I,J) == 0)
%AT GRID POINTS AWAY FROM THE RECTANGULAR BOUNDARY:
PSI(I,J) = 0.25 * ( PSI(I-1,J) + PSI(I+1,J) + PSI(I,J-1) + PSI(I,J+1) );
elseif (ID(I,J) == 1)
%AT POINTS NEIGHBORING RECTANGULAR BOUNDARY:
PSI(I,J) = ( PSI(I-1,J)/(A(I,J)*(A(I,J)+B(I,J))) + PSI(I+1,J)/(B(I,J)*(A(I,J)+B(I,J))) + PSI(I,J-1)/(C(I,J)*(C(I,J)+D(I,J))) + PSI(I,J+1)/(D(I,J)*(C(I,J)+D(I,J))) ) / ( 1./(A(I,J)*B(I,J)) + 1./(C(I,J)*D(I,J)) );
elseif (ID(I,J) == 2)
%AT GRID POINTS AWAY FROM THE CURVED BOUNDARY:
PSI(I,J) = PSI(2,J);
end
ERROR = ERROR + abs(PSI(I,J)-PSIPRV);
end
end
end
% ***** PLOT REPRESENTATIVE STREAMLINES *****
% [Q,P] = contour(PSI');
% clabel(Q);
% title(['Chanel flow past a circular cylinder: number of iterations = ' num2str(ITER)])
% xlabel('X - AXIS')
% ylabel('Y - AXIS')
% figure;
% for I=1:M
% for J=N:-1:1
% PSIF(I,N-J+1)=PSI(I,J);
% end
% end
% [QF,PF] = contour3(PSIF');
% clabel(QF);
% title('Stream function plot.')
% xlabel('X - AXIS')
% ylabel('Y - AXIS')
X(1) = 0.0;
for I=2:M
X(I) = X(I-1) + H;
end
Y(1) = 0.0;
for J=2:N
Y(J) = Y(J-1) + H;
end
hold on;
for K=1:length(CONST)
[XX,YY,KMAX] = SEARCH4(X,Y,PSI,CONST(K));
plot(XX(1:KMAX),YY(1:KMAX),SYMBOL(K))
end
hold off;
title(['Chanel flow past a rectangular cylinder: number of iterations = ' num2str(ITER)])
xlabel('X - AXIS')
ylabel('Y - AXIS')
legend(['PSI=' num2str(CONST(1))],['PSI=' num2str(CONST(2))],['PSI=' num2str(CONST(3))],['PSI=' num2str(CONST(4))],['PSI=' num2str(CONST(5))],['PSI=' num2str(CONST(6))],['PSI=' num2str(CONST(7))],['PSI=' num2str(CONST(8))],['PSI=' num2str(CONST(9))],['PSI=' num2str(CONST(10))],['PSI=' num2str(CONST(11))])
hold on;
for K=1:length(CONST)
[XX,YY,KMAX] = SEARCH4(X,Y,PSI,CONST(K));
plot(XX(1:KMAX),-YY(1:KMAX),SYMBOL(K))
end
hold off;
0 Comments
Answers (0)
See Also
Categories
Find more on Map Display 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!