Biturbowned
Triple Kill!
Joined: 29 Nov 2006
Posts: 267
Location: here 'n there
|
 Posted: Mon Mar 19, 2007 3:56 pm Post subject: |
 |
|
| PC Oibo wrote: | | that tells us nothing. I want to see code mister. |
you want it, you got it:
function flatDeltaWing1(Minf,Sweep,AoA,Nmax)
% This program calculates the change in Cp across a flat delta wing
% (bottom - top) for linearized supersonic flow. It is based on the
% publication "Numerical Methods for the Design and Analysis of Wings at
% Supersonic Speeds" by Carlson and Miller, NASA Technical Note TN D-7713.
% The inputs are the free stream Mach number (Minf), the sweep angle
% (Sweep) in degrees, the angle of attack (AoA) in degrees, and Nmax which
% is the number of additional spanwise elements to the wing tip. (The
% total number of spanwise elements will be 2*Nmax+1).
clc;clear all;
help flatDeltaWing1
% Convert the angle of attack to radians
AoA=input('Input angle of attack in degrees ');
alpha = AoA*pi/180;
% Find the beta parameter
Minf=1.4;
beta = sqrt(Minf*Minf-1);
% Locate the trailing edge and hence find the number of chordwise elements
% required.
Sweep1 = 60;
Sweep2 = 45;
Nmax=30;
scale=Nmax/(4*tand(90-Sweep1)+2*tan(90-Sweep2));
S=scale*(16*tan(90-Sweep1)+2);
xTE = (tan(Sweep1*pi/180)*(4/6)+(2/6)*tan(Sweep2*pi/180))*Nmax/beta;
Lmax = ceil(xTE);
% Dimension the required arrays
A = zeros(Lmax+1,2*Nmax+1);
B = zeros(Lmax+1,2*Nmax+1);
C = zeros(Lmax+1,2*Nmax+1);
delCp = zeros(Lmax+1,2*Nmax+1);
delCpa = zeros(2*Nmax+1,1);
delCpb = zeros(2*Nmax+1,1);
% Calculate the leading-edge-element weighting factors (A)
% Calculate the trailing-edge-element weighting factors (B)
% Calculate the wing-tip-element weighting factors (C)
for i = 1:Lmax+1
for j = 1:2*Nmax+1
N = j - Nmax - 1;
xle = xLE(N,scale,beta,Sweep1,Sweep2);
if i - xle <= 0
A(i,j) = 0;
end
if (0 < i - xle && i - xle < 1)
A(i,j) = i - xle;
end
if i - xle >= 1
A(i,j) = 1;
end
if i - xTE >= 1
B(i,j) = 0;
end
if (0 < i - xTE && i - xTE < 1)
B(i,j) = 1 - (i - xTE);
end
if i - xTE <= 0
B(i,j) = 1;
end
if abs(N) == Nmax
C(i,j) = 0.5;
else
C(i,j) = 1;
end
end
end
% Solve for the change in Cp
for i = 1:Lmax
% Find the preliminary delta Cp for the ith column
for j = 1:2*Nmax+1
N = j - Nmax - 1;
if A(i,j) == 0
delCpa(j)=0;
else
delCpa(j) = delCpf(i,N,Nmax,scale,beta,alpha,A,B,C,delCp);
end
end
% Add the preliminary Cp's to the already known ones to get ready to
% calculate the "correction" Cp
for j = 1:2*Nmax +1
delCp(i,j) = delCpa(j);
end
% Calculate the "correction" Cp. Note: Here we are icluding an
% element beyond the trailing edge so that we can compute the
% "correction" Cp and average.
for j = 1:2*Nmax +1
N = j - Nmax - 1;
if A(i+1,j) == 0
delCpb(j) = 0;
else
delCpb(j) = delCpf(i+1,N,Nmax,scale,beta,alpha,A,B,C,delCp);
end
end
% Do the appropriate averaging
for j = 1:2*Nmax +1
N = j - Nmax - 1;
xle = xLE(N,scale,beta,Sweep1,Sweep2);
if i - xle <= 1
delCp(i,j) = (1 + A(i,j)/(1 + A(i,j)))*delCp(i,j)/2 + (A(i,j)/(1+A(i,j)))*delCpb(j)/2;
end
if i - xle > 1
delCp(i,j) = 3*delCp(i,j)/4 + 1*delCpb(j)/4;
end
Cl(i,j,beta,S,delCp,A,B,C,xTE,xLE)
end
end
% surf(delCp)
end
% -------------------------------------------------------------------------
function a = xLE(N,scale,beta,Sweep1,Sweep2)
% This function locates the leading edge as a function of spanwise cell
% element index, N.
if abs(N)<=scale*4*tand(90-Sweep1)
a = abs(N)*tan(Sweep1*pi/180.)/beta;
elseif abs(N)>scale*4*tand(90-Sweep1)
a = abs(N)*(tan(Sweep2*pi/180))/beta+(scale*4*tand(90-Sweep1));
end
end
% -------------------------------------------------------------------------
function b = delCpf(Lstar,Nstar,Nmax,scale,beta,alpha,A,B,C,delCp)
% This function calculates the approximate Cauchy Principle Value integral
% as a finite sum over the elements contained within the forward Mach
% cone.
b2 = 0;
for j = 1:2*Nmax+1
N = j - Nmax -1;
Lmax = Lstar-abs(Nstar-N);
if Lmax >=1
for i = 1:Lmax
if (i<=4*scale)
dzdxu=.028/(4*scale); %slope of upper surface
dzdxl=-dzdxu;
elseif(4*scale<i<=6*scale)
dzdxu=-.028/(2*scale); %slope of lower surface
dzdxl=-dzdxu;
end
b1 = (4*alpha-2*(dzdxl-dzdxu))/beta; %account for upper and lower surfaces
b3 = (((Lstar - i + 0.5)^2 - (Nstar - N - 0.5)^2)^0.5)/((Lstar - i + 0.5)*(Nstar - N - 0.5));
b3 = b3 - (((Lstar - i + 0.5)^2 - (Nstar - N + 0.5)^2)^0.5)/((Lstar - i + 0.5)*(Nstar - N + 0.5));
b3 = b3*A(i,j)*B(i,j)*C(i,j)*delCp(i,j);
b2 = b2 + b3;
end
end
end
b = b1 + b2/pi;
end
% -------------------------------------------------------------------------
function c = Cl(Lstar,Nstar,beta,S,delCp,A,B,C,xTE,xLE)
sum=0;
for Nstar=0:1:Nmax
for Lstar=xLE:1:xTE
sum=sum+(3/4)*delCp(Lstar,Nstar)+(1/4)*delCp(Lstar+1,Nstar)*A(Lstar,Nstar)*B(Lstar,Nstar)*C(Lstar,Nstar);
end
end
c=2/(beta*S)*sum
end
_________________
PS3::XB360 |
|
[/img]
Don't mind me, I'm just jealous. If I could only learn to do that shit, I wouldn't need to get famous 