Skip to content

Instantly share code, notes, and snippets.

@ina111 ina111/TR797.py
Created Mar 24, 2018

Embed
What would you like to do?
# -*- coding: utf-8 -*-
# -----------
# 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 2018-03-24.
# Copyright (c) 2018 Takahiro Inagawa. All rights reserved.
# -----------
import numpy as np
from numpy import sin, cos, tan, pi, sqrt
import matplotlib.pyplot as plt
class Wing:
def __init__(self):
self.lift = 1000 # Lift [N]
self.Uinf = 7.2 # Aircraft speed [m/s]
self.rho = 1.2 # Air density [kg/m3]
self.span = 30 # Span width [m]
self.le = self.span / 2 # length from root to tip
self.N = 100 # Partition number
self.beta = 0.85 # Coefficient related to the bending moment
self.delta_S = np.ones(self.N) * self.le / self.N / 2 # hafl of partition, Equal distance
self.phi = np.zeros(self.N) # local dihedral angle [rad]
def calc(self):
N = self.N
le = self.le
delta_S = self.delta_S
lift = self.lift
Uinf = self.Uinf
rho = self.rho
span = self.span
beta = self.beta
# ---- Geometric Condition ----
# yz plane
# y axis is vertical, z axis is horizontal.
# phi is the local dihedral angel.
# y is the center of the partition.
# -----------------------------
y = np.linspace(delta_S[0], le - delta_S[N-1], N)
z = np.zeros(N)
phi = self.phi
ydash = np.zeros((N,N))
zdash = np.zeros((N,N))
y2dash = np.zeros((N,N))
z2dash = np.zeros((N,N))
R2puls = np.zeros((N,N))
R2minus = np.zeros((N,N))
Rdash2puls = np.zeros((N,N))
Rdash2minus = np.zeros((N,N))
Q = np.zeros((N,N))
q = np.zeros((N,N))
# ---- Geometric variable number ----
# Variables required to seek the Q.
for i in range(N):
for j in range(N):
ydash[i,j] = (y[i]-y[j]) * cos(phi[j]) + (z[i]-z[j]) * sin(phi[j])
zdash[i,j] = (y[i]-y[j]) * sin(phi[j]) + (z[i]-z[j]) * cos(phi[j])
y2dash[i,j] = (y[i]+y[j]) * cos(phi[j]) - (z[i]-z[j]) * sin(phi[j])
z2dash[i,j] = (y[i]+y[j]) * sin(phi[j]) + (z[i]-z[j]) * cos(phi[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;
for i in range(N):
for j in range(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(phi[i] - phi[j]) \
+ (zdash[i,j] / R2puls[i,j] - zdash[i,j] / R2minus[i,j]) * sin(phi[i] - phi[j]) \
+ ((y2dash[i,j] - delta_S[j]) / Rdash2minus[i,j] \
- (y2dash[i,j] + delta_S[j]) / Rdash2puls[i,j]) * cos(phi[i] + phi[j]) \
+ (z2dash[i,j] / Rdash2minus[i,j] - z2dash[i,j] / Rdash2puls[i,j]) \
* sin(phi[i] + phi[j]))
# ---- 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 in range(N):
for j in range(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(phi[i] - phi[j]) \
+ (zetadash[i,j] / gamma2puls[i,j] - zetadash[i,j] / gamma2minus[i,j]) \
* sin(phi[i] - phi[j]) \
+ ((eta2dash[i,j] - delta_sigma[j]) / gammadash2minus[i,j] \
- (eta2dash[i,j] + delta_sigma[j]) / gammadash2puls[i,j]) * cos(phi[i] + phi[j]) \
+ (zeta2dash[i,j] / gammadash2minus[i,j] - zeta2dash[i,j] / gammadash2puls[i,j]) \
* sin(phi[i] + phi[j]))
# ---- elliptic loading aerodynamic force ----
# Vn is Induced vertical velocisy.
# Vn is constant when elliptical circulation distribution.
bending_moment_elpl = 2 / 3 / pi * le * lift
induced_drag_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(phi) * delta_sigma
b = 3 * pi / 2 *(eta * cos(phi) + zeta * sin(phi)) * delta_sigma
A = np.zeros((N,N))
for i in range(N):
for j in range(N):
A[i,j] = pi * q[i,j] * delta_sigma[i]
# ---- solve optimization problem ----
AAA = A + A.T
ccc = -c
cccT = ccc.reshape(N,1)
ccc0 = np.append(ccc, np.zeros(2)).reshape(1,N+2)
bbb = -b
bbbT = bbb.reshape(N,1)
bbb0 = np.append(bbb, np.zeros(2)).reshape(1,N+2)
AAAcb = np.concatenate((AAA, cccT, bbbT), axis=1)
# import pdb; pdb.set_trace()
left_matrix = np.concatenate((AAAcb, ccc0, bbb0), axis=0)
right_matrix = np.append(np.zeros(N), np.array([-1, -beta]))
solve_matrix = np.linalg.solve(left_matrix, right_matrix)
g = solve_matrix[0:N]
mu = solve_matrix[N: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 = np.dot(np.dot(g,A),g)
efficiency = 1 / efficiency_inverse
Gamma = g * lift / (2 * le * rho * Uinf)
induced_drag = efficiency_inverse * induced_drag_elpl
local_lift = 4 * rho * Uinf * Gamma.T * cos(phi)
Vn = np.zeros(N);
for i in range(N):
for j in range(N):
Vn[i] += Q[i,j] * Gamma[j]
local_induced_drag = 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(phi[0]))
Gamma_elpl = Gamma0_elpl * sqrt(1 - (y / le)**2)
local_induced_drag_elpl = 2 * rho * Gamma_elpl * Vn_elpl
# ---- Bending Moment ----
local_bending_moment = np.zeros(N);
local_bending_moment_elpl = np.zeros(N);
for i in range(N):
tmp1 = 0;
tmp2 = 0;
for j in range(i,N):
tmp1 += local_lift[j] * (y[j] - y[i]);
tmp2 += lift_elpl[j] * (y[j] - y[i]);
local_bending_moment[i] = tmp1;
local_bending_moment_elpl[i] = tmp2;
# ---- Display Input and Output ----
print("==== Input ====")
print("Lift\t\t\t: %.1f [N]" % (lift))
print("Aircraft speed\t\t: %.2f [m/s]" % (Uinf))
print("Wing span width\t\t: %.1f [m]" % (span))
print("Partition number\t: %d" % (N))
print("beta\t\t\t: %.2f [-]" % (beta))
print("==== Output ====")
print("Induced Drag\t\t: %.2f [N]" % (induced_drag))
print("efficiency\t\t: %.1f [%%]" % (efficiency * 100))
print("==== cf. ellipse circulation distribution ====")
print("Induced Drag\t\t: %.2f [N]" % (induced_drag_elpl))
# ---- Plot ----
plt.ion()
plt.close("all")
plt.figure()
plt.plot(y, Gamma, label="beta=%.2f" % (beta))
plt.plot(y, Gamma_elpl, label="ellipse circulation distribution")
plt.xlabel("span [m]");plt.ylabel("Circulation")
plt.grid();plt.legend()
plt.title("Circulation")
plt.figure()
plt.plot(y, local_lift, label="beta=%.2f" % (beta))
plt.plot(y, lift_elpl, label="ellipse circulation distribution")
plt.xlabel("span [m]");plt.ylabel("Lift [N]")
plt.grid();plt.legend()
plt.title("Lift")
plt.figure()
plt.plot(y, local_bending_moment, label="beta=%.2f" % (beta))
plt.plot(y, local_bending_moment_elpl, label="ellipse circulation distribution")
plt.xlabel("span [m]");plt.ylabel("Bending Moment [Nm]")
plt.grid();plt.legend()
plt.title("Bending Moment")
plt.figure()
plt.plot(y, Vn, label="beta=%.2f" % (beta))
plt.plot(y, np.ones(N) * Vn_elpl, label="ellipse circulation distribution")
plt.xlabel("span [m]");plt.ylabel("velocity [m/s]")
plt.grid();plt.legend()
plt.title("Induced vertical velocity")
plt.figure()
plt.plot(y, local_induced_drag, label="beta=%.2f" % (beta))
plt.plot(y, local_induced_drag_elpl, label="ellipse circulation distribution")
plt.xlabel("span [m]");plt.ylabel("Drag [N]")
plt.grid();plt.legend()
plt.title("Induced Drag")
# plt.show()
if __name__ == '__main__':
wing = Wing()
wing.calc()
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.