Commits

Nathan Collier committed 7eb20c6

bug fixes

  • Participants
  • Parent commits c52113b

Comments (0)

Files changed (4)

File demo/flapping_wing.py

     tmp = np.dot(tmp,R)
 
     # pitch
-    a = 0./180.*np.pi
+    a = -5./180.*np.pi
     R = np.zeros((3,3))
     R[1,1] = 1
     R[0,[0,2]] = [ np.cos(a), np.sin(a)]
 
 # Compute the wing, regular wing in x-y directions
 Lr,Lc = 1.,3.
-nr,nc = 3,3
+nr,nc = 7,11
 rho   = 1.225
 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,rho,symaxis=-1,symcord=0.,maxwake=10)
+grid = UVLMGrid(wing,flap,rho,symaxis=0,symcord=0.,maxwake=15)
 
 # Temporal stuffs
 N_period = 40
     t += dt
 
     # Plot
-    #grid.Plot("simple",i)
+    grid.Plot("simple",i)
 
-data = np.fromfile('./demo/data.dat',sep=' ').reshape((-1,4))
+data = np.fromfile('./data.dat',sep=' ').reshape((-1,4))
 
 ts = np.linspace(0,t,N_steps)
 
 plt.plot((data[:,0]-1.)*dt,data[:,2],'k--')
 plt.subplot(313)
 plt.title('Power')
-plt.plot(ts,grid.power,'g-')
+plt.plot(ts,np.asarray(grid.power)/(rho*0.5*6.*10.0**3.),'g-')
 plt.plot((data[:,0]-1.)*dt,data[:,3],'k--')
 plt.show()
 

File demo/turbine.py

 blade_width    = (blade[-1,:,0]-blade[0,:,0]).max()
 
 # setup blades
-wake_history = 15
+wake_history = 30
+rho = 1.225
 grids = []
-grids.append(UVLMGrid(blade,spin,maxwake=wake_history))
-grids.append(UVLMGrid(spin(blade,1.,3.),spin,maxwake=wake_history))
-grids.append(UVLMGrid(spin(blade,2.,3.),spin,maxwake=wake_history))
-
-"""
-# setup base
-pole_height = 1.2*blade_diameter
-pole_radius = 0.2*blade_width
-pole = np.zeros((5,45,3))
-ang = np.linspace(1.5*np.pi,2.0*np.pi,5)
-pole[:,:,0] = np.outer(pole_radius*np.sin(ang),np.ones(45)) + blade[:,:,0].max()+pole_radius+0.1
-pole[:,:,2] = np.outer(pole_radius*np.cos(ang),np.ones(45))
-pole[:,:,1] = np.outer(np.ones(5),np.linspace(-pole_height,0,45))
-grids.append(UVLMGrid(pole,nada,maxwake=wake_history))
-ang = np.linspace(1.5*np.pi,np.pi,5)
-pole[:,:,0] = np.outer(pole_radius*np.sin(ang),np.ones(45)) + blade[:,:,0].max()+pole_radius+0.1
-pole[:,:,2] = np.outer(pole_radius*np.cos(ang),np.ones(45))
-pole[:,:,1] = np.outer(np.ones(5),np.linspace(-pole_height,0,45))
-grids.append(UVLMGrid(pole,nada,maxwake=wake_history))
-"""
+grids.append(UVLMGrid(blade,spin,rho,maxwake=wake_history))
+grids.append(UVLMGrid(spin(blade,1.,3.),spin,rho,maxwake=wake_history))
+grids.append(UVLMGrid(spin(blade,2.,3.),spin,rho,maxwake=wake_history))
 
 # temporal data
 advance_ratio = 1.
     left.append(right[-1])
 N = right[-1]
 
-step0 = -1
-t = -(wake_history+1)*dt
-for step in range(N_period+wake_history):
+t = 0
+for step in range(N_period):
 
     print "Step %d, t = %.3f" % (step,t)
 
             Ai,bi = grids[i].Influence(grids[j])
             A[left[i]:right[i],left[j]:right[j]] += Ai
             b[left[i]:right[i]] += bi[:,0]
-        b[left[i]:right[i]] += grids[i].RelativeMotion(dt)[:,0]
+        b[left[i]:right[i]] += grids[i].RelativeMotion(dt,fsv)[:,0]
     c = np.linalg.solve(A,-b)
 
     # Load solution into grid
 
     # Convect the wake
     for grid in grids:
-        grid.ConvectWake(grids,dt,fsv=fsv)
+        grid.ConvectWake(grids,t,t_period,dt,fsv=fsv)
     t += dt
 
     # Plot
-    if t > 0.:
-        if step0 == -1: step0 = step
-        for i in range(len(grids)):
-            grids[i].Plot("turbine%d" % i,step-step0)
+    for i in range(len(grids)):
+        grids[i].Plot("turbine%d" % i,step)
 
-    

File src/pyfly/grid.py

 from pyfly.uvlm import base,influence
 from pyvtk2 import *
+from sds import CatmullClark
 import numpy as np
 
 __all__ = ['UVLMGrid']
 
 class UVLMGrid():
-    def __init__(self,grid,transform,rho,symaxis=-1,symcord=0.,maxwake=50):
+    def __init__(self,grid,transform,rho,symaxis=-1,symcord=0.,maxwake=50,cutoff=1e-9):
         self.grid0 = np.copy(grid)
         self.wake = np.zeros((maxwake+1,grid.shape[1],3))
         delta = (self.grid0[-1,:,:]-self.grid0[-2,:,:])
         self.Cdrag = []
         self.Clift = []
         self.rho = rho
+        self.cutoff = cutoff
 
     def Update(self,t,t_period):
         self.gridp = np.copy(self.grid)
         return
 
     def Influence(self,grid):
-        A = influence.Grid2OnGrid1(self.grid,self.normal,grid.grid,grid.symaxis,grid.symcord)
-        b = influence.NormalVelocityAtCenter(grid.wake,grid.wake_circ,self.grid,self.normal,grid.symaxis,grid.symcord)
+        A = influence.Grid2OnGrid1(self.grid,self.normal,grid.grid,grid.symaxis,grid.symcord,cutoff=self.cutoff)
+        b = influence.NormalVelocityAtCenter(grid.wake,grid.wake_circ,self.grid,self.normal,grid.symaxis,grid.symcord,cutoff=self.cutoff)
         b = b.reshape((-1,1))
         return A,b
 
-    def RelativeMotion(self,dt):
+    def RelativeMotion(self,dt,fsv=[0,0,0]):
         v_grid = (self.grid-self.gridp)/dt
         vel  = 0.25*(v_grid[:-1,:-1,:]+v_grid[1:,:-1,:]+v_grid[:-1,1:,:]+v_grid[1:,1:,:])
+        vel += np.asarray(fsv)
         vel *= self.normal
         vel = np.apply_along_axis(np.sum,2,vel)
         vel = vel.reshape((-1,1))
         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.wake,self.wake_circ)
+                                                               vrel,dt,self.rho,vref,self.wake,self.wake_circ,cutoff=self.cutoff)
         self.power.append(power)
         self.Cdrag.append(Cdrag)
         self.Clift.append(Clift)
-        print pressure
 
         # compute velocity at wake grid-points
         vel = np.zeros(self.wake.shape)
         pnts = self.wake
         for grid in grids:
-            vel += influence.VelocityAtPoints(grid.grid,grid.circ,pnts,grid.symaxis,grid.symcord)
-            vel += influence.VelocityAtPoints(grid.wake,grid.wake_circ,pnts,grid.symaxis,grid.symcord)
+            vel += influence.VelocityAtPoints(grid.grid,grid.circ,pnts,grid.symaxis,grid.symcord,cutoff=self.cutoff)
+            vel += influence.VelocityAtPoints(grid.wake,grid.wake_circ,pnts,grid.symaxis,grid.symcord,cutoff=self.cutoff)
         vel += np.asarray(fsv)
 
         # convect wake points

File src/pyfly/uvlm.f90

         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
-        print *,force
         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