# -*- 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
import math
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec # サブプロットのレイアウト用のgridspecモジュール
import cmath
class VLM:
def __init__(self, alpha, V, rho, Sref=None, b=None, AR=None, taperRatio=None, cr=None, sweep025=0, dihedral=0, iw=0, washout=0, airfoil='0012'):
self.alpha = math.radians(alpha) # 迎角 [deg]
self.V = V # 対気速度 [m/s]
self.rho = rho # 空気密度 [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.tbarc = 0.12 # 翼厚比 [-]
self.cMAC = None # 平均空力翼弦長 [m]
self.yMAC = None
self.x25MAC = None
self.airfoil = airfoil # NACA 4-digit airfoil
self.nx = 16
self.ny = 10 # セミスパン方向分割数
self.nSTA = 2 * self.ny + 1
self.nPanel = 2 * self.ny * self.nx
self.chord = None
self.x = np.zeros((self.nSTA) * self.nx + self.nSTA)
self.y = np.zeros((self.nSTA) * self.nx + self.nSTA)
self.z = np.zeros((self.nSTA) * 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.normV = np.zeros((self.nPanel,3))
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)
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
xc = (1-np.cos(np.pi/2*(np.arange(self.nx+1)/self.nx)))
yc = np.zeros_like(xc)
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))
for i in range(self.nSTA):
n1 = (i+0)*(self.nx+1)
n2 = (i+1)*(self.nx+1)
rot = Rotation.from_euler('Y', theta[i], degrees=True)
coord = rot.apply(np.stack([xc-0.25, yc, zc],axis=1))
self.x[n1:n2] = coord[:,0] * self.chord[i] + np.abs(ySTA[i]) * np.sin(self.swp)
self.y[n1:n2] = coord[:,1] * self.chord[i] + ySTA[i] * np.cos(self.dh)
self.z[n1:n2] = coord[:,2] * self.chord[i] + np.abs(ySTA[i]) * np.sin(self.dh)
#self.x -= np.max(self.x)
sinD = math.sin(self.dh)
cosD = math.cos(self.dh)
sinE = np.sin(self.dzdx)
cosE = np.cos(self.dzdx)
self.normV = np.array([-sinE*cosD,-cosE*sinD,cosE*cosD]).T
coord = np.stack([self.x,self.y,self.z],axis=1)
coord_cp = coord[:-1] * 0.25 + coord[1:] * 0.75
coord_sp = coord[:-1] * 0.75 + coord[1:] * 0.25
coord_diff = np.diff(coord, axis=0)
coord_norm = np.sqrt(np.sum(coord_diff**2,axis=1))
for i in range(2 * self.ny):
for j in range(self.nx):
n1 = i * self.nx + j
m1 = (i+0) * (self.nx + 1) + j
m2 = (i+1) * (self.nx + 1) + j
self.cp0[n1] = (coord_cp[m1] + coord_cp[m2]) / 2
self.sp1[n1] = coord_sp[m1]
self.sp2[n1] = coord_sp[m2]
self.dzdx[n1] = math.atan2((coord_diff[m1,2] + coord_diff[m2,2]) / 2, (coord_diff[m1,0] + coord_diff[m2,0]) / 2)
self.dS[n1] = (coord_norm[m1]+coord_norm[m2]) * np.abs(coord[m1,1]-coord[m2,1])/2
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)
Aij = (dy1*dz2-dz1*dy2)/r21 * ((dx21*dx1+dy21*dy1+dz21*dz1)/r1-(dx21*dx2+dy21*dy2+dz21*dz2)/r2)
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)
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 Aij, Bij, Cij
def solve_equation(self, Aij,Bij,Cij):
cosD = math.cos(self.dh)
matrix1 = Aij * self.normV[:,0].T + Bij * self.normV[:,1].T + Cij * self.normV[:,2].T # * 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):
wall = self.cp0 # + surf * self.cMAC * 0.001 * self.normV
Aij, Bij, Cij = self.calculate_Cij(wall)
ui = (1 / (4 * math.pi)) * (Aij @ self.circulation)
vi = (1 / (4 * math.pi)) * (Bij @ self.circulation)
wi = (1 / (4 * math.pi)) * (Cij @ self.circulation)
Vx = -self.V*math.cos(self.alpha)
Vz = -self.V*math.sin(self.alpha)
Vwall = np.sqrt((Vx+ui)**2 + vi**2 + (Vz+wi)**2)
self.Cp = 1 - (Vwall/self.V)**2
pn = -self.Cp * self.dS * np.cos(self.dzdx)
pt = self.Cp * self.dS * np.sin(self.dzdx)
x = wall[:,0]
z = wall[:,2]
cl = 2 * np.sum(pn * math.cos(self.alpha) - pt * math.sin(self.alpha))
cm = 2 * -np.sum(pn * x + pt * z)
cd = 2 * np.sum(pn * math.sin(self.alpha) + pt * math.cos(self.alpha))
cl /= self.Sref
cd /= self.Sref
cm /= (self.Sref * self.cMAC)
print(cl,cm,cd)
return cl, cd, cm
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 run(self):
# VLMの計算を実行
#self.calculate_chord()
self.calculate_geom()
self.calculate_mesh()
Aij, Bij, Cij = self.calculate_Cij(self.cp0)
#print(Cij)
self.solve_equation(Aij,Bij,Cij)
cl, cd, cm = self.calculate_section_coeff()
self.plot_3d()
fig = plt.figure(figsize=(16/2, 9/2)) # 図のサイズを指定して作成
ax = fig.add_subplot(111) # 3Dプロット用のサブプロットを追加
#ax.scatter(x.reshape(self.nx+1,-1)[:,:-1], y.reshape(self.nx+1,-1)[:,:-1], Cp.T, color='blue') # 3Dプロット
#ax.scatter(np.concatenate([np.flip(x.reshape(2*self.ny+1,-1)[:-1,:]),x.reshape(2*self.ny+1,-1)[:-1,1:]],axis=1), np.concatenate([np.flip(y.reshape(2*self.ny+1,-1)[:-1,:]),y.reshape(2*self.ny+1,-1)[:-1,1:]],axis=1), Cp_list, color='blue') # 3Dプロット
#ax.plot(self.x.reshape(-1,self.nx+1)[:-1,:-1].T, self.Cp[0].reshape(-1,self.nx).T, '.-', color='blue')
#ax.plot(self.x.reshape(-1,self.nx+1).T, self.z.reshape(-1,self.nx+1).T, '.-', color='red')
ax.plot(self.cp0[:,0].reshape(-1,self.nx).T, self.Cp.reshape(-1,self.nx).T, '.-', color='red')
#ax.plot(self.cp0[:,0], self.cp0[:,2], '.-', color='blue')
#wall = self.cp0 + 1 * 0.1 * self.normV
#ax.plot(wall[:,0], wall[:,2], '.-', color='red')
#wall = self.cp0 -1 * 0.1 * self.normV
#ax.plot(wall[:,0], wall[:,2], '.-', color='green')
#ax.set_aspect('equal')
# print(cl)
if __name__ == "__main__":
vlm = VLM(alpha=1, Sref=10, b=10, dihedral=5, sweep025=30, taperRatio=1, airfoil='4412', iw=0,washout=0, V=30, rho=1.225)
vlm.run()
