Fourier Series for thin airfoil theory

I'm trying to use a Fourier Series Approximation, specifically that specified in thin airfoil theory, to write a function that will do this more universally than hard coding it. So far what I have is the following code. It runs, and my value for is correct, however none of my other values are even remotely close to correct. I believe the problem is coming from the way that I find dydx, but I don't understand why the same method works for and no the other A values.

clear all;close all;clc
alph = 0;
[ ~, ~, ~, xcam, ~ ]= NACA4Digit('2412',100,100,100);
[ ~, ~, ~, ~, ycam ]= NACA4Digit('2412',100,100,100);
ncoeff = 15;
% %% Error Messages
% if size(xcam) ~= size(ycam)
%     error('Camber line vector sizes do not match.');
%     return;
% else
% end

%% Variable Initialization
iter = length(xcam);
Uinf = 1; % Freestream Velocity, value is a placeholder as the final 
% results are non-dimensionalized. 
c = iter; % Chord Length, value is a placeholder as the final results are non-
% dimensionalized.
alpha = deg2rad(alph); % Converts the angle of attack to rad


%% Constant Functions
DF = gradient([xcam;ycam]); % Takes the gradient at the given points
[~,fore] = max(ycam); % Finds the point within the array where the camber is greatest, important for differentiating
zeroToFore = linspace(0,fore,iter); % Defines Integration Bounds
zeroToFore = acos(1 - ((2*zeroToFore)/c)); % Converts to theta from x
foreToAft = linspace(fore,pi,iter); % Defines Integration Bounds
foreToAft = acos(1 - ((2*foreToAft)/c)); % Converts to theta from x
A = zeros(1,ncoeff-1); % Intializes A 
A0 = alpha - (1/pi)*(trapz(zeroToFore,DF(1,:)) + trapz(foreToAft,DF(2,:))); % Defines A_0
U_theta = zeros(1,iter); % Intitalizes U_theta

%% Solve

for i = 1:iter
theta = acos(1-((2*xcam(i))/c)); % Defines Theta

for n = 1:1:ncoeff

trapzDataFore =  DF(1,:)*cos(n*theta) % Trapz Data from the given function for thin airfoil theory using Fourier Series Approx, fore
trapzDataAft = DF(1,:)*cos(n*theta) % Trapz Data from the given function for thin airfoil theory using Fouerier Series Approx, aft
An(n) = (2/pi) * (trapz(zeroToFore,trapzDataFore) + trapz(foreToAft,trapzDataAft)) % A_n value for the Fourier Series, there are ncoeff A values



end
Eqn1 = 0 % Intitializes Equation 1
for n = 1:ncoeff
Eqn1 = Eqn1 + An(n)*sin(n*theta) % Sums all of the A_n Values
end
gamma = 2*Uinf*(A0*((1 + cos(theta))/(sin(theta))) + Eqn1) % Eqn for gamma at a given theta

U_theta(i) =  gamma/2; % U at theta
end


Cp(1) = 0;
Cp(iter) = 0;
C_p = 1 - ((Uinf+U_theta)/Uinf).^2;
Cl = pi*(2*A0 + A(1));
Cm = (pi/4)*(A(2)-A(1));

[KSSV]

Read about gradient function. You are taking only one ouput from graidient. There might be another output which is along y-axes.

[Tyler Reohr]

KSSV,

When I utilize the other output from gradient, I then have 3 DF outputs:

2x100 DFX, which is equal to what I had with just the DF output, and

2x100 DFY, which is 2 equal 1x100 rows

I'm assuming that you're right and I do need to incorporate this data, but I'm not exactly sure how do to so with the way I have the code set up, how would I define the dydx for different iterations in the first nested for loop? Under trapzDataFore and trapzDataAft?

For reference, if you're not familiar with thin airfoil theory, the equation I am going from is derived from:

and can be simplified to:

and trapzDataFore and trapzDataAft are the respective equations within the integral.

Thanks!

Tyler

Answer 1

[Paul]

From the equation it looks to me that y is scalar function and the independent variable x is also a scalar.

Are xcam and ycam both 1-dimensional?

In which case I think this line needs to change:

%DF = gradient([xcam;ycam]); % Takes the gradient at the given points
DF = gradient(ycam,xcam);

I've used that form of the gradient command for a long time, but it doesn't seem to be documented anymore.

[Tyler Reohr]

Paul,

They're an array of x-y coordinates. So each one is technically a one dimensional array, but together they represent a point on the 2-D airfoil, which is why I combined them. I tried switching to this and my value for , which was correct, was no longer correct. To be fully transparent, I'm not exactly sure how to get from the gradient function in the way that's useful for my code. For the term, it's relatively simple and I essentially just plugged in different ways of referencing DF until I got the correct answer. For the , I can't seem to get this to work.

In this equation, the only error I can find is in , and I believe it's because of how I use gradient, I just can't seem to do it properly and, reference it properly after the fact.

[Paul]

If xcam and ycam are an array of x-y coordinates, then I don't understand what dy/dx means, or what the code is actually trying to do. Maybe the notation, including the bracket/paren pair means something unique to this application? Without more information about the math, I'm afraid I can't help with the code.

[Tyler Reohr]

The call to the function NACA4Digit in the first 2 lines returns several different things. All of them besides xcam and ycam aren't really necessary for this (they're signficiant using the VPM to get the same coefficients, but for this we are told specifically to use thin airfoil theory). NACA4Digit takes in the Airfoil Designation as a string. The other 3 values are the number of points being considered at the top, bottom, and camber line of the the airfoil. I chose 100 arbitrarily based on it being accurate enough to gage if my answers are correct (based on the other method I mentioned, VPM), but still runs in less than a second. It then returns, based on the NACA 4 Digit Airfoil Designation, 2 separate arrays. The first, xcam, is a 1x100 (because of what I input to NACA4Digit.m) array of 100 equally spaced points along the camber line of the airfoil. The second, ycam, is a 1x100 array of where the corresponding y-value is for each of the variables in xcam.

is the change of y with respect to x. To fully derive where it comes from, I'd have to get back to you when I get to my textbook for the class, however the important equations, are as follows:

Finding theta:

where x is the x coordinate (from xcam) and c is the chord length (in this case, because the answer is eventually non-dimensionalized, it is arbitrary, so I made the chord length equal to the length of xcam to make for easier calculations within the loop, so it is currently set to 100 units).

Then the vortex strength, γ, found using a Fourier Series, is:

where,

and

These equations are the real meat behind any calculations, after this is simply a matter of plugging them in to much simplier equations to find the Coefficients I'm looking for:

The velocity, Uat θ is:

and

Which are the final values that I'm looking for.

EDIT:

I missed your inquiry about the brackets for xcam and ycam. Assuming that you mean where I define DF = gradient([xcam;ycam]), the purpose of this was an attempt to create an array that actual stored the entire point, instead of just the x-coordinates and y-coordinates. It is definitely possible that I did not do that the right way though. The only thing I have to check my values with, i.e. the only value I have gotten a correct answer for, is , and even that could be merely a coincidence.

[Paul]

So you have a function y = f(x). The elements of xcam are values of x that go into the function and the elements of ycam are the corresponding values of y at each value of xcam. If this is the case, then the gradient command that I showed above estimates dy/dx at each value of xcam. And the values of xcam don't even have to be equally spaced. For example:

x = 0:.01:2*pi;
x = sort(x + rand(size(x))); % make x not equally spaced
y = sin(x)

y = 1×629

0.0943 0.0986 0.1111 0.1959 0.2322 0.2478 0.2514 0.3152 0.3212 0.3214 0.3350 0.3470 0.3674 0.3789 0.3913 0.4398 0.4581 0.4632 0.5189 0.5199 0.5257 0.5551 0.5567 0.5703 0.5761 0.6392 0.6501 0.6721 0.6844 0.7018

figure;
plot(x,y)

figure
plot(x,cos(x),x,gradient(y,x))

So the gradient command estimates the derivative. I'd suggest do somehting similar to the above, i.e.,

plot(xcam,ycam)
plot(xcam,gradient(ycam,xcam))

and do some spot checks and see if the output of gradient() makes sense for your data.

[Tyler Reohr]

Paul,

I did as you suggested and can see that you're definitely right about the gradient. I got the following. I'm not exactly sure why if I'm being honest, but for some reason the way that I currently have the gradient setup yields the third graph, which appears to be a scaled version of the xcam and ycam points.

% subplot(3,1,1)
% plot(xcam,ycam)
% subplot(3,1,2)
% plot(xcam,gradient(ycam,xcam))
% subplot(3,1,2)
% plot(xcam,ycam,DF(1,:),DF(2,:))

[Tyler Reohr]

EDIT:

I misinterpreted the graph, it is in fact, much weirder than a scaled version of the xcam and ycam values: