Commits

Nathan Collier committed b426b45

added post-processing

  • Participants
  • Parent commits 028b591

Comments (0)

Files changed (4)

File demo/simple.py

 from pyfly.uvlm import base,influence
 from pyfly.grid import UVLMGrid
 import numpy as np
+import pylab as plt
 
 def flap(grid,t,t_period):
     tmp = np.copy(grid)
 
 # Compute the wing, regular wing in x-y directions
 Lr,Lc = 1.,3.
-nr,nc = 3,3
+nr,nc = 5,15
+rho   = 1.0 
 wing = np.zeros((nr,nc,3))
 wing[:,:,0] = np.outer(np.linspace(0,Lr,nr),np.ones(nc))
 wing[:,:,1] = np.outer(np.ones(nr),np.linspace(0,Lc,nc))
-grid = UVLMGrid(wing,flap,symaxis=-1,symcord=0.,maxwake=3)
+grid = UVLMGrid(wing,flap,rho,symaxis=-1,symcord=0.,maxwake=3)
 
 # Temporal stuffs
 N_period = 40
 Omega = 2.
 t_period = (2.*np.pi/Omega)
 dt = t_period/N_period
-N_steps = 4
+N_steps = 40
 np.set_printoptions(precision=15)
 t = 0.
 for i in range(N_steps):
     # Plot
     #grid.Plot("simple",i)
 
-    
+ts = np.linspace(0,t,N_steps)
+plt.subplot(311)
+plt.plot(ts,grid.Cdrag,'r-')
+plt.subplot(312)
+plt.plot(ts,grid.Clift,'b-')
+plt.subplot(313)
+plt.plot(ts,grid.power,'g-')
+plt.show()
+

File src/pyfly/grid.py

 __all__ = ['UVLMGrid']
 
 class UVLMGrid():
-    def __init__(self,grid,transform,symaxis=-1,symcord=0.,maxwake=50):
+    def __init__(self,grid,transform,rho,symaxis=-1,symcord=0.,maxwake=50):
         self.grid0 = np.copy(grid)
         self.wake = np.zeros((maxwake+1,grid.shape[1],3))
         delta = (self.grid0[-1,:,:]-self.grid0[-2,:,:])
             self.wake[i,:,:] = self.grid0[-1,:,:]+i*delta
         self.wake_circ = np.zeros((maxwake,grid.shape[1]-1))
         self.circ = np.zeros((grid.shape[0]-1,grid.shape[1]-1))
+        self.circp = self.circ
         self.grid  = self.grid0
         self.gridp = self.grid
         self.normal,self.area = base.ComputeCellData(self.grid)
         self.symaxis = symaxis
         self.symcord = symcord
         self.transform = transform
+        self.power = []
+        self.Cdrag = []
+        self.Clift = []
+        self.rho = rho
 
     def Update(self,t,t_period):
         self.gridp = np.copy(self.grid)
 
     def ConvectWake(self,grids,t,t_period,dt,fsv=[0,0,0]):
 
+        v_grid = (self.grid-self.gridp)/dt
+        vrel = 0.25*(v_grid[:-1,:-1,:]+v_grid[1:,:-1,:]+v_grid[:-1,1:,:]+v_grid[1:,1:,:])
+        vref = abs(vrel[...,0].mean()-fsv[0])
+        pressure,power,Cdrag,Clift = influence.AerodynamicData(self.grid,self.normal,self.area,self.circ,self.circp,vrel,dt,self.rho,vref)
+        self.power.append(power)
+        self.Cdrag.append(Cdrag)
+        self.Clift.append(Clift)
+
         # compute velocity at wake grid-points
         vel = np.zeros(self.wake.shape)
         pnts = self.wake
         # move circulations        
         self.wake_circ[1:,:] = np.copy(self.wake_circ[:-1,:])
         self.wake_circ[0,:] = np.copy(self.circ[-1,:])
+
+        # this will be wrong for multiple grids
+        self.circp = np.copy(self.circ)
+
         return
 
     def Plot(self,fhead,num):

File src/pyfly/uvlm.f90

   r2 = pt-x1
   aa = CrossProduct(r1,r2)
   bb = DotProduct(aa,aa)
-  out = 0.
+  out = 0.0d0
   if(bb .gt. cutoff) then
      cc = r1/sqrt(DotProduct(r1,r1)) - r2/sqrt(DotProduct(r2,r2))
      out = DotProduct(r0,cc)/bb *aa*0.07957747154594767d0 ! 1/(4PI)
   end do
 end subroutine NormalVelocityAtCenter
 
+subroutine AerodynamicData(nr,nc,grid,normals,areas,circ,circ0,vrel,dt,rho,vref,pressure,power,Cdrag,Clift)
+  use base
+  implicit none
+  integer(kind=4), intent(in)  :: nr,nc
+  real   (kind=8), intent(in)  :: grid(3,nc+1,nr+1),normals(3,nc,nr),areas(nc,nr),circ(nc,nr),circ0(nc,nr),vrel(3,nc,nr),rho,dt,vref
+  real   (kind=8), intent(out) :: pressure(nc,nr),Clift,Cdrag,power
+  real   (kind=8)              :: force(3),dirL(3),dirD(3),tangent(3),T1(3,3),T2(3,3),X(3,3)
+  real   (kind=8)              :: chord,lift,drag,const,v_v,v_h,aoa,circ_local,dfdt,sig,sig0,span,W_ind
+  integer(kind=4)              :: i,j
+  Clift = 0.0d0
+  Cdrag = 0.0d0
+  power = 0.0d0
+  if(abs(vref) .lt. 1.0d-12) then
+     const = 0.0d0
+  else
+     const = 4.0d0/(rho*vref**2.0d0)
+  end if
+  T2 = 0.0d0
+  do i=1,nr
+     do j=1,nc
+        tangent = 0.5d0*(grid(:,j,i+1)+grid(:,j+1,i+1))-0.5d0*(grid(:,j,i)+grid(:,j+1,i))
+        chord = sqrt(DotProduct(tangent,tangent))
+        tangent = tangent/chord
+        v_h = DotProduct(vrel(:,j,i),tangent)
+        v_v = DotProduct(vrel(:,j,i),normals(:,j,i))
+        if(abs(v_h) .lt. 1.0d-12) then
+           aoa = 0.0d0
+        else
+           aoa = atan(v_v/v_h)
+        end if
+        T1(:,1) = tangent(:)
+        T1(:,2) = normals(:,j,i)
+        T1(:,3) = -CrossProduct(tangent,normals(:,j,i)) ! minus due to C-ordering
+
+        T2(1,1) =  cos(aoa)
+        T2(1,2) =  sin(aoa)
+        T2(2,1) = -T2(1,2)
+        T2(2,2) =  T2(1,1)
+        T2(3,3) =  1.0d0
+
+        X = matmul(transpose(T1),matmul(T2,T1))
+        dirL = matmul(tangent,X)
+        dirD = matmul(normals(:,j,i),X)
+
+        if(i .eq. 1) then
+           sig  = 0.5d0*circ(j,i)
+           sig0 = 0.5d0*circ0(j,i)
+           circ_local = circ(j,i)
+        else
+           sig  = 0.5d0*(circ(j,i-1)+circ(j,i))
+           sig0 = 0.5d0*(circ0(j,i-1)+circ0(j,i))
+           circ_local = circ(j,i)-circ(j,i-1)
+        end if
+
+        W_ind = 1.0d0
+        dfdt = (sig-sig0)/dt
+        span = areas(j,i)/chord
+        lift = rho*span*(sqrt(DotProduct(vrel(:,j,i),vrel(:,j,i)))*circ_local+chord*dfdt)*cos(aoa)
+        drag = rho*span*(-W_ind*circ_local+chord*dfdt*sin(aoa))
+        force = lift*dirL + drag*dirD
+        pressure(j,i) = DotProduct(force,normals(:,j,i))/areas(j,i)
+        power = power + 2.0d0*pressure(j,i)*areas(j,i)*v_v
+        Cdrag = Cdrag + force(1)*const
+        Clift = Clift + force(3)*const
+     end do
+  end do
+end subroutine AerodynamicData
+
 end module influence

File src/pyfly/uvlm.pyf

        real   (kind=8), intent(c,out) :: velocity(npx,npy)
      end subroutine NormalVelocityAtCenter
 
+     subroutine AerodynamicData(nr,nc,grid,normals,areas,circ,circ0,vrel,dt,rho,vref,pressure,power,Cdrag,Clift)
+       implicit none
+       integer(kind=4), intent(hide)  :: nr,nc
+       real   (kind=8), intent(c,in)  :: grid(nr+1,nc+1,3),normals(nr,nc,3),areas(nr,nc),circ(nr,nc),circ0(nr,nc),vrel(nr,nc,3)
+       real   (kind=8), intent(in)    :: dt,rho,vref
+       real   (kind=8), intent(c,out) :: pressure(nr,nc)
+       real   (kind=8), intent(out)   :: Clift,Cdrag,power
+     end subroutine AerodynamicData
    end module influence
 
      !