Troubleshooting Airfoil Calculation Using Joukowski Transformation

[danja]

I'm having a bit of trouble with a homework assignment for a fluid mechanics course. I'd like to ask if my solution method is appropriate. The goal of the assignment is to determine the pressure coefficient distribution along the upper and lower airfoil surfaces and plot it (I'm using MATLAB for this). To do so, we're supposed to use the joukowski transformation.

Allow me to define quantities of interest:

a = circle radius
L = airfoil chord length
c = L/4
h = camber
t = airfoil thickness
m*exp(i*delta) = -0.77*t*c/L + i*h/2 (parameter in the complex potential function).

Transforming from Zeta plane (circle) to Z plane (normal shape).

In order to obtain the pressure coefficient, I am using the following procedure:

1) Determine values for a, c, h, m*exp(i*delta) for the desired conditions
2) Calculate X and Y locations along upper and lower surfaces
3) Transform X and Y data into the Zeta plane using the quadratic solution to Z(zeta)
4) Take the derivative of F(zeta) and use the data from steps 1 and 3 to
determine W(zeta) for all X and Y
5) Transform W(zeta) into W(z) by applying the operation: W(z) =
W(zeta)/dZ/dZeta
6) Apply Bernoulli's equation while using the real parts of W(z) as the
local velocity to find P(x)
7) Calculate Cp(x) based on P(x).

Results appear to be correct (in trend anyways) for the case when camber is
zero. But when I add camber, I am getting a discontinuity along the lower
surface. I have already set the appropriate roots of the quadratic equation
used in step 3 so the plots are using the right pieces. I can post screen shots if it helps.

I'd really appreciate a little help on this. Is this an appropriate method?
Or is it better to try and do everything in the Z-plane instead?

I've already been through this thread:
https://www.physicsforums.com/showthread.php?t=365978&highlight=joukowski+pressure

But it doesn't help me too much.

Matlab code:

close all; clear; clc;

camb_per = 0.0002;
t_per = 0.12;
at = 3;

L = .5;
U = 50;
rho = 1.2;
Pinf = 101325;

h = camb_per*L;
t = t_per*L;
c = L/4;
M = 1i*h/2 - 0.77*t*c/L;
a = c + 0.77*t*c/L;
at = at*pi/180;

% Circulation
Gam = pi*U*L*(1 + 0.77*t/L)*sin(at + 2*h/L);
K = Gam/(2*pi);

airshapeU = @(x) sqrt((L^2/4)*(1 + L^2/(16*h^2) ) - x^2) - L^2/(8*h) + 0.385*t*(1 - 2*x/L)*sqrt(1 - (2*x/L)^2);
airshapeL = @(x) sqrt((L^2/4)*(1 + L^2/(16*h^2) ) - x^2) - L^2/(8*h) - 0.385*t*(1 - 2*x/L)*sqrt(1 - (2*x/L)^2);

%fplot(airshapeU,[-L/2 L/2]);
%hold on;
%fplot(airshapeL,[-L/2 L/2]);

[XU YU] = fplot(airshapeU,[-L/2 L/2],800);
[XL YL] = fplot(airshapeL,[-L/2 L/2],800);

% Determine complex coordinates
ZU = double(XU + 1i*YU);
ZL = double(XL + 1i*YL);
Z = ZU;

for i = 1 : 2

    if i == 2
        Z = ZL;
    end

% Map to Zeta space
JP = 0.5*Z + ((Z.^2)/4 - c^2).^0.5;
JN = 0.5*Z - ((Z.^2)/4 - c^2).^0.5;

% Evaluate Complex Potential
WTP = U*exp(-1i*at) - U*a^2*exp(1i*at)./(JP - M).^2 + 1i*K./(JP - M);
WTN = U*exp(-1i*at) - U*a^2*exp(1i*at)./(JN - M).^2 + 1i*K./(JN - M);

% Map back to Z space
WP = WTP./(1 - c^2./JP.^2);
WN = WTN./(1 - c^2./JN.^2);

% Calculate Pressure Distribution
PP = Pinf + 0.5*rho*(U^2 - real(WP.^2));
CPP = (PP - Pinf)/(0.5*rho*U^2);

PN = Pinf + 0.5*rho*(U^2 - real(WN.^2));
CPN = (PN - Pinf)/(0.5*rho*U^2);


% Set plotting information
for k = 1 : size(XU,1)
    if XU(k) < 0
        WU(k) = WN(k);
        XN(k) = (XU(k)+L/2)/L;
    else
        WU(k) = WP(k);
        XN(k) = (XU(k)+L/2)/L;
    end
end

P = Pinf + 0.5*rho*(U^2 - real(WU.^2));
CP(:,i) = (P - Pinf)/(0.5*rho*U^2);

end

plot(XN,CP(:,1));
hold all;
plot(XN,CP(:,2));
set(gca,'YDir','reverse')
axis([0 1 -1 1]);

 

[danja]

Screen shots of what's happening:

With no camber, 6 deg angle of attack:

With camber of 4%, 6 deg angle of attack:

 

[danja]

Well I might just be an idiot. I was trying to map X-Y coordinates into the Zeta plane to determine the pressures at those coordinates. What I should've been doing was just getting the pressures along the curve of the circle in the Zeta plane and then translating them back to X-Y space. I think it's working now. Will report back.

 

[danja]

Does this look more correct?

With no camber, 6 deg angle of attack:

With 4% camber, 6 deg angle of attack:

 

[samkumar]

hi danja,
Can you please post the modified correct code?
thanks in advance.