User:Eas4200c.f08.gator.edwards/Wk7 Matlab

function hw7

% My and Mz in N*m, A's in m^2

Vy = 10000;
Vz = 5000;

NACA = '2415';
Ab = .0002;
Ae = .0002;
Ah = .0001;
Af = .0001;
E = 72000000000;
nu = .33;
ns = 1000;
c = .5;
y1 = .25;
y2 = .75;
t1 = .002;
t2 = .003;
ycen = 0;
zcen = 0;
Ab = 2 * Ab;
Ae = 2 * Ae;
Ah = 2 * Ah;
Af = 2 * Af;
sigUlt = 300000000;

A = [Af Ab Ae Ah];

% Convert inputs to percentages
NACA = num2str(NACA);
m = str2num(NACA(1)) / 100;
p = str2num(NACA(2)) / 10;
ta = str2num(NACA(3:4)) / 100;

% Divide y-axis into ns segments
y = 0:1/ns:1;

% Calculation of mean camber line (zc) and dz/dy
for i = 1:length(y)
    if i <= p * length(y)
        zc(i) = (m / p^2) * (2 * p * y(i) - y(i)^2);
        dzdy(i) = (m / p^2) * (2 * p - 2 * y(i));
    else
        zc(i) = (m / (1 - p)^2) * ((1 - 2 * p) + 2 * p * y(i) - y(i)^2);
        dzdy(i) = (m / (1 - p)^2) * (2 * p - 2 * y(i));
    end
end

% Thickness distribution calculation
zt = 5 * ta .* (.2969 .* sqrt(y) - .1260 .* y - .3516 .* y.^2 + .2843 .* y.^3 - .1015 .* y.^4);

% Calculation of angle from y-axis to tangent to mean camber line
theta = atan(dzdy);

% Calculate coordinates of points on airfoil surface
yu = y - zt .* sin(theta);
zu = zc + zt .* cos(theta);
yl = y + zt .* sin(theta);
zl = zc - zt .* cos(theta);

% Scale to cord length
zc = zc * c;
y = y * c;
yu = yu * c;
zu = zu * c;
yl = yl * c;
zl = zl * c;
y1 = y1 * c;
y2 = y2 * c;

% Combine upper and lower airfoil outlines
ytot = [yu fliplr(yl)];
ztot = [zu fliplr(zl)];

% Find points on surface that correspond to partitions
% Upper portion, first partition
diff = 1;
i = 0;
while diff > 0
    i = i + 1;
    diff = y1 - yu(i);
end
if abs(y1 - yu(i)) < abs(y1 - yu(i-1))
    p1(1) = i;
else
    p1(1) = i - 1;
end
% Upper portion, second partition
diff = 1;
j = 0;
while diff > 0
    j = j + 1;
    diff = y1 - yl(j);
end
if abs(y1 - yl(j)) < abs(y1 - yl(j-1))
    p1(2) = j;
else
    p1(2) = j - 1;
end
% Lower portion, first partition
diff = 1;
i = 0;
while diff > 0
    i = i + 1;
    diff = y2 - yu(i);
end
if abs(y2 - yu(i)) < abs(y2 - yu(i-1))
    p2(1) = i;
else
    p2(1) = i - 1;
end
% Lower portion, second partition
diff = 1;
j = 0;
while diff > 0
    j = j + 1;
    diff = y2 - yl(j);
end
if abs(y2 - yl(j)) < abs(y2 - yl(j-1))
    p2(2) = j;
else
    p2(2) = j - 1;
end

for i = 1:p2(1) - p1(1) + 1
    BF(i,1) = yu(p1(1) + i - 1);
    BF(i,2) = zu(p1(1) + i - 1);
end
for i = 1:p2(2) - p1(2) + 1
    EH(i,1) = yl(p1(2) + i - 1);
    EH(i,2) = zl(p1(2) + i - 1);
end

Abar = 0;
for i = 1:length(ytot)-1
    r = [ytot(i) - ycen, ztot(i) - zcen, 0];
    dl = [ytot(i+1) - ytot(i), ztot(i+1) - ztot(i), 0];
    dAt = cross(dl, r) / 2;
    Abar = Abar + dAt;
end

% Begin HW 7

% Part 1
ystr = [(yu(p2(1)) + yl(p2(2))) / 2, (yu(p1(1)) + yl(p1(2))) / 2];
zstr = [(zu(p2(1)) + zu(p1(1))) / 2, (zl(p1(2)) + zl(p2(2))) / 2];
ystr(3) = ystr(2);
ystr(4) = ystr(1);
zstr(3) = zstr(2);
zstr(4) = zstr(2);
zstr(2) = zstr(1);

ybar = (ystr(1) * Af + ystr(2) * Ab) / (Af + Ab);
zbar = sum(zstr) / 4;

Iy = sum(A .* ystr.^2);
Iz = sum(A .* zstr.^2);
Iyz = sum(A .* ystr .* zstr);

ky = Iy / (Iy * Iz - Iyz^2);
kz = Iz / (Iy * Iz - Iyz^2);
kyz = Iyz / (Iy * Iz - Iyz^2);

% 1 --> 4-1
% 2 --> 1-2
% 3 --> 2-3
% 4 --> 3-4
for i = 1:4
    Qy = A(i) * zstr(i);
    Qz = A(i) * ystr(i);
    qtild(i) = -(ky * Vy - kyz * Vz) * Qz - (kz * Vz - kyz * Vy) * Qy;
end
qtild(1) = 0;

q0 = (-qtild(2) * zstr(2) + qtild(4) * zstr(4)) / (2 * Abar(3))

EDU>> hw7

q0 =

  6.5052e+005

% Buckling
sumA = 0;
ybarA = 0;
zbarA = 0;
for i = 1:length(BF) - 1
    dA = sqrt((BF(i,1) - BF(i+1,1))^2 + (BF(i,2) - BF(i+1,2))^2) * t1;
    dybar = BF(i,1) + (BF(i,1) - BF(i+1,1)) / 2;
    dzbar = BF(i,2) + (BF(i,2) - BF(i+1,2)) / 2;
    sumA = sumA + dA;
    ybarA = ybarA + dA * dybar;
    zbarA = zbarA + dA * dzbar;
end
for i = 1:length(EH) - 1
    dA = sqrt((EH(i,1) - EH(i+1,1))^2 + (EH(i,2) - EH(i+1,2))^2) * t1;
    dybar = EH(i,1) + (EH(i,1) - EH(i+1,1)) / 2;
    dzbar = EH(i,2) + (EH(i,2) - EH(i+1,2)) / 2;
    sumA = sumA + dA;
    ybarA = ybarA + dA * dybar;
    zbarA = zbarA + dA * dzbar;
end
plateCen(1) = ybarA / sumA;
plateCen(2) = zbarA / sumA;
plateTotal = [BF; EH];

for i = 1:length(BF) - 1
    dA = sqrt((BF(i,1) - BF(i+1,1))^2 + (BF(i,2) - BF(i+1,2))^2) * t1;
    Iy = Iy + dA * ((BF(i,2) - BF(i+1,2)) / 2 - plateCen(2))^2;
    Iz = Iz + dA * ((BF(i,1) - BF(i+1,1)) / 2 - plateCen(1))^2;
    Iy = Iy + dA * ((BF(i,2) - BF(i+1,2)) / 2 - plateCen(2)) * ((BF(i,1) - BF(i+1,1)) / 2 - plateCen(1));
end
for i = 1:length(EH) - 1
    dA = sqrt((EH(i,1) - EH(i+1,1))^2 + (EH(i,2) - EH(i+1,2))^2) * t1;
    Iy = Iy + dA * ((EH(i,2) - EH(i+1,2)) / 2 - plateCen(2))^2;
    Iz = Iz + dA * ((EH(i,1) - EH(i+1,1)) / 2 - plateCen(1))^2;
    Iyz = Iyz + dA * ((EH(i,2) - EH(i+1,2)) / 2 - plateCen(2)) * ((EH(i,1) - EH(i+1,1)) / 2 - plateCen(1));
end



b = [0 0];
for i = 1:length(BF) - 1
    b(1) = b(1) + sqrt((BF(i,1) - BF(i+1,1))^2 + (BF(i,2) - BF(i+1,2))^2);
end
for i = 1:length(EH) - 1
    b(2) = b(2) + sqrt((EH(i,1) - EH(i+1,1))^2 + (EH(i,2) - EH(i+1,2))^2);
end

AR = .5:.1:2;
m = 2;
D = E * t1^3 / (12 * (1 - nu^2));
for i = 1:length(AR)
    kc = (m / AR(i) + AR(i) / m)^2;
    sigCrit(i) = kc * pi^2 * D / (b(1)^2 * t1);
    lambda = (AR(i)^4 / (81 * (1 + AR(i)^2)^4) * (1 + 81 / 625 + (81 / 25) * ((1 + AR(i)^2) / (1 + 9 * AR(i)^2))^2 + (81 / 25) * ((1 + AR(i)^2) / (9 + AR(i)^2))^2))^.5;
    kc2 = pi^2 / (32 * lambda * AR(i));
    sigCrit2(i) = kc2 * pi^2 * D / (b(1)^2 * t1);
end
subplot(2, 1, 1)
hold off
plot(AR, sigCrit, 'kx', AR, sigCrit2, 'rx')
hold on

for i = 1:length(AR)
    kc = (m / AR(i) + AR(i) / m)^2;
    sigCrit(i) = kc * pi^2 * D / (b(2)^2 * t1);
    lambda = (AR(i)^4 / (81 * (1 + AR(i)^2)^4) * (1 + 81 / 625 + (81 / 25) * ((1 + AR(i)^2) / (1 + 9 * AR(i)^2))^2 + (81 / 25) * ((1 + AR(i)^2) / (9 + AR(i)^2))^2))^.5;
    kc2 = pi^2 / (32 * lambda * AR(i));
    sigCrit2(i) = kc2 * pi^2 * D / (b(2)^2 * t1);
end

subplot(2, 1, 2)
hold off
plot(AR, sigCrit, 'kx', AR, sigCrit2, 'rx')
hold on