Stokes stream function for spherical system

Can anybody see where I am going wrong? Trying to produce a contour plot of the Stokes stream function with radial and polar velocity of u1 and u2.

Number_of_grid_points = 500;
scale = 200;

xx0 = -1; %minimum
xx = 1; %maximum
dx = (xx-xx0)/Number_of_grid_points; %step

a = 1; %radius 
W = 1; %uniform velocity W along z axis

The = xx0:dx:xx;
R = The;
[the, r] = meshgrid(The,R);

u1 = (W * cos(the)).*(1 + (a^3)/(2*r^3) - (3*a)/(2*r)); %radial velocity 
u2 = -(W * sin(the)).*(1 - (a^3)/(4*r^3) - (3*a)/(4*r)); %polar velocity 
Z1 = u1 + u2 %velocity of fluid
contour(Z1,scale);

Answer 1

[AHMAD TALAEI]

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

[Harshit Anand]

Hi, can you tell me how to plot the streamlines vertically? In my problem statement ball is falling vertically. So the streamline and vector field shoud be in upward direction.

[AHMAD TALAEI]

Just contourplot z, which is the stream function. Also, notice that stream function is a scalar quantity, so it has no vertical component:

contour(x,y,real(z),'LineStyle','--','LineWidth',1.5,'LineColor',[0.5 0.5 0.5],'LevelStep',0.4);