danja
Mar5-11, 06:36 PM
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:
http://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]);
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:
http://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]);