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);

8 views (last 30 days)
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')

Answers (1)

VBBV
VBBV on 1 Oct 2022
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.
  3 Comments
VBBV
VBBV on 6 Oct 2022
Edited: VBBV on 6 Oct 2022
% 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

Sign in to comment.

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!