How to plot Rankine Half Body Streamlines

Hello,

I am attempting to plot the streamlines and velocity potential lines for a superposition of a source + uniform velocity field (Rankine Half Body). However, my code is producing the following streamlines:

They should look like the following:

My code is the following:

syms x y

U = 2.5; %Free Stream Velocity
m = 10/pi; %Source Strength
d = -m/U; %Stagnation Point

psi = U*sqrt(x^2+y^2)*sin(atan(y/x))+m*atan(y/x);
phi = U*sqrt(x^2+y^2)*cos(atan(y/x))+m*log(sqrt(x^2+y^2));

for i = 5:15

    hold on
    a = ezplot(psi == i);
    b = ezplot(phi == i);

    set(a,'Color','black',"LineStyle","-");
    set(b,"LineStyle","-.");

    if i == m*pi %Rankine Half Body Streamline
        set(a,'Color','blue',"LineStyle","-");
    end
    
end

Answer 1

[Star Strider]

I can’t run that because the relevant data aren’t provided. Tha atan calls could be the problem. Substituting‘atan(y/x)’ with ‘atan2(y,x)’ could be helpful. See the documentation for atan2 for details.

[Peter Di Palma]

Apologies, x and y are sym variables. I've updated the question.

[Star Strider]

The atan2 substitution produced an improved plot, however since I don¹t understand what you’re doing, I can make no further suggestions. It could be worth experimenting with fcontour or fimplicit if you are looking at a bivariate function and want to assess its contours specifically at various levels or at zero, instead of the ezplot calls.

syms x y

U = 2.5; %Free Stream Velocity
m = 10/pi; %Source Strength
d = -m/U; %Stagnation Point

psi = U*sqrt(x^2+y^2)*sin(atan2(y,x))+m*atan2(y,x);
phi = U*sqrt(x^2+y^2)*cos(atan2(y,x))+m*log(sqrt(x^2+y^2));

for i = 5:15

    hold on
    a = ezplot(psi == i);
    b = ezplot(phi == i);

    set(a,'Color','black',"LineStyle","-");
    set(b,"LineStyle","-.");

    if i == m*pi %Rankine Half Body Streamline
        set(a,'Color','blue',"LineStyle","-");
    end
    
end

.

[Peter Di Palma]

Thank you. I was unaware that atan only operated on the interval of -pi/2 pi/2.

[Star Strider]

My pleasure!

I’m not certain where you want to go from here.

[Peter Di Palma]

That was all. I was only confused why the plot was not showing up correctly.

Any idea as to why the plot is only in the range of [-6,6] for both x and y? If I try to extend the axes past this nothing is plotted.

[Star Strider]

Using fcontour seems to get closer to the image you posted.

I’ve plotted ‘phi’, ‘psi’ and the sum and product of both of them here, as well as overlaid on the same plot, since I’m not certain how they’re supposed to interact.

This seems closer to what you want —

syms x y

U = 2.5; %Free Stream Velocity
m = 10/pi; %Source Strength
d = -m/U; %Stagnation Point

psi = U*sqrt(x^2+y^2)*sin(atan2(y,x))+m*atan2(y,x);
phi = U*sqrt(x^2+y^2)*cos(atan2(y,x))+m*log(sqrt(x^2+y^2));

figure
hfc = fcontour(phi, [-1 1  -1 1]*2*pi, 'MeshDensity',150);
hfc.LevelList = linspace(-25, 25, 50);
colormap(turbo)

figure
hfc = fcontour(psi, [-1 1  -1 1]*2*pi, 'MeshDensity',150);
hfc.LevelList = linspace(-25, 25, 50);
colormap(turbo)

figure
hfc = fcontour(phi+psi, [-1 1  -1 1]*2*pi, 'MeshDensity',50);
hfc.LevelList = linspace(-50, 50, 50);
colormap(turbo)

figure
hfc = fcontour(phi*psi, [-1 1  -1 1]*2*pi, 'MeshDensity',50);
hfc.LevelList = linspace(-500, 500, 50);
colormap(turbo)

figure
hfc = fcontour(psi, [-1 1  -1 1]*2*pi, 'MeshDensity',150);
hfc.LevelList = linspace(-25, 25, 50);
hold on
hfc = fcontour(phi, [-1 1  -1 1]*2*pi, 'MeshDensity',150);
hfc.LevelList = linspace(-25, 25, 50);
hold off

Make appropriate changes to get the result you want.

.