How to convert plot a contour plot or streamline plot with this code?
6 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
% ***** 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 comentarios
Respuestas (0)
Ver también
Categorías
Más información sobre Map Display en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!