gibson:teaching:spring-2016:math445:lecture:cylinderflow_code
- cylinderflow.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inviscid, irrotational flow around a cylinder
% Make quiver plot of velocity and contour plot of pressure for
% Vx(r,theta) = V0 (1 - a^2/R^2 cos(2 theta))
% Vy(r,theta) = -V0 a^2/r^2 sin(2 theta)
% p(r,theta) = 2 a^2/r^2 cos(2 theta) - a^4/r^4;
% for a=1, V0=1, on domain -4 <= x <= 4, -4 <= x <= 4
%%%%%%%%% Define constants %%%%%%%%%
V0 = 1;
a = 1;
%%%%%%%%% First draw the cylinder %%%%%%%%%
theta = linspace(0, 2*pi, 100);
plot(a*cos(theta), a*sin(theta), 'k-')
hold on
%%%%%%%%% Make the quiver plot of velocity %%%%%%%%%
x = linspace(-4,4,21); % quiver plots work better on coarse grids
y = linspace(-4,4,21);
[X,Y] = meshgrid(x,y);
R = sqrt(X.^2 + Y.^2);
Theta = atan2(Y,X);
% evaluate the formula for the velocity field
Vx = V0*(1 - a^2./R.^2 .* cos(2*Theta));
Vy = - V0*a^2./R.^2 .* sin(2*Theta);
% set Vx, Vy to zero inside the cylinder
Vx = (R > a) .* Vx;
Vy = (R > a) .* Vy;
% draw the quiver plot
quiver(x,y, Vx, Vy);
hold off
axis equal
axis tight
xlabel('x');
ylabel('y');
title('inviscid 2D flow around a cylinder')
gibson/teaching/spring-2016/math445/lecture/cylinderflow_code.txt · Last modified: 2016/04/14 09:38 by gibson