In [1]:
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as sf
from IPython import display
import math as mt
import matplotlib.animation as animation
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['animation.embed_limit']=60
plt.rcParams['animation.html'] = 'jshtml'
Drift Wave Instability¶
Kono&Miyashita1988Phys.Fluids
Governing Equation¶
The nonlinear equations of the collisional drift wave instability can be obtained as follows,
$$\left[\left(1-\kappa\partial_t\right)\left(\partial_t+\partial_y\right)-\left(\partial_t+\partial_y-\nu\delta\nabla^2\right)\nabla^2\right]\phi+\partial_zv_z=\left[\left({\bf z}\times\nabla_\perp\phi\right)\cdot\nabla_\perp\right]\nabla^2_\perp\phi;$$$$\partial_tv_z+\partial_z\phi=\left[\left({\bf z}\times\nabla_\perp\phi\right)\cdot\nabla_\perp\right]v_z.$$Futher assuming that the ion motion in the z-direction is ignored ($k_{\parallel}\gg k_\perp$) and $\kappa\partial_t$ can be replaced by $\kappa\partial_y$, the two equations reduce to the following equation,
$${\partial _t}\left( {1 - \nabla _ \bot ^2 + \kappa {\partial _y}} \right)\phi + \left[ {{\partial _y} + \kappa \partial _y^2 - \delta \left( {{\partial _y} - \nu \nabla _ \bot ^2} \right)\nabla _ \bot ^2} \right]\phi = \left[ {\left( {{\bf{\hat z}} \times {\nabla _ \bot }\phi } \right) \cdot {\nabla _ \bot }} \right]\nabla _ \bot ^2\phi. $$Linear Dispersion Relation¶
Linearizing the last equation, we obtain the dispersion relation of the collisional drift wave instability for the real and imaginary frequency,
$$\omega=\frac{k_y\left[(1+k^2)(1+\delta k^2)+\kappa(\kappa k_y^2-\nu\delta k^4)\right]}{(1+k^2)^2+\kappa^2k_y^2};$$$$\gamma=\frac{k^2\left[\kappa(1-\delta)k_y^2-\nu\delta k^2(1+k^2)\right]}{(1+k^2)^2+\kappa^2k_y^2}.$$Here is an example of the real (solid) and imaginary frequency (dashed) with parameters, $k_0=\pi/16$, $k_x=4k_0$, $\kappa=\nu=0.3$, and $\delta=0.1$. Note that at $k_y/k_0=8$ the growth rate becomes $\sim0.05$.
In [2]:
k0=np.pi/16
kap=0.3; nu=0.3; delt=0.1
ky=np.linspace(0,16,100)*k0
kx=4*k0
kk=np.sqrt(ky**2+kx**2)
wr=ky*((1+kk**2)*(1+delt*kk**2)+kap*(kap*ky**2-nu*delt*kk**4))
wr=wr/((1+kk**2)**2+kap**2*ky**2)
wi=kk**2*(kap*(1-delt)*ky**2-nu*delt*kk**2*(1+kk**2))
wi=wi/((1+kk**2)**2+kap**2*ky**2)
plt.plot(ky/k0,wr,'k',label='$\omega_r$')
plt.plot(ky/k0,wi*4,'k--',label=r'$4\times\omega_i$')
plt.legend()
plt.xlabel('$k_y/k_0$')
plt.ylabel('$\omega$')
plt.xlim(0,16)
plt.ylim(0,0.6)
plt.show()
Simulation of the Drift Wave Instability¶
In [17]:
def DWI(nx,ny,tmax,dx,dy,kap,nu,delt,phi,tsav):
global KX,KY,KX2,KY2,KXD,KYD
### define grids ###
kx =2*np.pi/(nx*dx)*np.r_[np.arange(nx/2),np.arange(-nx/2,0)]
ky =2*np.pi/(ny*dy)*np.r_[np.arange(ny/2),np.arange(-ny/2,0)]
kxd=np.zeros(nx)
kyd=np.zeros(ny)
kxd[np.abs(kx)<=2/3*np.max(abs(kx))]=1.0 ### 2/3-rule for dealiasing
kyd[np.abs(ky)<=2/3*np.max(abs(ky))]=1.0 ### 2/3-rule for dealiasing
kx2=kx**2; ky2=ky**2
KX ,KY =np.meshgrid(kx ,ky )
KX2,KY2=np.meshgrid(kx2,ky2)
KXD,KYD=np.meshgrid(kxd,kyd)
KD=np.zeros((nx,ny))
for i in range(nx):
for j in range(ny):
if(np.sqrt((i-nx//2)**2+(j-ny//2)**2)>1/3*(nx//2)):
KD[j,i]=1
phif=sf.fft2(phi)
phihst=np.zeros((1,nx,ny))
phihst[0,:,:]=np.real(sf.ifft2(phif))
nt=0
t=0
isav=1
dt=0.005
while t < tmax:
#if(np.max(np.abs(np.real(sf.ifft2(phif))))>1):
# dt=0.05*dx/(np.max(np.abs(np.real(sf.ifft2(phif))))/1)
#else:
# dt=0.05*dx
#---4th-order Runge-Kutta---#
#phif=np.exp(-(1j*KY-kap*KY2+delt*(1j*KY+nu*(KX2+KY2))*(KX2+KY2))/(1+(KX2+KY2)+kap*1j*KY)*dt)*phif
#gw1=adv(phif )
#gw2=adv(phif+0.5*dt*gw1)
#gw3=adv(phif+0.5*dt*gw2)
#gw4=adv(phif+ dt*gw3)
#phif=phif+dt*(gw1+2*gw2+2*gw3+gw4)/6.0
##---4th-order Runge-Kutta---#
#A=-(1j*KY-kap*KY2+delt*(1j*KY+nu*(KX2+KY2))*(KX2+KY2))/(1+(KX2+KY2)+kap*1j*KY)*0.5*dt
#gw1=adv(phif )
#gw2=adv(np.exp(A)*(phif +0.5*dt*gw1))
#gw3=adv(np.exp(A)*(phif)+0.5*dt*gw2)
#gw4=adv(np.exp(A)*(np.exp(A)*phif+dt*gw3))
#phif=np.exp(A)*(np.exp(A)*(phif+dt*gw1/6)+dt*(gw2+gw3)/3)+gw4/6
phif=np.exp(-(1j*KY-kap*KY2+delt*(1j*KY+nu*(KX2+KY2))*(KX2+KY2))/(1+(KX2+KY2)+kap*1j*KY)*dt)*phif
flx=adv(phif)
phiff = phif + dt*flx
flx=adv(phiff)
phiff = 3/4*phif + 1/4*(phiff + dt*flx)
flx=adv(phiff)
phif = 1/3*phif + 2/3*(phiff + dt*flx)
if(t>isav*tsav):
isav=isav+1
phi=np.real(sf.ifft2(phif))
phihst=np.vstack((phihst,[phi]))
#print(nt,t,np.max(phi))
nt=nt+1
t =t +dt
if(np.max(np.real(sf.ifft2(phif)))>1000): break
phi=np.real(sf.ifft2(phif))
phihst=np.vstack((phihst,[phi]))
#print(nt,t)
return locals()
def adv(phif):
phixf = 1j*KX*phif ;phix =np.real(sf.ifft2(phixf *KXD*KYD))
phiyf = 1j*KY*phif ;phiy =np.real(sf.ifft2(phiyf *KXD*KYD))
phi2xf = -1j*KX*(KX2+KY2)*phif;phi2x =np.real(sf.ifft2(phi2xf*KXD*KYD))
phi2yf = -1j*KY*(KX2+KY2)*phif;phi2y =np.real(sf.ifft2(phi2yf*KXD*KYD))
### Nonlinear part ###
nonl=sf.fft2(-phiy*phi2x+phix*phi2y)/(1+(KX2+KY2)+kap*1j*KY)
return nonl
Initial Condition¶
For the initial conditions, we choose $$\phi=\phi_0\sin(k_xx+k_yy), \phi_0=0.1$$ $$(k_x,k_y)=(4k_0,8k_0), k_0=\pi/16.$$ Here, we use $\kappa=\nu=0.3$ and $\delta=0.1$. The grid points and the system size are $N_x\times N_y=128\times 128$ and $L_x\times L_y=32\times32$.
In [18]:
%%time
nx=128; ny=128; tmax=200; tsav=1
kap=0.3; nu=0.3; delt=0.1
lx=32; ly=32
dx=lx/nx; dy=ly/ny
k0=2*np.pi/lx
x =np.arange(nx)*dx
y =np.arange(ny)*dy
X,Y=np.meshgrid(x,y)
phi=0.1*np.sin(4*k0*X+8*k0*Y)
data=DWI(nx,ny,tmax,dx,dy,kap,nu,delt,phi,tsav)
locals().update(data)
Nonlinear Time Evolution of $\phi$¶
In [19]:
def update_anim(it):
fig.clf()
ax1 = fig.add_subplot(111)
ax1.clear()
im1=ax1.imshow(phihst[it,:,:]/np.max(phihst[it,:,:]),extent=[0,nx*dx,0,ny*dy],aspect='auto',origin='lower',cmap='gray')
#im1=ax1.imshow(phihst[it,:,:],extent=[0,nx*dx,0,ny*dy],aspect='auto',origin='lower',cmap='gray')
fig.colorbar(im1, ax=ax1,label='$\phi/\phi_{max,t}$')
ax1.set_title(r'$t=$%d' %(it*tsav))
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
In [20]:
fig=plt.figure(figsize=(6,5))
anim=animation.FuncAnimation(fig,update_anim,frames=np.size(phihst,0))
plt.close()
anim
Out[20]: