Commits

Nathan Collier committed 88b73b0

added demos, full turbine example

Comments (0)

Files changed (3)

+from pyfly.uvlm import base,influence
+from pyfly.grid import UVLMGrid
+import numpy as np
+
+def flap(grid,t,t_period):
+    tmp = np.copy(grid)
+
+    # flap angle
+    a = np.sin(2.*np.pi*t/t_period)*np.pi*0.25
+    R = np.zeros((3,3))
+    R[0,0] = 1
+    R[1,[1,2]] = [ np.cos(a), np.sin(a)]
+    R[2,[1,2]] = [-np.sin(a), np.cos(a)]
+    tmp = np.dot(tmp,R)
+
+    # pitch
+    a = 5./180.*np.pi
+    R = np.zeros((3,3))
+    R[1,1] = 1
+    R[0,[0,2]] = [ np.cos(a), np.sin(a)]
+    R[2,[0,2]] = [-np.sin(a), np.cos(a)]
+    tmp = np.dot(tmp,R)
+
+    # forward flight @ 10 m/s
+    tmp[:,:,0] -= 10.0*t
+    return tmp
+
+# Compute the wing, regular wing in x-y directions
+Lr,Lc = 2.,6.
+nr,nc = 6,10
+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=10)
+
+# Temporal stuffs
+N_period = 50
+h = np.sqrt((Lr/nr)**2+(Lc/nc)**2)
+dt = h/10.
+t_period = N_period*dt
+
+t = 1e-2
+for i in range(N_period):
+    
+    print "Step %d:" % i
+
+    # Setup grid
+    grid.Update(t,t_period)
+
+    # Setup/solve the linear system
+    A,b = grid.Influence(grid)
+    b += grid.RelativeMotion(dt)
+    c = np.linalg.solve(A,-b)
+    grid.circ = c.reshape(grid.circ.shape)
+
+    # Convect the wake
+    grid.ConvectWake([grid],dt)
+    t += dt
+
+    # Plot
+    grid.Plot("simple",i)
+
+    
+from pyfly.uvlm import base,influence
+from pyfly.grid import UVLMGrid
+from pickle import load
+import numpy as np
+
+def spin(grid,t,t_period):
+    tmp = np.copy(grid)
+    a = t/t_period*(2.*np.pi)
+    R = np.zeros((3,3))
+    R[0,0] = 1
+    R[1,[1,2]] = [ np.cos(a), np.sin(a)]
+    R[2,[1,2]] = [-np.sin(a), np.cos(a)]
+    tmp = np.dot(tmp,R)
+    return tmp
+
+def nada(grid,t,t_period):
+    return np.copy(grid)
+
+# read in blade
+f = open('blade.dat','rb')
+blade = load(f)
+f.close()
+blade_diameter = 2*(blade[:,-1,1]-blade[:,0,1]).max()
+blade_width    = (blade[-1,:,0]-blade[0,:,0]).max()
+
+# setup blades
+grids = []
+grids.append(UVLMGrid(blade,spin))
+grids.append(UVLMGrid(spin(blade,1.,3.),spin))
+grids.append(UVLMGrid(spin(blade,2.,3.),spin))
+
+# setup base
+pole_height = 1.2*blade_diameter
+pole_radius = 0.2*blade_width
+pole = np.zeros((10,30,3))
+ang = np.linspace(1.5*np.pi,2.0*np.pi,10)
+pole[:,:,0] = np.outer(pole_radius*np.sin(ang),np.ones(30)) + pole_radius + 0.5*(blade[-1,0,0]-blade[0,0,0])
+pole[:,:,2] = np.outer(pole_radius*np.cos(ang),np.ones(30))
+pole[:,:,1] = np.outer(np.ones(10),np.linspace(-pole_height,0,30))
+grids.append(UVLMGrid(pole,nada))
+ang = np.linspace(1.5*np.pi,np.pi,10)
+pole[:,:,0] = np.outer(pole_radius*np.sin(ang),np.ones(30)) + pole_radius + 0.5*(blade[-1,0,0]-blade[0,0,0])
+pole[:,:,2] = np.outer(pole_radius*np.cos(ang),np.ones(30))
+pole[:,:,1] = np.outer(np.ones(10),np.linspace(-pole_height,0,30))
+grids.append(UVLMGrid(pole,nada))
+
+# temporal data
+fsv = np.asarray([10,0,0])
+N_period = 75
+junk,area = base.ComputeCellData(grids[0].grid)
+h = np.sqrt(area.mean())
+dt = h/np.linalg.norm(fsv)
+t_period = N_period*dt
+wake_history = 15
+t = -wake_history*dt
+
+# solver setup
+left = [0]
+right = []
+for grid in grids:
+    right.append(left[-1] + grid.circ.shape[0]*grid.circ.shape[1])
+    left.append(right[-1])
+
+N = right[-1]
+step0 = -1
+for step in range(N_period+wake_history):
+
+    print "Step %d..." % step
+
+    # Setup grid
+    for grid in grids:
+        grid.Update(t,t_period)
+
+    # Setup/solve the linear system
+    A = np.zeros((N,N))
+    b = np.zeros(N)
+    for i in range(len(grids)):
+        for j in range(len(grids)):
+            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]
+    c = np.linalg.solve(A,-b)
+
+    # Load solution into grid
+    for i in range(len(grids)):
+        grids[i].circ = c[left[i]:right[i]].reshape(grids[i].circ.shape)
+
+    # Convect the wake
+    for grid in grids:
+        grid.ConvectWake(grids,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)
+
+    

test/simple.py

-from pyfly.uvlm import base,influence
-from pyfly.grid import UVLMGrid
-import numpy as np
-from time import time
-
-def flap(grid,t,t_period):
-    tmp = np.copy(grid)
-
-    # flap angle
-    a = np.sin(2.*np.pi*t/t_period)*np.pi*0.25
-    R = np.zeros((3,3))
-    R[0,0] = 1
-    R[1,[1,2]] = [ np.cos(a), np.sin(a)]
-    R[2,[1,2]] = [-np.sin(a), np.cos(a)]
-    tmp = np.dot(tmp,R)
-
-    # pitch
-    a = 5./180.*np.pi
-    R = np.zeros((3,3))
-    R[1,1] = 1
-    R[0,[0,2]] = [ np.cos(a), np.sin(a)]
-    R[2,[0,2]] = [-np.sin(a), np.cos(a)]
-    tmp = np.dot(tmp,R)
-
-    # forward flight @ 10 m/s
-    tmp[:,:,0] -= 10.0*t
-    return tmp
-
-# Compute the wing, regular wing in x-y directions
-Lr,Lc = 2.,6.
-nr,nc = 4,12
-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.)
-
-# Temporal stuffs
-N_period = 50
-h = np.sqrt((Lr/nr)**2+(Lc/nc)**2)
-dt = h/10.
-t_period = N_period*dt
-
-t = 1e-2
-for i in range(N_period):
-    
-    print "Step %d:" % i
-
-    # Setup grid
-    tic = time()
-    grid.Update(t,t_period)
-    print "\tUpdate: %.3e" % (time()-tic)
-
-    # Setup/solve the linear system
-    tic = time()
-    A,b = grid.Influence(grid)
-    b += grid.RelativeMotion(dt)
-    print "\tInfluence: %.3e" % (time()-tic)
-    tic = time()
-    c = np.linalg.solve(A,-b)
-    grid.circ = c.reshape(grid.circ.shape)
-    print "\tSolve: %.3e" % (time()-tic)
-
-    # Convect the wake
-    tic = time()
-    grid.ConvectWake([grid],dt)
-    t += dt
-    print "\tConvect: %.3e" % (time()-tic)
-
-    # Plot
-    tic = time()
-    grid.Plot("simple",i)
-    print "\tPlot: %.3e" % (time()-tic)
-