tmp

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



質問・感想・意見などあれば気軽にTwitterのDMかコメントお願いします!
スポンサーリンク