In [1]:
## THIS script calculates the interpolation points for the delta wing
In [ ]:
#!!!!! Ideal treatment for hangover
In [2]:
## initialize the env
In [4]:
%pylab inline
import pylab as pl
import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
## coefficients. There are differrent coefficients for different r/c 9=(percentage)
# r/c = 0.0
# number of sampling points
i_Nx = 200
i_x = 0.15
# define coefficients
a = 0.11547005383792
d = 0.07003739672345
b = 0.12350964979191
c = -0.19567843344062
# define the initials coordinates of the leading edge
x_0 = 0.0
x_le = x_0
x_l = 0.15
# define the number of sampling points to be considered
x = np.zeros(i_Nx+1)
x[0] = 0.0
for i in range(1,i_Nx+1):
x[i] = x[i-1] + i_x/i_Nx
#print ('the vector x = '), x
# calculate psi values
psi = np.zeros(i_Nx+1)
psi[0] = 0.0
for i in range(1,i_Nx+1):
psi[i] = ( x[i]- x_0 ) / ( x_l )
#print ('The sampling points for the quatinty psi are = '), psi
#calculate the psi function for the upper surface
phi_of_psi_u = np.zeros(i_Nx+1)
for i in range(0,i_Nx+1):
phi_of_psi_u[i] = x_l * (a * math.sqrt(psi[i]) + b * psi[i] + c * psi[i]*psi[i] + d * psi[i] * psi[i] * psi[i] )
#print ('the function psi is = '),phi_of_psi_u
# mirror the upper surface
phi_of_psi_l = -1.0 * phi_of_psi_u
#print ('The lower surface phi of psi is = '),phi_of_psi_l
#________________________________________________________
# CREATE THE FLAT PLATE
#________________________________________________________
x_0_plate = 0.15
x_l_plate = 0.9 - x_0_plate
i_x_flat = x_l_plate
x_plate = np.zeros(i_Nx)
x_plate[0] = x[i_Nx] + i_x_flat/i_Nx
for i in range(1,i_Nx):
x_plate[i] = x_plate[i-1] + i_x_flat/i_Nx
#print x_plate
omega_psi_plate_u = np.zeros(i_Nx)
for i in range(0,i_Nx):
omega_psi_plate_u[i] = (d + (a * 3/8)) * x_l
omega_psi_plate_l = -1.0 * omega_psi_plate_u
#print omega_psi_plate_u
#print omega_psi_plate_l
#_______________________________________________________
#CONSTRUCT THE TRAILING EDGE CLOSURE REGION
#_______________________________________________________
#coefficients fom table A2
a_tr_edge = 0.0
d_tr_edge = 0.17000800036901
b_tr_edge = 3.0 * d_tr_edge
c_tr_edge = -1.0 * b_tr_edge
x_0_tr_edge = 0.9
x_l_tr_edge = 1.0 - x_0_tr_edge
i_x_tr_edge = x_l_tr_edge
x_tr_edge = np.zeros(i_Nx)
x_tr_edge[0] = x_plate[i_Nx-1] + i_x_tr_edge/i_Nx
for i in range(1,i_Nx):
x_tr_edge[i] = x_tr_edge[i-1] + i_x_tr_edge/i_Nx
#print x_tr_edge
psi_tr_edge = np.zeros(i_Nx)
psi_tr_edge[0] = 0.0
for i in range(1,i_Nx):
psi_tr_edge[i] = (x_tr_edge[i] - x_0_tr_edge)/x_l_tr_edge
#print ('the function psi for the trailing edge is = ') ,psi_tr_edge
phi_of_psi_tr_edge_u = np.zeros(i_Nx)
for i in range(0,i_Nx):
phi_of_psi_tr_edge_u[i] = x_l_tr_edge * (a_tr_edge * math.sqrt(psi_tr_edge[i_Nx-1-i]) + b_tr_edge * psi_tr_edge[i_Nx-1-i] + c_tr_edge * psi_tr_edge[i_Nx-1-i]*psi_tr_edge[i_Nx-1-i] + d_tr_edge * psi_tr_edge[i_Nx-1-i] * psi_tr_edge[i_Nx-1-i] * psi_tr_edge[i_Nx-1-i] )
#print ('The function phi of psi for the trailing edge upper surface is = '), phi_of_psi_tr_edge_u
phi_of_psi_tr_edge_l = -1.0 * phi_of_psi_tr_edge_u
#print ('The function phi of psi for the trailing edge lower surface is = '), phi_of_psi_tr_edge_l
#_______________________________________________________
#CONSTRUCT THE STING FAIRING
#coefficients for the stig fairing.....table A2
a_sting = 0.10040234847327
b_sting = 0.33279822819157
c_sting = -0.39554969598736
d_sting = 0.13603332984884
x_0_sting = 0.61057051
x_l_sting = 0.9797 - x_0_sting
i_x_sting = x_l_sting
x_sting = np.zeros(i_Nx+1)
x_sting[0] = x_0_sting #+ i_x_sting_f/i_Nx
for i in range(1,i_Nx+1):
x_sting[i] = x_sting[i-1] + i_x_sting/i_Nx
#print ('The coordinate x for the sting is = '), x_sting
psi_sting = np.zeros(i_Nx + 1)
for i in range(0,i_Nx+1):
psi_sting[i] = (x_sting[i]-x_0_sting)/x_l_sting
#print ('the function psi for the sting fairing is = '), psi_sting
phi_of_psi_sting_u = np.zeros(i_Nx+1)
phi_of_psi_sting_u[0] = 0.0
phi_of_psi_sting_u[i_Nx] = 0.06412
for i in range(1,i_Nx):
phi_of_psi_sting_u[i] = x_l_sting * ( a_sting * math.sqrt(psi_sting[i]) + b_sting * psi_sting[i] + c_sting * psi_sting[i] * psi_sting[i] + d_sting * psi_sting[i] * psi_sting[i] * psi_sting[i] )
phi_of_psi_sting_l = (-1.0)* phi_of_psi_sting_u
#print ('the function phi_of_psi_sting is = '), phi_of_psi_sting_u
#_______________________________________________________
#CONSTRUCT THE FORE FAIRING
## the profile for the fore fairing is the same at psi=0.06412 unitl the end of the sting from x/c = 0.9797 to x/c =1.758
## therefore it can be assumed that the fore fairing is a straight line. !!!!! Ideal treatment for hangover
x_0_fore = 0.9797
x_l_fore = 1.758 - x_0_fore
i_x_fore = x_l_fore
i_Nx_fore = 500
x_fore = np.zeros(i_Nx_fore+1)
x_fore[0] = x_0_fore
for i in range(1,i_Nx_fore+1):
x_fore[i] = x_fore[i-1] + i_x_fore/i_Nx_fore
#print ('the coordinate x for the fore fairing is = '), x_fore
phi_of_psi_fore_u = np.zeros(i_Nx_fore + 1)
for i in range(0,i_Nx_fore+1):
phi_of_psi_fore_u[i] = 0.06412
phi_of_psi_fore_l = (-1.0) * phi_of_psi_fore_u
#print ('the function psi for fore fairing is = '), phi_of_psi_fore
# CREATE ANOTHER PROFILE AND COPY IT ALONG THE SPAN
#_______________create the leading edge at x/c = 0.5 _________________________________________
x_0_cp = 0.5
x_le_cp = x_0
x_l_cp = 0.15
x_cp = np.zeros(i_Nx+1)
i_x_cp = x_l_cp
x_cp[0] = 5.0e-01
for i in range(1,i_Nx+1):
x_cp[i] = x_cp[i-1] + i_x_cp/i_Nx
#print ('the vector x_cp = '), x_cp
# calculate psi values
psi_cp = np.zeros(i_Nx+1)
for i in range(1,i_Nx+1):
psi_cp[i] = ( x_cp[i]- x_0_cp ) / ( x_l_cp )
#print ('The sampling points for the quatinty psi are = '), psi_cp
#calculate the psi function for the upper surface
phi_of_psi_cp_u = np.zeros(i_Nx+1)
for i in range(0,i_Nx+1):
phi_of_psi_cp_u[i] = x_l_cp * (a * math.sqrt(psi_cp[i]) + b * psi_cp[i] + c * psi_cp[i]*psi_cp[i] + d * psi_cp[i] * psi_cp[i] * psi_cp[i] )
#print ('the function psi is = '),phi_of_psi_cp_u
# mirror the upper surface
phi_of_psi_cp_l = -1.0 * phi_of_psi_cp_u
# ____________create the flat plate __________
x_0_plate_cp = 0.65
x_l_plate_cp = 0.9 - x_0_plate_cp
i_x_flat_cp = x_l_plate_cp
x_plate_cp = np.zeros(i_Nx)
x_plate_cp[0] = x_cp[i_Nx] + i_x_flat_cp/i_Nx
for i in range(1,i_Nx):
x_plate_cp[i] = x_plate_cp[i-1] + i_x_flat_cp/i_Nx
#print ('the coordinate x for the flat plate is '), x_plate_cp
omega_psi_plate_cp_u = np.zeros(i_Nx)
for i in range(0,i_Nx):
omega_psi_plate_cp_u[i] = (d + (a * 3/8)) * x_l
omega_psi_plate_cp_l = -1.0 * omega_psi_plate_cp_u
#print ('The function omega of psi for the flat plate is = '), omega_psi_plate_u
#print omega_psi_plate_l
#___________________________________________________
# ____________ create the trailing edge_____________
#coefficients fom table A2
a_tr_edge_cp = 0.0
d_tr_edge_cp = 0.17000800036901
b_tr_edge_cp = 3.0 * d_tr_edge
c_tr_edge_cp = -1.0 * b_tr_edge
x_0_tr_edge_cp = 0.9
x_l_tr_edge_cp = 1.0 - x_0_tr_edge
i_x_tr_edge_cp = x_l_tr_edge
x_tr_edge_cp = np.zeros(i_Nx)
x_tr_edge_cp[0] = x_plate_cp[i_Nx-1] + i_x_tr_edge_cp/i_Nx
for i in range(1,i_Nx):
x_tr_edge_cp[i] = x_tr_edge_cp[i-1] + i_x_tr_edge_cp/i_Nx
#print ('the x coordinate for the copied is profile is '), x_tr_edge_cp
psi_tr_edge_cp = np.zeros(i_Nx)
psi_tr_edge_cp[0] = 0.0
for i in range(1,i_Nx):
psi_tr_edge_cp[i] = (x_tr_edge_cp[i] - x_0_tr_edge_cp)/x_l_tr_edge_cp
#print ('the function psi for the trailing edge is = ') ,psi_tr_edge_cp
phi_of_psi_tr_edge_cp_u = np.zeros(i_Nx)
for i in range(0,i_Nx):
phi_of_psi_tr_edge_cp_u[i] = x_l_tr_edge_cp * (a_tr_edge_cp * math.sqrt(psi_tr_edge_cp[i_Nx-1-i]) + b_tr_edge_cp * psi_tr_edge_cp[i_Nx-1-i] + c_tr_edge_cp * psi_tr_edge_cp[i_Nx-1-i]*psi_tr_edge_cp[i_Nx-1-i] + d_tr_edge_cp * psi_tr_edge_cp[i_Nx-1-i] * psi_tr_edge_cp[i_Nx-1-i] * psi_tr_edge_cp[i_Nx-1-i] )
#print ('The function phi of psi for the trailing edge upper surface is = '), phi_of_psi_tr_edge_cp_u
phi_of_psi_tr_edge_cp_l = -1.0 * phi_of_psi_tr_edge_cp_u
#print ('The function phi of psi for the trailing edge lower surface is = '), phi_of_psi_tr_edge_cp_l
#_______________________________________________________
# WRITE THE DATA TO FILES
#_______________________________________________________write in matrix format for rhinoceros
#----------symmetry plane
# leading edge
LE = np.zeros(shape=(i_Nx+1,3))
for i in range(0,i_Nx+1):
LE[i,0] = x[i]
LE[i,1] = phi_of_psi_u[i]
LE[i,2] = 0.0
#print LE_u
#print x
# FLAT PLATE
FLAT_PLATE = np.zeros(shape=(i_Nx,3))
for i in range(0,i_Nx):
FLAT_PLATE[i,0] = x_plate[i]
FLAT_PLATE[i,1] = omega_psi_plate_u[i]
FLAT_PLATE[i,2] = 0.0
#print FLAT_PLATE
# TRAILING EDGE
TE = np.zeros(shape=(i_Nx,3))
for i in range(0,i_Nx):
TE[i,0] = x_tr_edge[i]
TE[i,1] = phi_of_psi_tr_edge_u[i]
TE[i,2] = 0.0
#print TE
# STING
STING = np.zeros(shape=(i_Nx+1,3))
for i in range(0,i_Nx+1):
STING[i,0] = x_sting[i]
STING[i,1] = phi_of_psi_sting_u[i]
STING[i,2] = 0.0
#print STING
# FORE
FORE = np.zeros(shape=(i_Nx_fore+1,3))
for i in range(0,i_Nx_fore+1):
FORE[i,0] = x_fore[i]
FORE[i,1] = phi_of_psi_fore_u[i]
FORE[i,2] = 0.0
#print FORE
# write in matrix format the profile copied at x/c 0.5 and exposed to y/c 0.233
# leading edge
LE_CP = np.zeros(shape=(i_Nx+1,3))
for i in range(0,i_Nx+1):
LE_CP[i,0] = x_cp[i]
LE_CP[i,1] = phi_of_psi_cp_u[i]
LE_CP[i,2] = 0.233
#print LE_CP
# flat plate
FLAT_PLATE_CP = np.zeros(shape=(i_Nx,3))
for i in range(0,i_Nx):
FLAT_PLATE_CP[i,0] = x_plate_cp[i]
FLAT_PLATE_CP[i,1] = omega_psi_plate_cp_u[i]
FLAT_PLATE_CP[i,2] = 0.233
#print FLAT_PLATE_CP
# trailing edge
TE_CP = np.zeros(shape=(i_Nx,3))
for i in range(0,i_Nx):
TE_CP[i,0] = x_tr_edge_cp[i]
TE_CP[i,1] = phi_of_psi_tr_edge_cp_u[i]
TE_CP[i,2] = 0.233
#print TE_CP
#-------------------SAVE THE DATA FOR RHINOCEROS----------------------
np.savetxt('LE_m.txt',LE,delimiter=' ')
np.savetxt('TE_m.txt',TE,delimiter=' ')
np.savetxt('FLAT_PLATE_m.txt',FLAT_PLATE,delimiter=' ')
np.savetxt('FORE_m.txt',FORE,delimiter=' ')
np.savetxt('STING_m.txt',STING,delimiter=' ')
np.savetxt('LE_CP_m.txt',LE_CP,delimiter=' ')
np.savetxt('TE_CP_m.txt',TE_CP,delimiter=' ')
np.savetxt('FLAT_PLATE_CP_m.txt',FLAT_PLATE_CP,delimiter=' ')
#_______________________________________________________
# PLOT THE POINTS
#_______________________________________________________
fig = plt.figure(figsize=(17, 2),dpi=1200, facecolor='w', edgecolor='k')
plt.plot(x,phi_of_psi_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=1)
#plt.plot(x,phi_of_psi_l,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_plate,omega_psi_plate_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_plate,omega_psi_plate_l,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_tr_edge,phi_of_psi_tr_edge_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_tr_edge,phi_of_psi_tr_edge_l,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_sting, phi_of_psi_sting_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_sting, phi_of_psi_sting_l,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_fore,phi_of_psi_fore_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_fore,phi_of_psi_fore_l,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_cp,phi_of_psi_cp_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_plate_cp,omega_psi_plate_cp_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_tr_edge_cp,phi_of_psi_tr_edge_cp_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#print x
#print psi
#print phi_of_psi_u
#print omega_psi_plate_u
plt.show()
print LE_CP
In [3]:
In [ ]: