Use the following code. I found it online but did a few modifications to plot spherical stream function. For more info look at: https://www.mathworks.com/matlabcentral/fileexchange/4674-streamlines-magnus-and-cp-aroud-a-cylinder-section
close all
clear all
clc
V_0 = 1; % flow veloicty (m/s)
a = 1; % radius (m)
b = a*2;
c = -a*2;
n = a*15; % number of intervals
% the spherical coordinate solution for a falling sphere looking only at x-y plane (looking from side not top)
[x,y]=meshgrid((c:(b-c)/n:b),(c:(b-c)/n:b)');
% Preliminar DATA & purification
for i = 1:length(x)
for k = 1:length(x)
if sqrt(x(i,k).^2 + y(i,k).^2) < 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 the streamline function
z = V_0.*sin(theta).*r.*(1-(a.^2./(r.^2))); % laminar stream function in polar coordinate
% stream line Plot
figure(1), contour(x,y,real(z),25,'Linewidth',2); hold on; axis square;
% circle plot
r = ones(1,n+1)*a;
t = (0:2*pi/n:2*pi);
polar(t,r,'-k');
% creation of vectors around the circle
figure(2), contour(x,y,abs(z),10); axis square;
% recreation of x,y
n = a*15; % number of intervals
[x,y] = meshgrid((c:(b-c)/n:b),(c:(b-c)/n:b)');
for i = 1:length(x)
for k = 1:length(x)
if sqrt(x(i,k).^2 + y(i,k).^2) < a
x(i,k) = 0;
y(i,k) = 0;
end
end
end
r = sqrt(x.^2 + y.^2);
theta = atan2(y,x);
u_r = V_0 .* (1 - (3/2).*(a./r) + (1/2).*(a./r).^3) .* cos(theta); % laminar stream function in spherical coordinate
u_theta = - V_0 .* (1 - (3/4).*(a./r) - (1/4).*(a./r).^3) .* sin(theta); % laminar stream function in spherical coordinate
% velocities in x-y cordinates
u = u_r .* cos(theta) - u_theta .* sin(theta);
v = u_r .* sin(theta) + u_theta .* cos(theta);
% vectors plotting
figure(2); hold on
quiver(x,y,u,v)
% creating The Filled Circle
t_r = (0:.1:2*pi);
xxx = a*cos(t_r);
yyy = a*sin(t_r);
fill(xxx,yyy,'y'); axis square
