# Homework 6¶

due 10/31/2018 at 11:59pm

## Problem 1¶

When a fluid is incompressible, $\rho(x,y,z,t)=\rho_0$. Using the conservation of mass and momentum equations with no source of particles show that

1. $\vec\nabla\cdot\vec u=0$
2. Also show that $$\rho_0\frac{\partial\vec u}{\partial t}+\rho_0\big(\vec u\cdot \nabla\big)\vec u=-\vec\nabla p$$

## Problem 2¶

1. Using a version of Debye_length_check from HW04, write a new function that evolves particles in time. In this function, the electric field of all particles is computed on a square mesh. Then the function uses the closest mesh points of any particles to compute the local electric field felt by the particle, using some sort of interpolation (e.g. bilinear). Once the field at each particle location is known, the function "pushes" all the particles in time using a first order scheme.
2. Using this function, plot the time average electric field for a plasma where $n_i=n_e=10^{13}cm^{-3}$, $Z_i=1$ and $T_e=0.2eV$.

NB: If you can make it work, then you just wrote a particle in cell (or PIC) code!

## Problem 1¶

Since the fluid is incompressible $\rho(x,y,z,t)=\rho_0$ and the conservation of mass gives $\vec\nabla\cdot\vec u=0$. From the conservation of momentum, $$\partial_t\rho\vec u+\vec\nabla\cdot(\rho\vec u\vec u)=-\vec\nabla p$$ we have $$\rho_0\partial_t\vec u+\rho_0\vec\nabla\cdot(\vec u\vec u)=-\vec\nabla p$$ Since $\vec\nabla\cdot(\vec u\vec u)|_x=\partial_x u_xu_x+\partial_y u_xu_y+\partial_z u_xu_z=u_x\vec\nabla\cdot\vec u+u_x\partial_xu_x+u_y\partial_yu_x+u_z\partial_zu_x$. Again because $\vec\nabla\cdot\vec u=0$ this is simply $\vec\nabla\cdot(\vec u\vec u)|_x=u_x\partial_xu_x+u_y\partial_yu_x+u_z\partial_zu_x$. So we get $$\rho_0\frac{\partial\vec u}{\partial t}+\rho_0\big(\vec u\cdot \nabla\big)\vec u=-\vec\nabla p$$

## Problem 2¶

We import our modules

In [3]:
import math
from math import sqrt as sqrt
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
%matplotlib inline


We define the particle class

In [4]:
class charged_particle:
m=9.10938356e-31 #mass
q=-1.60217662e-19 #charge
x=np.zeros(3) #location
v=np.zeros(3) #speed


We compute $\vec E$ using the function from HW04

In [5]:
def comp_E(rs,q,rc,r_cutoff=0):
R=rc-rs
r=np.linalg.norm(R)
if (r>r_cutoff):
E=q*R/(4*math.pi*8.85e-12*r**3)
else:
E=np.zeros(3) #avoids infinite electric fields
return E;


This is the same function than the one used in HW04 with some modifications. We replace E_average with E_on_grid that is zeroed out at every time step. We sum all the field from all the particles. Then we compute the force on each particles using charged_particles[j].v+=F*dt/charged_particles[j].m. Note that we use a very coarse interpolation to get the field to and from the grid.

In [6]:
def Debye_length_check(charged_particles,box,time_frame,n_steps,r_cutoff,box_E,N_E):
N=len(charged_particles)
loc=np.zeros((3,n_steps,N))
dt=(time_frame[1]-time_frame[0])/n_steps
lattice_x=np.linspace(box_E[0],box_E[1],N_E)
lattice_y=np.linspace(box_E[2],box_E[3],N_E)
rds=.05
E_ave=np.zeros((3,N_E,N_E))
for i in range(n_steps):
E_on_grid=np.zeros((3,N_E,N_E))
for j in range(N):
for m in range(N_E):
for n in range(N_E):
x1=np.array([lattice_x[m],lattice_y[n],0])
E_on_grid[:,m,n]+=comp_E(charged_particles[j].x,charged_particles[j].q,x1,r_cutoff)
for j in range(N):
k=int((charged_particles[j].x[0]-box_E[0])/(box_E[1]-box_E[0])*(N_E-1))
l=int((charged_particles[j].x[1]-box_E[2])/(box_E[3]-box_E[2])*(N_E-1))
k=max(0,k)
k=min(N_E-1,k)
l=max(0,l)
l=min(N_E-1,l)
F=E_on_grid[:,k,l]*charged_particles[j].q
charged_particles[j].v+=F*dt/charged_particles[j].m
charged_particles[j].x+=charged_particles[j].v*dt # we compute all the velocities then we compute the locations
if (charged_particles[j].x[0]<box[0]): #so randomness in reflections at the wall is not really required
charged_particles[j].x[0]=box[0]-(charged_particles[j].x[0]-box[0])#but it helps the numerics
rd=0.5-np.random.random_sample()/2
charged_particles[j].v[0]*=-1*math.cos(math.pi/4*(1+rds*rd))#note that this way to assign random
charged_particles[j].v[1]*=1*math.sin(math.pi/4*(1+rds*rd))#reflections angles conserve momentum
if (charged_particles[j].x[0]>box[1]):
charged_particles[j].x[0]=box[1]-(charged_particles[j].x[0]-box[1])
rd=rd=np.random.random_sample()
charged_particles[j].v[0]*=-1*math.cos(math.pi/4*(1+rds*rd))
charged_particles[j].v[1]*=1*math.sin(math.pi/4*(1+rds*rd))
if (charged_particles[j].x[1]<box[2]):
charged_particles[j].x[1]=box[2]-(charged_particles[j].x[1]-box[2])
rd=rd=np.random.random_sample()
charged_particles[j].v[0]*=1*math.cos(math.pi/4*(1+rds*rd))
charged_particles[j].v[1]*=-1*math.sin(math.pi/4*(1+rds*rd))
if (charged_particles[j].x[1]>box[3]):
charged_particles[j].x[1]=box[3]-(charged_particles[j].x[1]-box[3])
rd=rd=np.random.random_sample()
charged_particles[j].v[0]*=1*math.cos(math.pi/4*(1+rds*rd))
charged_particles[j].v[1]*=-1*math.sin(math.pi/4*(1+rds*rd))
loc[:,i,j]=charged_particles[j].x
E_ave+=E_on_grid*dt
E_ave/=time_frame[1]-time_frame[0]
return loc,E_ave


Our beloved plot_2D function is defined here.

In [7]:
def plot_2D(r,n_interpolated=1,plot_size=3,plot_axes='axes'):
n = len(r[0,:,0])
k = n_interpolated
fig, ax = plt.subplots(1, 1, figsize=(plot_size, plot_size))
# Now, we draw our points with a gradient of colors.
for i in range(len(r[0,0,:])):
x_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, r[0,:,i])
y_interpolated = np.interp(np.arange(n * k), np.arange(n) * k, r[1,:,i])
ax.scatter(x_interpolated, y_interpolated, c=range(n*k), linewidths=0,
marker='o', s=2*plot_size, cmap=plt.cm.jet,)
#compute autoscale parameters
xc=(r[0,:,:].max()+r[0,:,:].min())/2.
x_low=xc-(r[0,:,:].max()-r[0,:,:].min())/2.*1.1
x_high=xc+(r[0,:,:].max()-r[0,:,:].min())/2.*1.1
yc=(r[1,:,:].max()+r[1,:,:].min())/2.
y_low=yc-(r[1,:,:].max()-r[1,:,:].min())/2.*1.1
y_high=yc+(r[1,:,:].max()-r[1,:,:].min())/2.*1.1
#set autoscale parameters
ax.set_xlim(min(x_low,y_low),max(x_high,y_high))
ax.set_ylim(min(x_low,y_low),max(x_high,y_high))
if (plot_axes!="axes"):
ax.set_axis_off()
else:
plt.xlabel("x (m)")
plt.ylabel("y (m)")


Our initial parameters, similar to HW04.

In [8]:
T_e=.2
n_e=1e19
m_e=9e-31
e=1.6e-19
v_e=sqrt(2*e*T_e/m_e)
omega_e=sqrt(n_e*e**2/(8.85e-12*m_e))
lambda_D=v_e/omega_e
length=8*lambda_D
dt=length/v_e/4
T=2*3.1416/omega_e
N=int(n_e**(1/3)*length)
time_frame=[0,2*T]
n_steps=int((time_frame[1]-time_frame[0])/dt)
print('v_e=',v_e,'m/s')
print('lambda_D=',lambda_D,'m')
print('T=',T,'s')
print('number of particles in the plane:',N,'x',N)
print('time step',dt)
print('number of cycles',n_steps)

v_e= 266666.6666666666 m/s
lambda_D= 1.4874474780643516e-06 m
T= 3.5047237478152256e-11 s
number of particles in the plane: 25 x 25
time step 1.1155856085482638e-11
number of cycles 6


Initialization of the particles and the grid. Nothing different from HW04 except foe the fact that now the grid for E lies outside of the box where the particles are trapped in.

In [9]:
delta_x=length/N/10
delta_y=length/N/10
box=np.array([-length,length,-length,length])
box_E=np.copy(box)+np.array([length/N,-length/N,length/N,-length/N])
box*=.75
lattice_x=np.linspace(-length,length,N)
lattice_y=np.linspace(-length,length,N)
N_E=2*N
cp=np.empty(0,dtype=object)
l=0
for i in range (N):
for j in range (N):
cp=np.append(cp,charged_particle())
cp[l].x=np.array([lattice_x[i],lattice_y[j],0])
cp[l].m=1.67e-27
cp[l].q*=-1
cp[l].v=np.array([0.,0,0])
l+=1

for i in range (N):
for j in range (N):
cp=np.append(cp,charged_particle())
rd=np.random.random_sample()
vx=v_e/sqrt(2)*math.cos(2*math.pi*rd)
vy=v_e/sqrt(2)*math.sin(2*math.pi*rd)
cp[l].x=np.array([lattice_x[i]+delta_x,lattice_y[j]+delta_y,0])
cp[l].v=np.array([vx,vy,0])
l+=1


Let's run it.

In [10]:
r,Ev=Debye_length_check(cp,box,time_frame,n_steps,lambda_D/10,box_E,N_E)


We plot the trajectory

In [11]:
plot_2D(r,n_interpolated=10,plot_size=7)


And we also plot the electric field norm at this time.

In [12]:
fig, ax = plt.subplots(figsize=(10, 10))

im = plt.imshow(np.flip(np.flip(np.hypot(Ev[0,:,:],Ev[1,:,:]),0),1), cmap=plt.cm.jet, extent=box)
im.set_interpolation('bilinear')

cb = fig.colorbar(im)
plt.xlabel('x', rotation=0)
plt.ylabel('y', rotation=90)
cb.ax.set_ylabel('E', rotation=-90)
cb.ax.yaxis.set_label_coords(6, 0.5)
plt.show()