Index exceeds the number of array elements. Index must not exceed 20. Error in lumped_vor​tex_parabo​lic_camber (line 54) dx = xcollocation(i) - xvortex(k);

After runnung this code on lumped vortex model it gave an error on line 52(Bold section) as follows.

%### Index exceeds the number of array elements. Index must not exceed 20.

Error in lumped_vortex_parabolic_camber (line 54)

dx = xcollocation(i) - xvortex(k);###%

Someone please help.

% Lumped Vortex Panel Code

% This code is for a thin cambered plate

clear all;

%close all;

% Input Data

chord = 1;

Uinfinity = 1;

alpha = 3 * pi / 180;

Npanel = 20;

% Location of Lumped Vortex, correlation point, panel start and end points

% Assuming constant panel length - which can be changed

xpanel = chord / Npanel;

for i = 1 : 1 : Npanel

xstart(i) = (i - 1) * xpanel;

zstart(i) = 0.4 * xstart(i) * (1 - xstart(i));

zstart(i) = -0.4 * xstart(i) * (1 - xstart(i));

xend(i) = (i) * xpanel;

zend(i) = 0.4 * xend(i) * (1 - xend(i));

zend(i) = -0.4 * xend(i) * (1 - xend(i));

xvortex(i) = 1 / 4 * xpanel + (i - 1) * xpanel;

zvortex(i) = 0.4 * xvortex(i) * (1 - xvortex(i));

zvortex(i) = -0.4 * xvortex(i) * (1 - xvortex(i));

xcollocation(i) = 3 / 4 * xpanel + (i - 1) * xpanel;

zcollocation(i) = 0.4 * xcollocation(i) * (1 - xcollocation(i));

zcollocation(i) = -0.4 * xcollocation(i) * (1 - xcollocation(i));

end

% Find panel angles

for i = 1 : 1 : Npanel

dx = xend(i) - xstart(i);

dz = zend(i) - zstart(i);

th(i) = atan2(dz,dx);

end

th = [th,-th];

% Find normals and tangents

nx = -sin(th);

nz = cos(th);

tx = nz;

tz = -nx;

nx = sin(th);

nz = -cos(th);

tx = nz;

tz = -nx;

% Build the influence coefficients

for i = 1 : 1 : 2*Npanel % Loop over the collocation points

for k = 1 : 1 :2*Npanel % Loop over the influence from vortex

dx = xcollocation(i) - xvortex(k);

dz = zcollocation(i) - zvortex(k);

r2 = dx ^ 2 + dz ^2;

velocity_influence = 1 / 2 / pi / r2 * [0 1;-1 0] * [dx;dz];

% You need to calculate the normal of the velocity induced at the

% panel point from the singularity. This is done by taking the dot

% product as shown below

a(i,k) = velocity_influence' * [nx(i); nz(i)];

end

% The velocity from the free stream. The do product is to calculate the

% normal

b(i,1) = - Uinfinity *[cos(alpha) sin(alpha)] * [nx(i); nz(i)];

end

% Determine singularity strength

Gamma = inv(a) * b;

% Plotting results

figure; plot(xvortex, 2*Gamma / xpanel, 'Linewidth', 2)

% Analytic Solution

dcp = 4 * sqrt((chord - xcollocation)./(xcollocation)) * alpha + ...

32 * 0.1 * sqrt(xcollocation.*(1-xcollocation));

hold; plot(xvortex, dcp, 'or')

xlabel('x/c'); ylabel('\Delta c_p')

legend('Numerical','Analytical')

% Lumped Vortex Code Finished ... Fancy Processing Below

% If you want to build the velocity vectors and streamlines, just build a

% grid and calculate the velocity at grid points

xmin = -2;

xmax = 3;

zmin = -1;

zmax = 1;

Nx = 20;

Nz = 20;

dx = (xmax - xmin)/Nx;

dz = (zmax - zmin)/Nz;

% Starting points for the streamlines

startx1 = xmin : dx : xmax;

startz1 = ones(size(startx1)) * zmin;

startz2 = zmin : dz : zmax;

startx2 = ones(size(startz2)) * xmin;

[x,z] = meshgrid( xmin : dx : xmax , zmin : dz : zmax);

u = zeros(Nz + 1, Nx + 1);

w = zeros(Nz + 1, Nx + 1);

% Calculate the velocity at the grid points

for k = 1 : 1 : Nz + 1

for i = 1 : 1 : Nx + 1

for j = 1 : 1 : Npanel

dx = x(k, i) - xvortex(j);

dz = z(k, i) - zvortex(j);

r2 = dx ^ 2 + dz ^2;

velocity_influence = 1 / 2 / pi / r2 * Gamma(j) * ...

[0 1;-1 0] * [dx;dz];

u(k,i) = u(k,i) + velocity_influence(1);

w(k,i) = w(k,i) + velocity_influence(2);

end

u(k,i) = u(k,i) + Uinfinity * cos(alpha);

w(k,i) = w(k,i) + Uinfinity * sin(alpha);

end

end

figure; hold on;

streamline(x,z,u,w,startx1,startz1);

streamline(x,z,u,w,startx2,startz2)

% Plot the airfoil surface

for j = 1 : 1 : Npanel

xtemp = [xstart(j) xend(j)];

ztemp = [zstart(j) zend(j)];

plot(xtemp, ztemp,'k','Linewidth',2)

end

axis image; box on

title('Streamlines'); xlabel('x/c'); ylabel('z/c')

Answer 1

[VBBV]

for i = 1 : 1 : Npanel % change

In the first case you used this. But look at the 2nd for loop why you changed to 2*Npanel.

[VBBV]

For influence coefficients loop, the number of panels are higher 2 times what they are computed earlier. So,size of xcollocation and xvortex are same as Npanel. When you try to access the elements greater than its size you get such error

[Chukwuemeka Edward]

I am writing the code for two panels(inverted parabola) with exact same x-axis , hence the 2* Npanel. So there is something wrong with line 54, i am not sure how to fix.

[VBBV]

% Lumped Vortex Panel Code
% This code is for a thin cambered plate
clear all;
%close all;
% Input Data
chord = 1;
Uinfinity = 1;
alpha = 3 * pi / 180;
Npanel = 20;
% Location of Lumped Vortex, correlation point, panel start and end points
% Assuming constant panel length - which can be changed
xpanel = chord / Npanel;
for i = 1 : 1 : 2*Npanel
  xstart(i) = (i - 1) * xpanel;
  zstart(i) = 0.4 * xstart(i) * (1 - xstart(i));
  zstart(i) = -0.4 * xstart(i) * (1 - xstart(i));
  xend(i) = (i) * xpanel;
  zend(i) = 0.4 * xend(i) * (1 - xend(i));
  zend(i) = -0.4 * xend(i) * (1 - xend(i));
  xvortex(i) = 1 / 4 * xpanel + (i - 1) * xpanel;
  zvortex(i) = 0.4 * xvortex(i) * (1 - xvortex(i));
  zvortex(i) = -0.4 * xvortex(i) * (1 - xvortex(i));
  xcollocation(i) = 3 / 4 * xpanel + (i - 1) * xpanel;
  zcollocation(i) = 0.4 * xcollocation(i) * (1 - xcollocation(i));
  zcollocation(i) = -0.4 * xcollocation(i) * (1 - xcollocation(i));   
end
% Find panel angles
for i = 1 : 1 : 2*Npanel
  dx = xend(i) - xstart(i);
  dz = zend(i) - zstart(i);
  th(i) = atan2(dz,dx);
end

Then you need to change the Npanel for above two loops to 2*Npanel to make same size vectors