In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import h5py
from importlib import reload
import PIVutils
import PODutils
import math
In [2]:
BoxSize = 10
X, Y = np.meshgrid(np.arange(-1*BoxSize, BoxSize+1), np.arange(-1*BoxSize, BoxSize+1))
U = np.zeros([2*BoxSize+1,2*BoxSize+1])
V = U.copy()
R = np.hypot(X, Y)
T = np.arctan2(Y,X)
#for i in range(X.shape[0]):
# T[i] = math.atan2(Y[i],X[i])
In [3]:
#Create Rankine Vortex
Circ = 30
r = 3
Ut = Circ*R/(2*np.pi*r**2)
#Ut[R>=r] = Circ/(2*np.pi*R[R>=r])
Ut[R>=r] = Circ/(2*np.pi*r) #make velocities constant outside core
#Now convert back to cartesian velocities
Uvort = Ut*np.sin(T)
Vvort = -1*Ut*np.cos(T)
In [4]:
plt.figure()
M = np.hypot(Uvort, Vvort)
Q = plt.quiver(X, Y, Uvort, Vvort, M, units='x', pivot='tip',headwidth=5, scale=1.5)
plt.axis('scaled')
Out[4]:
In [17]:
Rot = 45*np.pi/180*2 #Degrees of rotation
#position of stagnation point in polarcoords
rs = 5
Ts = 5/4*np.pi
xs = rs*np.cos(Ts) #shift in stagnation point in x
ys = rs*np.sin(Ts) #shift in stagnation point in y
StagStren = 2;
Xs = X-xs
Ys = Y-ys;
Ts = np.arctan2(Ys,Xs)
M = np.hypot(Xs, Ys)
U = M*np.cos(Ts-Rot)
V = -1*M*np.sin(Ts-Rot)
M = np.hypot(U, V)
Ustag = U/M*StagStren
Vstag = V/M*StagStren
Ustag[np.isnan(U)] = 0
Vstag[np.isnan(V)] = 0
In [18]:
plt.figure()
M = np.hypot(Ustag, Vstag)
Q = plt.quiver(X, Y, Ustag, Vstag, M, units='x', pivot='tip',headwidth=5, scale=1.5)
plt.axis('scaled')
Out[18]:
In [433]:
#Combine fields
Gvort_x = 5 #Radius of gaussian weighting function for vortex field
Gvort_y = Gvort_x
Gstag = 5 #Radius of gaussian weighting function for stagnation point field
Wvort = np.exp(-((X**2)/(2*Gvort_x)+(Y**2)/(2*Gvort_y)))
Wvort_inv = -1*Wvort+1 #invert the weightings so that only vortex appears at vortex location
Rstag = np.hypot(Xs, Ys)
Wstag = np.exp(-((X-xs)**2/(2*Gstag)+(Y-ys)**2/(2*Gstag)))
Wstag_inv = -1*Wstag+1
#U = (Wvort*Uvort+Wstag*Ustag)/(Wvort+Wstag)
#V = (Wvort*Vvort+Wstag*Vstag)/(Wvort+Wstag)
U = (Wvort*Wstag_inv*Uvort+Wstag*Wvort_inv*Ustag)/(Wvort*Wstag_inv+Wstag*Wvort_inv)
V = (Wvort*Wstag_inv*Vvort+Wstag*Wvort_inv*Vstag)/(Wvort*Wstag_inv+Wstag*Wvort_inv)
#U = (Wstag_inv*Uvort+Wvort_inv*Ustag)/(Wvort_inv+Wstag_inv)
#V = (Wstag_inv*Vvort+Wvort_inv*Vstag)/(Wvort_inv+Wstag_inv)
In [406]:
Uextra = Uvort
Uextra[:] = 0
Uextra[Y>0] = 0.5
U = U+Uextra
In [ ]:
In [434]:
plt.figure()
M = np.hypot(U, V)
Q = plt.quiver(X, Y, U, V, M, units='x', pivot='tip',headwidth=5, scale=1.5)
plt.axis('scaled')
Out[434]:
In [59]:
reload(PIVutils)
BoxSize = 10
Circ = 30
r = 3
rs = 7
Ts = 200
Rot = 45
StagStren = 2
Gvort = 5
Gstag = 5
Conv = 0
[U, V] = PIVutils.genHairpinField(BoxSize,Circ,r,rs,Ts,Rot,StagStren,Gvort,Gstag,Conv)
In [60]:
plt.figure()
M = np.hypot(U, V)
Q = plt.quiver(X, Y, U, V, M, units='x', pivot='tip',headwidth=5, scale=1.5)
plt.axis('scaled')
Out[60]:
In [103]:
reload(PIVutils)
x = np.arange(-0.05, 0.055,0.005)
y = np.arange(-0.05, 0.055,0.005)
X2, Y2 = np.meshgrid(x, y)
BoxSize = 10
Circ = 0.15
r = 0.01
rs = 0.04
Ts = 200
Rot = 45
StagStren = 3
Gvort = 0.015
Gstag = 0.015
Conv = 0
[U, V] = PIVutils.genHairpinField(BoxSize,Circ,r,rs,Ts,Rot,StagStren,Gvort,Gstag,Conv,x=x,y=y)
In [104]:
plt.figure()
M = np.hypot(U, V)
Q = plt.quiver(X2, Y2, U, V, M, units='x', pivot='tip',headwidth=5, scale=300)
plt.axis('scaled')
Out[104]:
In [97]:
from scipy.optimize import minimize
minimize?
In [38]:
x = np.arange(-0.05, 0.055,0.005)
(x.shape[0]-1)/2
print(x)
In [36]:
x
Out[36]:
In [ ]:
In [ ]: