Question

Develop a MATLAB (or programming language of your choice) code to show the flow field for lifting cylinder flow using the following specifications.

Set your free-stream velocity, cylinder radius, and vortex strength in the code. There are a total of four cases for which you will run your code (or have all four in one script). You will only be changing the vortex strength for the different cases. Therefore, set ?∞ = 1 and ? = 1 for all cases. Set your vortex strength to the following values for each case, Case 1: ? = 0, Case 2: ? = 8, Case 3: ? = 16, Case 4: Find the value for vortex strength that produces one stagnation point on the surface of the cylinder.

Set V = R = 1; Set Γ = value desired for specific case; Use mesh grid to create my domain, but do this in polar coordinates; Calculate x and y coordinates based on mesh grid; Calculate stream function as a function of r and theta; Calculate Vr and Vtheta as a function of r and theta; Convert Vr and Vtheta to u and v; Plot contours of stream function; Plot velocity field using quiver function;

Answer 1

U_i = 20;            % Ambient velocity
a = 4;               % cylinder radius
c = -a*5;            % starting coordinate (x)
b = a*10;            % ending coordinate (x)
d = -60;             % starting coordinate (y)
e = 60;              % ending coordinate (y)
n = a*50;            % number of intervals (step size in grid)

[x,y] = meshgrid([c:(b-c)/n:b],[d:(e-d)/n:e]');

for i = 1:length(x)
    for k = 1:length(x);
        f = sqrt(x(i,k).^2 + y(i,k).^2);
        if f < a
            x(i,k) = 0;
            y(i,k) = 0;
        end
    end
end

% Definition of polar variables
r = sqrt(x.^2+y.^2);
theta = atan2(y,x);    
               
%% Creation of Streamline function
z = U_i.*r.*(1-a^2./r.^2).*sin(theta);%- G*log(r)/(2*pi);
%% Creation of Figure
m = 100;
s = ones(1,m+1)*a;
t = [0:2*pi/m:2*pi];

%% Streamline plot
contour(x,y,z,50);
hold on 
polar(t,s,'-k');

% axis square
title('Stream Lines');
grid off