Updated and marked up code.

This commit is contained in:
Anne de Jong 2018-12-16 13:08:34 +01:00
parent 51637605b7
commit f0ec52f51f
2 changed files with 104 additions and 101 deletions

117
blflow.py Executable file → Normal file
View File

@ -1,72 +1,83 @@
#!/usr/bin/python #!/usr/bin/python
# Boundary layer flow # Boundary layer flow
from numpy import * import numpy as np
import time import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('TkAgg')
# from matplotlib.pylab import *
import pylab as p
# import matplotlib.animation as animation
def K(t): #Forcing function
return (1-exp(-0.1*t))*cos(t)
s=10 def K(t):
#Define domain """
n=50 #Number of gridpoints Forcing function as a function of time.
y=linspace(0,1,n) """
return (1-np.exp(-0.1*t))*np.cos(t)
dy=y[1]-y[0]
dt=0.0005
l=(dt/(s**2*dy**2)) s = 10 # Shear wave number
n = 50 # Number of gridpoints in vertical direction
hnu=exp(-sqrt(1j)*s*y) # Grid
y = np.linspace(0, 1, n)
dy = y[1]-y[0]
dt = 0.0005
l = (dt/(s**2*dy**2))
hnu = np.exp(-np.sqrt(1j)*s*y)
# fnu=(1-1j)/s # fnu=(1-1j)/s
fnu=0 fnu = 0
def u_ex(tn): def u_ex(tn):
return (((1-hnu)/(1-fnu))*exp(1j*(tn))/1j).real """
Exact harmonic solution of the velocity field
"""
return (((1-hnu)/(1-fnu))*np.exp(1j*(tn))/1j).real
def u_np1(un,tn,dt):
Kn=K(tn) def u_np1(un, tn, dt):
unp1=un Kn = K(tn)
unp1[0]=0 #Velocity zero ver here unp1 = un
unp1[1:-1]=un[1:-1]+dt*Kn+l*(un[0:-2]-2*un[1:-1]+un[2:]) unp1[0] = 0 # Velocity zero ver here
unp1[-1]=unp1[-2] #Approximate 'infinity' bc unp1[1:-1] = un[1:-1]+dt*Kn+l*(un[0:-2]-2*un[1:-1]+un[2:])
unp1[-1] = unp1[-2] # Approximate 'infinity' bc
return unp1 return unp1
un0=zeros(n,float)
t=0 un0 = np.zeros(n, float)
un=un0 t = 0
un = un0
# un.append(un0) # un.append(un0)
# Make the plot # Make the plot
p.ion() # plt.ioff()
linefd, = p.plot(un0,y) f = plt.figure(1, figsize=(12, 6))
linee, = p.plot(un0,y)
p.legend(('Finite difference','Periodic exact'))
p.ylim(0,1)
p.xlim(-1.5,1.5)
p.ylabel('y')
p.xlabel('u')
p.grid('on')
i=0
uold=un
while(True):
t+=dt
uold=un
un=u_np1(uold,t,dt)
if(i%100==0):
linefd.set_xdata(un)
linee.set_xdata(u_ex(t))
p.draw()
# print("Time:",t)
i+=1
linefd, = plt.plot(un0, y)
linee, = plt.plot(un0, y)
plt.legend(('Finite difference', 'Periodic exact'))
plt.ylim(0, 1)
plt.xlim(-1.5, 1.5)
plt.ylabel('Distance from wall [m]')
plt.xlabel('Velocity [m/s]')
plt.grid(True)
plt.show(block=False)
i = 0
uold = un
try:
while(True):
t += dt
uold = un
un = u_np1(uold, t, dt)
if(i % 100 == 0):
linefd.set_xdata(un)
linee.set_xdata(u_ex(t))
f.canvas.draw_idle()
plt.pause(.00001)
if not plt.fignum_exists(1):
break
i += 1
except Exception:
pass
exit(0)

88
parplateflow.py Executable file → Normal file
View File

@ -1,74 +1,66 @@
#!/usr/bin/python #!/usr/bin/python
# Boundary layer flow # Boundary layer flow
from numpy import * from numpy import zeros, linspace, exp, cos, cosh, sqrt, tanh
import time import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('TkAgg')
# from matplotlib.pylab import *
import pylab as p
# import matplotlib.animation as animation # import matplotlib.animation as animation
def K(t): #Forcing function def K(t): # Forcing function
return (1-exp(-0.1*t))*cos(t) return (1-exp(-0.1*t))*cos(t)
s=10 s = 10
#Define domain
n=50 #Number of gridpoints
y=linspace(-1,1,n)
dy=y[1]-y[0] # Define domain
dt=0.0005 n = 50 # Number of gridpoints
y = linspace(-1, 1, n)
l=(dt/(s**2*dy**2)) dy = y[1]-y[0]
dt = 0.0005
hnu=cosh(sqrt(1j)*s*y)/cosh(sqrt(1j)*s) l = (dt/(s**2*dy**2))
fnu=tanh(sqrt(1j)*s)/(sqrt(1j)*s)
hnu = cosh(sqrt(1j)*s*y)/cosh(sqrt(1j)*s)
fnu = tanh(sqrt(1j)*s)/(sqrt(1j)*s)
def u_ex(tn): def u_ex(tn):
return (((1-hnu)/(1-fnu))*exp(1j*(tn))/1j).real return (((1-hnu)/(1-fnu))*exp(1j*(tn))/1j).real
def u_np1(un,tn,dt):
Kn=K(tn) def u_np1(un, tn, dt):
unp1=un Kn = K(tn)
unp1[0]=0 #Velocity zero ver here unp1 = un
unp1[1:-1]=un[1:-1]+dt*Kn+l*(un[0:-2]-2*un[1:-1]+un[2:]) unp1[0] = 0 # Velocity zero ver here
unp1[-1]=0 #Boundary other side unp1[1:-1] = un[1:-1]+dt*Kn+l*(un[0:-2]-2*un[1:-1]+un[2:])
unp1[-1] = 0 # Boundary other side
return unp1 return unp1
un0=zeros(n,float)
t=0 un0 = zeros(n, float)
un=un0 t = 0
un = un0
# un.append(un0) # un.append(un0)
# Make the plot # Make the plot
p.ion() linefd, = plt.plot(un0, y)
linefd, = p.plot(un0,y) linee, = plt.plot(un0, y)
linee, = p.plot(un0,y) plt.legend(('Finite difference', 'Periodic exact'))
p.legend(('Finite difference','Periodic exact')) plt.ylim(-1, 1)
p.ylim(-1,1) plt.xlim(-1.5, 1.5)
p.xlim(-1.5,1.5) plt.ylabel('y')
p.ylabel('y') plt.xlabel('u')
p.xlabel('u') plt.grid('on')
p.grid('on') i = 0
i=0 uold = un
uold=un
while(True): while(True):
t+=dt t += dt
uold=un uold = un
un=u_np1(uold,t,dt) un = u_np1(uold, t, dt)
if(i%100==0): if(i % 100 == 0):
linefd.set_xdata(un) linefd.set_xdata(un)
linee.set_xdata(u_ex(t)) linee.set_xdata(u_ex(t))
p.draw() plt.draw()
# print("Time:",t) # print("Time:",t)
i+=1 i += 1