April 2020
Overview
For this project I created a calculator using MATLAB which uses the Vortex Lattice Method to find the lift and coefficient of lift for a given wing geometry.
Theory
The Vortex Lattice Method (VLM) is a simple aerodynamic model and an extension to the Prandtl Lifting Line, that is used in computational fluid dynamics. The VLM models the lifting surfaces, such as a wing, of an aircraft as an infinitely thin sheet of discrete vortices to compute lift and induced drag. The influence of the thickness and viscosity is neglected.
VLM is based around the Kutta Joukowski theorem that lift can be related to circulation (Gamma), which is essentially the swirl created in the flow when lift is created.
If we solve for this circulation on all of the panels that we created, then we will be able to find the lift created. Since we know that the lift must go to zero at the tips of the wings, we have a boundary condition where we can start our solve.
In summary, VLM splits the wing into chordwise sections and solves for circulation on each of these panels. This circulation distribution can then be used to find the Lift of the wing.
Solution
On the wing geometry the A and B points show where I solved for circulation, while the K points show the boundary conditions defined where velocity into the wing is zero.
The definition of the wing geometry is shown in the first lines of the code. For this wing I found a Lift of 57N and a Lift coefficient of 0.1548
Code
%% Inputs
syms y
N = 20; %Number of Panels on one half of wing
b = 6; %Wing Span (meter)
cr = 1; %Wing root chordlength (meter)
lambda = 1; %taper ratio lambda=ct/cr
sweep = 20*pi/180; (deg to rad)
vinf = 10; %Bulk wind velocity (m/s)
alfa = 2*pi/180; %angle of attack
yaw = 0; %yaw angle
rho = 1.225; %Density of Air
s=cr*(1+lambda)*b/2; %Planform Area
vmat=[vinf*cos(alfa),vinf*sin(yaw),vinf*sin(alfa)];
%% Define Equations and Mesh
A=[b]; %initialize with sym so you can add syms later
B=[b];
K=[b];
%Position equations
xle=abs(y)*tan(sweep); %Leading edge x position
c=cr+2/b*cr*(lambda-1)*abs(y); %chord length
xqrt=xle+c/4; %Quarter chord x position
x3qrt=xle+3*c/4; %Three quarter chord x posistion
for n=1:2*N %Get coordinates for all horseshoe filament vertices and BCs
A(n,2)=b*(n-1)/(2*N)-b/2;
A(n,1)=subs(xqrt,y,A(n,2));
B(n,2)=b*(n)/(2*N)-b/2;
B(n,1)=subs(xqrt,y,B(n,2));
K(n,2)=(A(n,2)+B(n,2))/2;
K(n,1)=subs(x3qrt,y,K(n,2));
A(n,3)=0;B(n,3)=0;K(n,3)=0;
end
%Plot filaments vertices and boundary points
figure(1)
scatter(A(:,1),A(:,2))
hold on
scatter(B(:,1),B(:,2)),scatter(K(:,1),K(:,2))
fplot(finverse(xle),sort([0,b/2*tan(sweep)]),'k')
fplot(-finverse(xle),sort([0,b/2*tan(sweep)]),'k')
fplot(finverse(xle+c),sort([cr,b/2*tan(sweep)+cr*lambda]),'k')
fplot(-finverse(xle+c),sort([cr,b/2*tan(sweep)+cr*lambda]),'k')
legend('A','B','K','Wing Edges')
%% Solve for Ciruclation Gamma, and Lift Force
%Create Matrix to Solve for Gamma
hmat=ones(N,N);
for i=1:N
for j=1:N
hmat(i,j)=dot(horseshoe(A(j,:),B(j,:),K(i,:))...
+horseshoe(A(2*N+1-j,:),B(2*N+1-j,:),K(i,:)),[0,0,1]);
end
end
%Solve For Gamma
G= hmat^-1*(dot(-vmat,[0,0,1]).*ones(N,1));
%Calculate Cl
CL=4*sum(G)/(N*vinf*cr*(lambda+1))
L=CL/2*rho*vinf^2*s
%Plot Gamma vs span(y)
ybodypoints=transpose(K(:,2));
gamgam= [transpose(G),fliplr(transpose(G))];
figure(2)
plot(ybodypoints,gamgam),title('Gamma vs span'),xlabel('Span'),ylabel('Gamma')
function V = horseshoe(A,B,K)
%This function takes in vertices of horseshoe filament A,B
% point of influence K, and circulation G
v1=BSvelocity([99999999,A(2),A(3)],A,K);
v2=BSvelocity(A,B,K);
v3=BSvelocity(B,[9999999,B(2),B(3)],K);
V=v1+v2+v3;