| # -*- 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