# -*- coding: utf-8 -*-
"""VLM.ipynb
Automatically generated by Colab.
Original file is located at
https://colab.research.google.com/drive/1M5P0QMl-ksPZJ7MXuaPVcVsL7hJz84VI
"""
import numpy as np
import pandas as pd
from scipy.linalg import solve
from scipy.spatial.transform import Rotation
from scipy.optimize import minimize
import math
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec # サブプロットのレイアウト用のgridspecモジュール
from mpl_toolkits.mplot3d import Axes3D
class VLM:
def __init__(self, nx=16, ny=8, Sref=None, b=None, AR=None, taperRatio=None, cr=None, sweep025=0, dihedral=0, iw=0, washout=0, airfoil='0012'):
self.alpha = None # 迎角 [deg]
self.V = None # 対気速度 [m/s]
self.rho = None # 空気密度 [kg/m^3]
self.Sref = Sref # 面積 [m^2]
self.b = b # スパン [m]
self.AR = AR # アスペクト比 [m]
self.swp = math.radians(sweep025) # 1/4コードライン後退角
self.tr = taperRatio # テーパー比
self.dh = math.radians(dihedral) #上反角
self.cr = cr # 翼根コード長 [m]
self.iw = iw # 翼根取付角 [deg]
self.washout = washout # ねじり下げ [deg]
self.cMAC = None # 平均空力翼弦長 [m]
self.yMAC = None
self.x25MAC = None
self.chord = None
self.airfoil = airfoil # NACA 4-digit airfoil
self.CL = None
self.nx = nx
self.ny = ny # セミスパン方向分割数
self.nSTA = 2 * self.ny + 1
self.nPanel = 2 * self.ny * 2 * self.nx
self.x = np.zeros((self.nSTA) * 2 * self.nx + self.nSTA)
self.y = np.zeros((self.nSTA) * 2 * self.nx + self.nSTA)
self.z = np.zeros((self.nSTA) * 2 * self.nx + self.nSTA)
self.cp0 = np.zeros((self.nPanel,3))
self.sp1 = np.zeros((self.nPanel,3))
self.sp2 = np.zeros((self.nPanel,3))
self.dzdx = np.zeros(self.nPanel)
self.dS = np.zeros(self.nPanel)
self.Aij = None
self.Bij = None
self.Cij = None
self.dy = None
self.normV = None
self.circulation = None
self.Cp = None
def calculate_geom(self):
if self.Sref is None:
self.Sref = self.AR * self.b
if self.b is None:
self.b = np.sqrt(self.AR * self.S)
if self.AR is None:
self.AR = self.b ** 2 /self.Sref
if self.cr is None:
self.cr = self.Sref / (1 + self.tr) / (self.b / 2)
if self.tr is None:
self.tr = self.Sref / self.cr / (self.b / 2) - 1
self.cMAC = (2/3)*self.cr*(1+self.tr+self.tr**2)/(1+self.tr)
self.yMAC = (1/3)*(1+2*self.tr)/(1+self.tr)*(self.b/2)
self.x25MAC = self.cr * 0.25 + self.yMAC * np.tan(self.swp)
ySTA = self.b * (np.arange(self.nSTA)-self.ny)/(self.nSTA-1)
self.chord = self.cr * (1 + (self.tr - 1) * abs(ySTA)/(self.b/2))
return
def calculate_mesh(self):
ySTA = self.b * (np.arange(self.nSTA)-self.ny)/(self.nSTA-1)
theta = (self.iw + self.washout * abs(ySTA / (self.b / 2)))
m = float(self.airfoil[0])/100
p = float(self.airfoil[1])/10
t = float(self.airfoil[2:])/100
xc = (1-np.cos(np.pi/2*(np.arange(self.nx+1)/self.nx)))
zc = np.where(xc < p, m / max(p**2, 1e-5) * (2*p*xc - xc**2), m / ((1-p)**2) *((1-2*p)+2*p*xc-xc**2))
zt = 5 * t * (0.2969 * np.sqrt(xc) - 0.1260 * xc - 0.3516 * (xc**2) + 0.2843 * (xc**3) - 0.1015 * (xc**4))
dzcdx = np.where(xc < p, 2*m / max(p**2, 1e-5) * (p - xc), 2*m / ((1-p)**2) *(p-xc))
phi = np.arctan(dzcdx)
coords = []
for surf in [1,-1]:
xs = xc - surf * zt * np.sin(phi)
ys = np.zeros_like(xs)
zs = zc + surf * zt * np.cos(phi)
coord = []
for i in range(self.nSTA):
n1 = (i + 0) * (2 * self.nx + 1)
n2 = (i + 1) * (2 * self.nx + 1)
rot = Rotation.from_euler('Y', theta[i], degrees=True)
airfoil = rot.apply(np.stack([xs - 0.25, ys, zs], axis = 1))
airfoil[:,0] = airfoil[:,0] * self.chord[i] + np.abs(ySTA[i]) * np.sin(self.swp)
airfoil[:,1] = airfoil[:,1] * self.chord[i] + ySTA[i] * np.cos(self.dh)
airfoil[:,2] = airfoil[:,2] * self.chord[i] + np.abs(ySTA[i]) * np.sin(self.dh)
coord.append(airfoil)
coords.append(coord)
coords = np.array(coords)
coords[:,:,:,0] -= np.min(coords[:,:,:,0])
self.x, self.y, self.z = coords.reshape(-1,3).T
coord_cp = coords[:,:,:-1,:] * 0.25 + coords[:,:,1:,:] * 0.75
coord_sp = coords[:,:,:-1,:] * 0.75 + coords[:,:,1:,:] * 0.25
coord_diff = np.diff(coords, axis=2)
coord_norm = np.sqrt(np.sum(coord_diff ** 2, axis=3))
self.cp0 = ((coord_cp[:,1:,:,:] + coord_cp[:,:-1,:,:]) / 2).reshape(-1,3)
self.sp1 = coord_sp[:,:-1,:,:].reshape(-1,3)
self.sp2 = coord_sp[:,1:,:,:].reshape(-1,3)
self.dzdx = np.arctan((coord_diff[:,1:,:,2] + coord_diff[:,:-1,:,2]) / (coord_diff[:,1:,:,0] + coord_diff[:,:-1,:,0])).reshape(-1)
#self.dS = ((coord_norm[:,1:,:]+coord_norm[:,:-1,:])/2 * np.abs(coords[:,1:,0,1]-coords[:,:-1,0,1]).reshape(2,-1,1)).reshape(-1)
self.dy = self.sp2[:,1]-self.sp1[:,1]
sinD = np.sin(np.sign(self.cp0[:,1])*self.dh)
cosD = np.cos(self.dh)
sinE = np.sin(self.dzdx)
cosE = np.cos(self.dzdx)
self.normV = np.array([-sinE*cosD,-cosE*sinD,cosE*cosD]).T
#fig = plt.figure() # 図のサイズを指定して作成
#ax = fig.add_subplot(111) # 3Dプロット用のサブプロットを追加
#ax.plot(self.cp0[:,0].reshape(2*self.ny,-1), self.cp0[:,1].reshape(2*self.ny,-1), color='blue') # 3Dプロット
#ax.plot(self.cp0[:,0].reshape(2*self.ny,-1).T, self.dzdx.reshape(2*self.ny,-1).T, color='blue') # 3Dプロット
#ax.set_aspect('equal')
#fig.show()
return
def calculate_Cij(self,cp):
csp1 = cp.reshape(-1,1,3) - self.sp1.reshape(1,-1,3)
csp2 = cp.reshape(-1,1,3) - self.sp2.reshape(1,-1,3)
sp21 = self.sp2 - self.sp1
dx1, dy1, dz1 = csp1[:,:,0], csp1[:,:,1], csp1[:,:,2]
dx2, dy2, dz2 = csp2[:,:,0], csp2[:,:,1], csp2[:,:,2]
dx21, dy21, dz21 = sp21[:,0], sp21[:,1], sp21[:,2]
r1 = np.sqrt(dx1 ** 2 + dy1 ** 2 + dz1 ** 2)
r2 = np.sqrt(dx2 ** 2 + dy2 ** 2 + dz2 ** 2)
r21 = ((dy1*dz2-dz1*dy2)**2+(dz1*dx2-dx1*dz2)**2+(dx1*dy2-dy1*dx2)**2)
self.Aij = (dy1*dz2-dz1*dy2)/r21 * ((dx21*dx1+dy21*dy1+dz21*dz1)/r1-(dx21*dx2+dy21*dy2+dz21*dz2)/r2)
self.Bij = (dz1*dx2-dx1*dz2)/r21 * ((dx21*dx1+dy21*dy1+dz21*dz1)/r1-(dx21*dx2+dy21*dy2+dz21*dz2)/r2) + dz1/(dz1**2+dy1**2)*(1+dx1/r1) - dz2/(dz2**2+dy2**2)*(1+dx2/r2)
self.Cij = (dx1*dy2-dy1*dx2)/r21 * ((dx21*dx1+dy21*dy1+dz21*dz1)/r1-(dx21*dx2+dy21*dy2+dz21*dz2)/r2) - dy1/(dz1**2+dy1**2)*(1+dx1/r1) + dy2/(dz2**2+dy2**2)*(1+dx2/r2)
return
def solve_equation(self):
cosD = np.cos(self.dh)
matrix1 = self.Aij * self.normV[:,0].reshape(-1,1) + self.Bij * self.normV[:,1].reshape(-1,1) + self.Cij * self.normV[:,2].reshape(-1,1) # * self.cMAC
matrix2 = -4 * np.pi * self.V * np.sin(self.alpha - self.dzdx) * cosD
self.circulation = solve(matrix1, matrix2)
return
def calculate_section_coeff(self):
dh = np.sign(self.cp0[:,1]) * self.dh
ui = (1 / (4 * math.pi)) * (self.Aij @ self.circulation)
vi = (1 / (4 * math.pi)) * (self.Bij @ self.circulation)
wi = (1 / (4 * math.pi)) * (self.Cij @ self.circulation)
Vx = self.V * np.cos(self.alpha)
Vz = self.V * np.sin(self.alpha)
Vwall = np.sqrt((Vx + ui) ** 2 + vi ** 2 + (Vz + wi) ** 2)
self.Cp = 1 - (Vwall/self.V)**2
pn = self.rho * Vwall * self.circulation * self.dy * np.cos(self.dzdx) * np.cos(dh)
pt = -self.rho * Vwall * self.circulation * self.dy * np.sin(self.dzdx)
x = self.cp0[:,0] - self.x25MAC
z = self.cp0[:,2]
CL = np.sum(pn * np.cos(self.alpha) - pt * np.sin(self.alpha))
Cm025MAC = -np.sum(pn * x + pt * z)
CD = np.sum(pn * np.sin(self.alpha) + pt * np.cos(self.alpha))
CL /= (0.5 * self.rho * self.V * self.V * self.Sref)
CD /= (0.5 * self.rho * self.V * self.V * self.Sref)
Cm025MAC /= (0.5 * self.rho * self.V * self.V * self.Sref * self.cMAC)
fig = plt.figure() # 図のサイズを指定して作成
ax = fig.add_subplot(311) # 3Dプロット用のサブプロットを追加
ax.plot(np.arange(self.dzdx.shape[0]).reshape(4*self.ny,-1).T, Vwall.reshape(4*self.ny,-1).T, color='black', label='Vwall') # 3Dプロット
ax.plot(np.arange(self.dzdx.shape[0]).reshape(4*self.ny,-1).T, ui.reshape(4*self.ny,-1).T, color='blue', label='ui') # 3Dプロット
ax.plot(np.arange(self.dzdx.shape[0]).reshape(4*self.ny,-1).T, wi.reshape(4*self.ny,-1).T, color='red', label='wi') # 3Dプロット
ax = fig.add_subplot(312) # 3Dプロット用のサブプロットを追加
ax.plot(np.arange(self.dzdx.shape[0]).reshape(4*self.ny,-1).T, np.arctan((Vz+wi)/(Vx+ui)).reshape(4*self.ny,-1).T, color='blue', label='Vwall') # 3Dプロット
ax.plot(np.arange(self.dzdx.shape[0]).reshape(4*self.ny,-1).T, self.dzdx.reshape(4*self.ny,-1).T, color='red', label='dzdx') # 3Dプロット
ax = fig.add_subplot(313) # 3Dプロット用のサブプロットを追加
ax.plot(np.arange(self.dzdx.shape[0]).reshape(4*self.ny,-1).T, self.circulation.reshape(4*self.ny,-1).T, color='blue', label='circulation') # 3Dプロット
#ax.plot(np.arange(self.dzdx.shape[0]).reshape(4*self.ny,-1).T, Vwall.reshape(4*self.ny,-1).T, color='blue') # 3Dプロット
#ax.plot(self.cp0[:,0].reshape(2*self.ny,-1).T, self.dzdx.reshape(2*self.ny,-1).T, color='blue') # 3Dプロット
#ax.set_aspect('equal')
fig.show()
return CL, CD, Cm025MAC
def plot_3d(self):
fig = plt.figure() # 図のサイズを指定して作成
ax = fig.add_subplot(111, projection='3d') # 3Dプロット用のサブプロットを追加
ax.plot_wireframe(self.cp0[:,0].reshape(2*self.ny,-1), self.cp0[:,1].reshape(2*self.ny,-1), self.cp0[:,2].reshape(2*self.ny,-1), color='blue') # 3Dプロット
#ax.plot_wireframe(self.cp0[:,0].reshape(2*self.ny,-1), self.cp0[:,1].reshape(2*self.ny,-1), self.circulation.reshape(2*self.ny,-1), color='red') # 3Dプロット
#ax.plot_wireframe(self.cp0[:,0].reshape(2*self.ny,-1), self.cp0[:,1].reshape(2*self.ny,-1), self.Cp.reshape(2*self.ny,-1), color='red') # 3Dプロット
ax.plot_wireframe(self.x.reshape(2*self.ny+1,-1), self.y.reshape(2*self.ny+1,-1), self.z.reshape(2*self.ny+1,-1), color='black') # 3Dプロット
ax.set_aspect('equal')
def plot_2d(self):
xmax = np.ceil(np.max(self.x))
ymax = np.ceil(np.max(self.y))
zmax = np.ceil(np.max(self.z))
fig = plt.figure() # 図を作成
gs = gridspec.GridSpec(nrows=2, ncols=2, width_ratios=[2*ymax+2, zmax+2], height_ratios=[xmax+2, zmax+2]) # サブプロットのレイアウトを設定
# 右側のサブプロットを設定
ax0 = fig.add_subplot(gs[0])
ax0.set_aspect('equal') # アスペクト比を1:10に設定
# X
ax0.set_xlabel('Y') # X軸のラベルを設定
ax0.set_xlim(-ymax-1, ymax+1) # X軸の範囲を設定
ax0.set_xticks(np.arange(-ymax-1, ymax+2)) # X軸の目盛りを設定
#ax0.invert_xaxis() # X軸を反転
# Y
ax0.set_ylabel('X') # X軸のラベルを設定
ax0.set_ylim(-1, xmax+1) # X軸の範囲を設定
ax0.set_yticks(np.arange(-1, xmax+2)) # X軸の目盛りを設定
# 右側のサブプロットを設定
ax2 = fig.add_subplot(gs[2])
ax2.set_aspect(aspect=1) # アスペクト比を1:10に設定
# X
ax2.set_xlabel('Y') # X軸のラベルを設定
ax2.set_xlim(-ymax-1, ymax+1) # X軸の範囲を設定
ax2.set_xticks(np.arange(-ymax-1, ymax+2)) # X軸の目盛りを設定
# Y
ax2.set_ylabel('Z') # X軸のラベルを設定
ax2.set_ylim(-1, zmax+1) # X軸の範囲を設定
ax2.set_yticks(np.arange(-1, zmax+2)) # X軸の目盛りを設定
# 右側のサブプロットを設定
ax1 = fig.add_subplot(gs[1])
ax1.set_aspect(aspect=1) # アスペクト比を1:10に設定
# X
ax1.set_ylabel('X') # X軸のラベルを設定
ax1.set_ylim(-1, xmax+1) # X軸の範囲を設定
ax1.set_yticks(np.arange(-1, xmax+2)) # X軸の目盛りを設定
# Y
ax1.set_xlabel('Z') # X軸のラベルを設定
ax1.set_xlim(-1, zmax+1) # X軸の範囲を設定
ax1.set_xticks(np.arange(-1, zmax+2)) # X軸の目盛りを設定
ax1.invert_xaxis() # X軸を反転
ax0.plot(self.y.reshape(2,-1,self.nx+1)[0].T,self.x.reshape(2,-1,self.nx+1)[0].T,'-',color='black')
ax0.plot(self.y.reshape(2,-1,self.nx+1)[0] ,self.x.reshape(2,-1,self.nx+1)[0] ,'-',color='black')
ax0.plot(self.y.reshape(2,-1,self.nx+1)[1].T,self.x.reshape(2,-1,self.nx+1)[1].T,'-',color='black')
ax0.plot(self.y.reshape(2,-1,self.nx+1)[1] ,self.x.reshape(2,-1,self.nx+1)[1] ,'-',color='black')
ax0.plot(0,self.x25MAC,'*',color='red')
ax1.plot(self.z.reshape(2,-1,self.nx+1)[0].T,self.x.reshape(2,-1,self.nx+1)[0].T,'-',color='black')
ax1.plot(self.z.reshape(2,-1,self.nx+1)[0] ,self.x.reshape(2,-1,self.nx+1)[0] ,'-',color='black')
ax1.plot(self.z.reshape(2,-1,self.nx+1)[1].T,self.x.reshape(2,-1,self.nx+1)[1].T,'-',color='black')
ax1.plot(self.z.reshape(2,-1,self.nx+1)[1] ,self.x.reshape(2,-1,self.nx+1)[1] ,'-',color='black')
ax2.plot(self.y.reshape(2,-1,self.nx+1)[0].T,self.z.reshape(2,-1,self.nx+1)[0].T,'-',color='black')
ax2.plot(self.y.reshape(2,-1,self.nx+1)[0] ,self.z.reshape(2,-1,self.nx+1)[0] ,'-',color='black')
ax2.plot(self.y.reshape(2,-1,self.nx+1)[1].T,self.z.reshape(2,-1,self.nx+1)[1].T,'-',color='black')
ax2.plot(self.y.reshape(2,-1,self.nx+1)[1] ,self.z.reshape(2,-1,self.nx+1)[1] ,'-',color='black')
# レイアウトの調整
ax0.grid() # グリッドを表示
ax1.grid() # グリッドを表示
ax2.grid() # グリッドを表示
plt.tight_layout() # レイアウトを自動調整
plt.savefig('graph_position.png')
plt.show() # 図を表示
def _objfunc(self, x):
self.alpha = np.radians(x)
self.solve_equation()
CL, CD, Cm025MAC = self.calculate_section_coeff()
return np.abs(CL-self.CL)
def pre(self):
# VLMの計算を実行
#self.calculate_chord()
self.calculate_geom()
self.calculate_mesh()
self.calculate_Cij(self.cp0)
#print(Cij)
def run(self, V, rho=1.225, alpha=None, CL=None):
self.V = V
self.rho = rho
self.CL = CL
if alpha is not None:
self.alpha = np.radians(alpha)
else:
initial_guess = 0
result = minimize(fun=self._objfunc, x0=initial_guess)
self.alpha = np.radians(result.x[0])
self.solve_equation()
CL, CD, Cm025MAC = self.calculate_section_coeff()
return np.degrees(self.alpha), CL, CD, Cm025MAC
if __name__ == "__main__":
vlm = VLM(nx=16, ny=4, Sref=10, b=5, dihedral=0, sweep025=30, taperRatio=1, airfoil='4412', iw=0,washout=0)
vlm.pre()
alpha, CL, CD, Cm = vlm.run(V=30, alpha=0.001)
vlm.plot_2d()
#vlm.plot_3d()
print(alpha,CL,CD,Cm)
print(CL*CL/math.pi/vlm.AR)
