Instantly share code, notes, and snippets.

@ina111 /TR797.m
Last active Mar 24, 2018

Embed
What would you like to do?
% -----------
% Optimum Design of Nonplanar wings-Minimum Induced Drag
% for A Given Lift and Wing Root Bending Moment (NAL TR-797)
%
% Created by Takahiro Inagawa on 2013-2-28.
% Copyright (c) 2013 Takahiro Inagawa. All rights reserved.
% -----------
clear all
tic; %Computation time measurement
% ---- Constant ----
Lift = 1000; % Lift[N]
Uinf = 7.2; % Aircraft speed[m/s]
rho = 1.2; % Air density[kg/m3]
span = 30; % Span width[m]
le = span / 2; % length from root to tip
N = 100; % Partition number
beta = 0.85; % Coefficient related to the bendin moment
delta_S = linspace(le / N / 2, le / N / 2, N); % half of partition, Equal distance
% ---- Geometric Condition ----
% yz plane
% y axis is vertical, z axis is horizontal.
% fai is the local dihedral angel.
% y is the center of the partition.
% -----------------------------
y = linspace(delta_S(1), le - delta_S(N), N);
z = linspace(0, 0, N);
fai = linspace(0, 0, N);
% ---- Geometric variable number ----
% Variables required to seek the Q.
for i = 1:N
for j = 1:N
ydash(i,j) = (y(i) - y(j)) .* cos(fai(j)) + (z(i) - z(j)) .* sin(fai(j));
zdash(i,j) = - (y(i) - y(j)) .* sin(fai(j)) + (z(i) - z(j)) .* cos(fai(j));
y2dash(i,j) = (y(i) + y(j)) .* cos(fai(j)) - (z(i) - z(j)) .* sin(fai(j));
z2dash(i,j) = (y(i) + y(j)) .* sin(fai(j)) + (z(i) - z(j)) .* cos(fai(j));
R2puls(i,j) = (ydash(i,j) - delta_S(j)).^2 + zdash(i,j).^2;
R2minus(i,j) = (ydash(i,j) + delta_S(j)).^2 + zdash(i,j).^2;
Rdash2puls(i,j) = (y2dash(i,j) + delta_S(j)).^2 + z2dash(i,j).^2;
Rdash2minus(i,j) = (y2dash(i,j) - delta_S(j)).^2 + z2dash(i,j).^2;
end
end
for i = 1:N
for j = 1:N
Q(i,j) = - 1 / (2 * pi) *(((ydash(i,j) - delta_S(j)) / R2puls(i,j)
- (ydash(i,j) + delta_S(j)) / R2minus(i,j)) * cos(fai(i) - fai(j))
+ (zdash(i,j) / R2puls(i,j) - zdash(i,j) / R2minus(i,j)) * sin(fai(i) - fai(j))
+ ((y2dash(i,j) - delta_S(j)) / Rdash2minus(i,j)
- (y2dash(i,j) + delta_S(j)) / Rdash2puls(i,j)) * cos(fai(i) + fai(j))
+ (z2dash(i,j) / Rdash2minus(i,j) - z2dash(i,j) / Rdash2puls(i,j))
* sin(fai(i) + fai(j)));
end
end
% ---- Normalization ----
% Variables required to seek q.
delta_sigma = delta_S ./ le;
eta = y ./ le;
etadash = ydash ./ le;
eta2dash = y2dash ./ le;
zeta = z ./ le;
zetadash = zdash ./ le;
zeta2dash = z2dash ./ le;
gamma2puls = R2puls / (le^2);
gamma2minus = R2minus / (le^2);
gammadash2puls = Rdash2puls / (le^2);
gammadash2minus = Rdash2minus / (le^2);
for i = 1:N
for j = 1:N
q(i,j) = - 1 / (2 * pi) *(((etadash(i,j) - delta_sigma(j)) / gamma2puls(i,j)
- (etadash(i,j) + delta_sigma(j)) / gamma2minus(i,j)) * cos(fai(i) - fai(j))
+ (zetadash(i,j) / gamma2puls(i,j) - zetadash(i,j) / gamma2minus(i,j))
* sin(fai(i) - fai(j))
+ ((eta2dash(i,j) - delta_sigma(j)) / gammadash2minus(i,j)
- (eta2dash(i,j) + delta_sigma(j)) / gammadash2puls(i,j)) * cos(fai(i) + fai(j))
+ (zeta2dash(i,j) / gammadash2minus(i,j) - zeta2dash(i,j) / gammadash2puls(i,j))
* sin(fai(i) + fai(j)));
end
end
% ---- elliptic loading aerodynamic force ----
% Vn is Induced vertical velocisy.
% Vn is constant when elliptical circulation distribution.
BendingMoment_elpl = 2 / 3 / pi * le * Lift;
InducedDrag_elpl = Lift^2 / (2 * pi * rho * Uinf^2 * le^2);
Vn_elpl = Lift / (2 * pi * rho * Uinf * le^2);
% ---- Creating the optimization equation ----
c = 2 * cos(fai) .* delta_sigma;
b = 3 * pi / 2 *(eta .* cos(fai) + zeta .* sin(fai)) .* delta_sigma;
for i = 1:N
for j = 1:N
A(i,j) = pi * q(i,j) .* delta_sigma(i);
end
end
% ---- solve optimization problem ----
AAA = A + A';
ccc = -c;
bbb = -b;
LeftMat = vertcat([AAA ccc' bbb'], [ccc 0 0], [bbb 0 0]);
RightMat = vertcat(linspace(0,0,N)', -1, -beta);
SolveMat = LeftMat \ RightMat;
g = SolveMat(1:N);
mu = SolveMat(N+1:N+2); % Lagrange multiplier
% ---- After Solve ----
% efficient : span efficiency
% Gamma : Local circulation
% InducedDrag : Sum of induced drag
% Vn : Local induced vertical velocisy
% Lift0_elpl : Lift at root culcurated by area of the ellipse circulation distribution
% Gamma0_elpl : Circulation at root
% Gamma_elpl : Analytical ellipse circulation distribution @ beta = 1
% ---------------------
efficiency_inverse = g' * A * g;
efficiency = 1 / efficiency_inverse;
Gamma = g * Lift / (2 * le * rho * Uinf);
InducedDrag = efficiency_inverse * InducedDrag_elpl;
Local_Lift = 4 * rho * Uinf * Gamma' .* cos(fai);
Vn = zeros(1,N);
for i = 1:N
for j = 1:N
Vn(i) += Q(i,j) .* Gamma(j);
end
end
Local_InducedDrag = rho * Gamma' .* Vn;
% ---- Aerodynamic force when Elliptical cierculation distribution----
Lift0_elpl = 2 * Lift / pi / le;
Lift_elpl = 4 * Lift0_elpl * sqrt(1 - (y / le).^2);
Gamma0_elpl = Lift0_elpl / (rho * Uinf * cos(fai(1)));
Gamma_elpl = Gamma0_elpl * sqrt(1 - (y / le).^2);
Local_InducedDrag_elpl = 2 * rho * Gamma_elpl * Vn_elpl;
% ---- Bending Moment ----
Local_BendingMoment = zeros(1,N);
Local_BendingMoment_elpl = zeros(1,N);
for i = 1:N
tmp1 = 0;
tmp2 = 0;
for j = i:N
tmp1 += Local_Lift(j) * (y(j) - y(i));
tmp2 += Lift_elpl(j) * (y(j) - y(i));
end
Local_BendingMoment(i) = tmp1;
Local_BendingMoment_elpl(i) = tmp2;
end
% ---- Plot ----
% figure 1 : Circulation
% 2 : Lift
% 3 : Benging Moment
% 4 : Induced vertical velocity
% 5 : Induced Drag
figure(1)
plot(y, Gamma, y, Gamma_elpl);
xlabel('position [m]');ylabel('Gamma');
str = sprintf('\beta = %.2f', beta);
legend(str, 'ellipse circulation distribution');
title('Circulation');
xlim([0 ,max(y) * 1.05]);
grid on;
print("Circulation.jpg","-djpg","-S800,500","-r300");
figure(2)
plot(y, Local_Lift, y, Lift_elpl);
xlabel('position [m]'); ylabel('Lift[N]');
str = sprintf('\beta = %.2f', beta);
legend(str, 'ellipse circulation distribution');
title('Lift');
xlim([0 ,max(y) * 1.05]);
grid on;
print("Lift.jpg","-djpg","-S800,500","-r300");
figure(3)
plot(y, Local_BendingMoment, y, Local_BendingMoment_elpl);
xlabel('position [m]'); ylabel('Bending Moment [Nm]');
str = sprintf('beta = %.2f', beta);
legend(str, 'ellipse circulation distribution');
title('Bending Moment');
xlim([0, max(y) * 1.05]);
grid on;
print("Bending Moment.jpg","-djpg","-S800,500","-r300");
figure(4)
plot(y, linspace(Vn_elpl, Vn_elpl, N), y, Vn);
xlabel('position [m]'); ylabel('Induced vertical velocity [m/s]');
str = sprintf('beta = %.2f', beta);
legend(str, 'ellipse circulation distribution');
title('Induced vertical velocity');
xlim([0, max(y) * 1.05]);
grid on;
print("Induced vertical velocity.jpg","-djpg","-S800,500","-r300");
figure(5)
plot(y, Local_InducedDrag, y, Local_InducedDrag_elpl);
xlabel('position [m]'); ylabel('Induced Drag [N]');
str = sprintf('beta = %.2f', beta);
legend(str, 'ellipse circulation distribution');
title('Induced Drag');
xlim([0, max(y) * 1.05]);
grid on;
print("Induced Drag.jpg","-djpg","-S800,500","-r300");
% ---- Display Input and Output ----
disp('---- Input ----')
str = sprintf('Lift : %.1f', Lift);
disp(str)
str = sprintf('Aircraft speed : %.2f', Uinf);
disp(str)
str = sprintf('Wing span width : %.2f', span);
disp(str)
str = sprintf('Partition number : %d', N);
disp(str)
str = sprintf('beta : %.3f', beta);
disp(str)
disp('---- Output ----')
str = sprintf('Induced Drag : %.4f', InducedDrag);
disp(str)
str = sprintf('efficiency : %.4f', efficiency);
disp(str)
str = sprintf('');
disp(str)
disp('---- End ----')
toc
@ina111

This comment has been minimized.

Copy link
Owner

ina111 commented Sep 28, 2017

https://twitter.com/kikatyan/status/314421508313317379

・改行しても式が終わらないところは ... をつける.
・ダブルクォーテション "  -> シングルクオーテション '
・matlabで+=演算子はつかえない?ので a = a+xとなおす.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment