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

2 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
Chukwuemeka Edward
Chukwuemeka Edward on 2 Oct 2022
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
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.

Categories

Find more on 2-D and 3-D Plots in Help Center and File Exchange

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!