Skip to content

Instantly share code, notes, and snippets.

@ina111 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 Author

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
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.