Rankine Leading Edge Graph

m = 1;
V = 1;
disp('V     m')
disp([V m])
disp('Velocity Potential:')
disp('Stream function:')
disp('psi=V*y+(m/2/pi * atan2(y,x)')
disp('The (x,y) components of velocity (u,v):')
disp('u = V + m/2pi*x/(x^2+y^2)')
disp('v = m/2pi * y/(x^2+y^2)')
xstg = -m/2/pi/V;
ystg = 0; %stagnation point locations
N = 1000;
xinf = 3;
xd = xstg:xinf/N:xinf;
for n = 1:length(xd)
    if n ==1
        yd(1) = 0;
    else
        yd(n) = m/2/V;
        for it = 1:2000
            yd(n) = (m/2/V)*(1-1/pi*atan2(yd(n),xd(n)));
        end
    end
end
xL(1) = xd(end);
yL = -yd(end);
for nn = 2:length(xd)-1
    xL(nn) = xd(end-nn);
    yL(nn) = -yd(end-nn);
end
plot([xd xL],[yd yL],'k',[-1 3],[0 0],'k'), axis([-1 3 -1 1])
u = V +(m/2*pi)*(xd./(xd.^2+yd.^2));
v = m/(2*pi)*(yd./(xd.^2+yd.^2));
Cp = 1 - (u.^2+v.^2)/V^2;
hold on 
plot(xd,Cp); axis([-1 3 -1 1])
plot(0,m/V/4,'o')
plot(xstg,ystg,'o')
plot([1 3],[m/2/V],'--k')
[Cpmin ixd] = min(Cp);
xmin = xd(ixd);
ymin = yd(ixd);
plot(xmin,ymin,'+r')
Cpmin
phi = V*xd + (m/4/pi).*log(xd.^2+yd.^2);
dx = diff(xd);
dy = diff (yd);
ds = sqrt(dx.^2+dy);
dph = diff(phi);
ut = dph./ds;
xm = xd(1:end-1)+dx/2;
psi= V*yd + (m/2/pi).*atan2(yd,xd);
plot(xm,1-ut.^2/V^2,'r')
th = 0:pi/25:2*pi;
r = (m/2/pi/V)*(pi-th)./sin(th);
xb=r.*cos(th);
yb=r.*sin(th);
plot(xb,yb,'om')
thm = 0;
for nit = 1:1000
    thm = atan2(pi-thm,pi-thm-1);
end
thdegrees = thm*180/pi;
rm=(m/2/pi/V)*(pi-thm)/sin(thm)
xm = rm*cos(thm);
ym = rm*sin(thm);
plot(xm,rm,'dk')
um = V +m/2/pi*xmin/(xmin^2+ymin^2);
vm = m/2/pi*ymin/(xmin^2+ymin^2);
Cpm = 1 - (um^2+vm^2)/V^2
title('Rankine Nose Outline')

Answer 1

[William Rose]

@Roxanne Willisamson,

What is the question?

This line is missing something:

Cp = 1 - (u.^2+v.^2)V^2;

We must add * or / or + or - or something. I choose to add "/" to keep the dimensions consistent:

Cp = 1 - (u.^2+v.^2)/V^2;

After that, I get an error at the for loop:

zd = xstg:xinf/N:xinf;

for n = 1:length(xd)...

which is fixed if I change zd to xd:

xd = xstg:xinf/N:xinf;

for n = 1:length(xd)...

Then I get a "dot indexing is not supported" error at this line

yd(n) = (m/2/V)*(1-1/pi*atan2(yd(n).xd(n)));

so I replace the . with ,

yd(n) = (m/2/V)*(1-1/pi*atan2(yd(n),xd(n)));

Then I get a lot of numbers on the console, so I add a semicolon to the end of a line

yL(nn) = -yd(end-nn)

Then I get another "dot indexing is not supported" error at line

plot(xd.Cp).axis([-1 3 -1 1])

so I change it to

plot(xd,Cp); axis([-1 3 -1 1])

The I get error "Unable to resolve the name xstf.ystg." at line

plot(xstf.ystg,'o')

so I change it to

plot(xstf,ystg,'o')

and now the error at same line is "Unrecognized function or variable 'xstf'.", so I change it again to

plot(xstg,ystg,'o')

Then it is "Unrecognized function or variable 'dx'." at line

xmin = dx(ixd);

so I change it to

xmin = xd(ixd);

Then I get error "Unrecognized function or variable 'dy'." at line

dy - diff (yd);

so I change it to

dy = diff (yd);

Then an error at line

xm = xd(1:end01)+dx/2;

so I change it to

xm = xd(1:end-1)+dx/2;

The dot indexing problem at line

psi= V*yd + (m/2/pi).*atan2(yd.xd);

so I change it to

psi= V*yd + (m/2/pi).*atan2(yd,xd);

Then "Unrecognized function or variable 'ym'." at line

plot(xm,ym,'dk')

so I change it so

plot(xm,rm,'dk')

Then the script runs to completion without error. It produces a plot and displays some results on the console. See below.

%roxanneWilliamsonCp.m
%By Roxanne Williamson with help from W. Rose
m = 1;
V = 1;
disp('V     m')

V     m

disp([V m])

     1     1

disp('Velocity Potential:')

Velocity Potential:

disp('phi = V*x+(m/4/pi)*log(x^2+y^2)')

phi = V*x+(m/4/pi)*log(x^2+y^2)

disp('Stream function')

Stream function

disp('psi=V*y+(m/2/pi * atan2(y,x)')

psi=V*y+(m/2/pi * atan2(y,x)

disp('The (x,y) components of velocity (u,v):')

The (x,y) components of velocity (u,v):

disp('u = V + m/2/pi*xc/(x^2+y^2)')

u = V + m/2/pi*xc/(x^2+y^2)

disp('v = m/2/pi * y/(x^2+y^2)')

v = m/2/pi * y/(x^2+y^2)

xstg = -m/2/pi/V;
ystg = 0; %stagnation point locations
N = 1000;
xinf = 3;
xd = xstg:xinf/N:xinf;
for n = 1:length(xd)
    if n ==1
        yd(1) = 0;
    else
        yd(n) = m/2/V;
        for it = 1:2000
            yd(n) = (m/2/V)*(1-1/pi*atan2(yd(n),xd(n)));
        end
    end
end
xL(1) = xd(end);
yL = -yd(end);
for nn = 2:length(xd)-1
    xL(nn) = xd(end-nn);
    yL(nn) = -yd(end-nn);
end
plot([xd xL],[yd yL],'k',[-1 3],[0 0],'k'), axis([-1 3 -1 1])
u = V +m/2/pi*yd./(xd.^2+yd.^2);
v = m/2/pi*yd./(xd.^2)/V^2;
Cp = 1 - (u.^2+v.^2)/V^2;
hold on 
plot(xd,Cp); axis([-1 3 -1 1])
plot(0,m/V/4,'o')
plot(xstg,ystg,'o')
plot([1 3],[m/2/V],'--k')
[Cpmin ixd] = min(Cp);
xmin = xd(ixd);
ymin = yd(ixd);
plot(xmin,ymin,'+r')
Cpmin

Cpmin = -2.7447e+12

phi = V*xd + (m/4/pi).*log(xd.^2+yd.^2);
dx = diff(xd);
dy = diff (yd);
ds = sqrt(dx.^2+dy);
dph = diff(phi);
ut = dph./ds;
xm = xd(1:end-1)+dx/2;
psi= V*yd + (m/2/pi).*atan2(yd,xd);
plot(xm,1-ut.^2/V^2,'r')
th = 0:pi/25:2*pi;
r = (m/2/pi/V)*(pi-th)./sin(th);
xb=r.*cos(th);
yb=r.*sin(th);
plot(xb,yb,'om')
thm = 0;
for nit = 1:1000
    thm = atan2(pi-thm,pi-thm-1);
end
thdegrees = thm*180/pi;
rm=(m/2/pi/V)*(pi-thm)/sin(thm)

rm = 0.3650

xm = rm*cos(thm);
plot(xm,rm,'dk')

um=V+m/2/pi*xmin/(xmin^2+ymin^2);
vm = m/2/pi*ymin/(xmin^2+ymin^2);
Cpm = 1-(um^2+vm^2)/V^2

Cpm = -0.4048

One adresses errors one at a time until there are none left. Good luck.

[Roxanne Williamson]

Thats much better thank you! Do you perhaps know how to incorperate the quiver and streamline function to get flow arrows?

[William Rose]

@Roxanne Williamson

I am familiar with quiver and streamline plots and have used both. You need to have four 2-D arrays of data, all of same size, for quiver(). You have four 1-D arrays: u, v, xd, yd. The common practice is to construct arrays X and Y with meshgrid(). The values of u and v must be calculated at each point in the 2D grid represented by X and Y.

You could make a grid with

[X,Y]=meshgrid(xd,yd);

but X and Y would each have over a million points, which is not useful or practical for a quiver plot. You could downsample vectors xd and yd to make a mesh grid at 50x50 points. Then you would compute u an v at the 2500 mesh points. Thus u and v would also be 50x50. Your code at present only evaluates u and v along a 1-dimensional curve.

streamline() requires the same 2D input arrays (X,Y,U,V) as quiver(). streamline() also requires vectors (startx, starty) of 1 or more starting points for the streamlines. The example in the help for streamline() is good. It shows how to use quiver() and streamline() together to make a nice plot.

[William Rose]

@Roxanne Williamson,

Your code includes

disp('The (x,y) components of velocity (u,v):')
disp('u = V + m/2/pi*xc/(x^2+y^2)')   %Should xc be x?
disp('v = m/2/pi * y/(x^2+y^2)')

and later you have

u = V +m/2/pi*yd./(xd.^2+yd.^2);    %Should yd be xd?
v = m/2/pi*yd./(xd.^2)/V^2;         %Looks wrong.

which have multiple inconsistencies.

Also, for the disp() equations to be correct, m must have dimensions of area/time. Is that your understanding of m? It is a good practice to add a comment with the units, after each constant you define. This helps you understand the problem, check your work, and find errors.

[William Rose]

@Roxanne Williamson,

Your equations are starting to make sense as I think about them. You have a planar system. A fluid emanates from point 0,0 and spreads out radially. As it spreads, it slows down, as you would expect if the 2D fluid is incompressible. There is an external "wind" blowing across the system with velocity V, in the +x direction. Constant m is the "volume flow rate" out of the source, except its units are area/time, rather than volume/time, because the system is 2-dimensional, not 3D. The velocity at the source point will be undefined (divide by zero error).

[William Rose]

@Roxanne Williamson,

Here is a script that plots velocity vectors and streamlines. See comments in the code for explanation. I made some changes to the calculations, where I thought there were errors in the original, as I described in my previous post. The difference between the two plots is the region covered and the locations of the streamline starting points.

[Roxanne Williamson]

That is awesome! Thank you so much. Perfect visualisation!

[William Rose]

@Roxanne Williamson, you're welcome.