Homework 4

due 10/10/2018 at 11:59pm

Problem 1

Compute the collisions frequencies $\nu_{ee}$, $\nu_{ei}$ and $\nu_{ii}$ for the following hydrogen plasmas:

  1. Earth ionosphere: $n_e=10^{11}m^{-3}$ and $T_e=10meV$
  2. Solar corona plasma: $n_e=10^{15}m^{-3}$ and $T_e=100eV$
  3. Solar wind plasma: $n_e=10^{7}m^{-3}$ and $T_e=10eV$
  4. A fusion reactor: $n_e=10^{23}m^{-3}$ and $T_e=10keV$
  5. A Z-pinch: $n_e=10^{25}m^{-3}$ and $T_e=500eV$
  6. A laser plasma: $n_e=10^{28}m^{-3}$ and $T_e=500eV$

Note that we use hydrogen, so $n_i=n_e$ and $Z=1$. Where necessary we will suppose that $T_i=T_e$

Problem 2

  1. Compute the Debye length $\lambda_D$ for several densities between $10^{19}m^{-3}$ and $10^{20}m^{-3}$ and for several temperatures between $0.1eV$ and $1eV$.

Using a version of binary_Coulomb_collisions that can handle any arbitrary number of particles, generate a hydrogen plasma inside a 2-dimensional box, using the previous densities and temperatures, supposing that $n_i=n_e$ and that the ions do not move. Make sure that the box is much larger than the large Debye length compute in question 1. After several time steps:

  1. Plot the average electric field inside the box for each density and temperature (Hint: sum the electric fields over time as vectors, then take the strength of this average electric field)
  2. Compute the Debye length $\lambda_D$ using the information you got from this electric field
  3. Compare it to the theoretical values computed in 1.
In [1]:
import math
from math import sqrt as sqrt
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
%matplotlib inline 

Problem 1

This is a pretty straightforward problem. First we define our three functions

In [2]:
def nu_ei(n,T):
    return sqrt(2)*n*1**2*1.6e-19**4*12/(12*math.pi**(3/2)*8.85e-12**2*9e-31**(1/2)*(1.6e-19*T)**(3/2))
def nu_ii(n,T):
    return sqrt(2)*n*1**4*1.6e-19**4*12/(12*math.pi**(3/2)*8.85e-12**2*1.6e-27**(1/2)*(1.6e-19*T)**(3/2))
def nu_ee(n,T):
    return n*1**2*1.6e-19**4*12/(8.85e-12**2*9e-31**(1/2)*(1.6e-19*T)**(3/2))

Note that we could have added all the parameters like $m_e$ or $T_i$ as inputs. We used $\ln\Lambda=10$

In [3]:
property_list=np.array([1e11,1e15,1e7,1e23,1e25,1e28,1e-2,100,10,1e4,500,500])
property_list.shape=(2,6)
answers=np.zeros((3,6))
for i in range (len(answers[0,:])):
    answers[0,i]=nu_ee(property_list[0,i],property_list[1,i])
    answers[1,i]=nu_ei(property_list[0,i],property_list[1,i])
    answers[2,i]=nu_ii(property_list[0,i],property_list[1,i])
for i in range (len(answers[0,:])):
    print('nu_ee for a plasma of density',property_list[0,i],'m^-3 and temperature',property_list[1,i],'V is',answers[0,i],'Hz')
nu_ee for a plasma of density 100000000000.0 m^-3 and temperature 0.01 V is 165376.3515726602 Hz
nu_ee for a plasma of density 1000000000000000.0 m^-3 and temperature 100.0 V is 1653.763515726602 Hz
nu_ee for a plasma of density 10000000.0 m^-3 and temperature 10.0 V is 0.0005229659420983752 Hz
nu_ee for a plasma of density 1e+23 m^-3 and temperature 10000.0 V is 165376351.57266018 Hz
nu_ee for a plasma of density 1e+25 m^-3 and temperature 500.0 V is 1479171055949.49 Hz
nu_ee for a plasma of density 1e+28 m^-3 and temperature 500.0 V is 1479171055949489.8 Hz

Problem 2

We first define our charge particle class.

In [52]:
def Debye_length(n,T):
    omega_pe=sqrt(n*1.6e-19**2/(8.85e-12*9e-31))
    v_Te=sqrt(2*1.6e-19*T/9e-31)
    return v_Te/omega_pe
In [56]:
property_list=np.array([1e19,5e19,1e20,0.1,0.5,1])
property_list.shape=(2,3)
answers=np.zeros((3,3))
for i in range (len(answers[0,:])):
    for j in range (len(answers[0,:])):
        answers[i,j]=Debye_length(property_list[0,i],property_list[1,j])
print(answers)
[[1.05178420e-06 2.35186097e-06 3.32603367e-06]
 [4.70372193e-07 1.05178420e-06 1.48744748e-06]
 [3.32603367e-07 7.43723739e-07 1.05178420e-06]]
In [57]:
class charged_particle:
    m=9.10938356e-31 #charge 
    q=-1.60217662e-19 #mass
    x=np.zeros(3) #location
    v=np.zeros(3) #speed

And our electric field function

In [13]:
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;

We rewrite binary_Coulomb_collisons as Debye_length_check.

In [63]:
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)
    E_average=np.zeros((3,N_E,N_E))
    rds=.05
    for i in range(n_steps):
        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_average[:,m,n]+=comp_E(charged_particles[j].x,charged_particles[j].q,x1,r_cutoff)*dt
        for j in range(N):
            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_average/=time_frame[1]-time_frame[0]    
    return loc,E_average

We have made two modifications here. First, we do not compute the electric of each particles on the other. The main reason is to follow Vlasov and realize that the Debye length is coming from screening of the ions by the moving electrons. Not from Coulomb collisions. While taking into account Coulomb collisions is more physical, the computation of trajectories would require a huge amount of small steps.

What we have done instead is adding random motion every time a particle hits the edge of the box. It is important to note that we change the components of the vector velocity in such a way that the strength of the vector stays the same. So our random transformations converse kinetic energy.

Every time a particle arrives at the edge of the box, we use reflecting boundary conditions (e.g. $v_x$ turns into $-v_x$ for the boundaries located at $x=-length$ or $x=length$. The random motion comes from the fact that there is some slippage on the box walls. For instance, let's look at the two vertical wall for constant $x$: $$\begin{align} &v_{x_{after}}=-v_{x_{before}}\cos\big(\frac{\pi}{4}[1+rd]\big)\\ &v_{y_{after}}=v_{y_{before}}\sin\big(\frac{\pi}{4}[1+rd]\big) \end{align}$$ where $rd$ is a random number between $[-0.5,0.5]$.

This randomness is necessary to approximate "thermal" motion. This process fixed out boundary conditions for both the $x$ and $y$ walls. Note that for the $y$ wall we simply switched the $-$ from the $x$ component to the $y$ component of the velocity.

Then we use our usual plot_2D function.

In [64]:
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)")

Now we define our parameters ...

In [70]:
T_e=.1
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,25*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= 188561.80831641264 m/s
lambda_D= 1.0517841983981314e-06 m
T= 3.5047237478152256e-11 s
number of particles in the plane: 18 x 18
time step 1.1155856085482638e-11
number of cycles 78

... and our initial conditions. Note that the velocity is "thermal" only in randomness direction but not in velocity strength. This is not a problem here.

In [71]:
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*=1.5
lattice_x=np.linspace(-length,length,N)
lattice_y=np.linspace(-length,length,N)
N_E=N-1
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

We put the grid in between the ions to avoid the large E that would results from a grid point falling exactly where an ion is. We now run the code.

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

We plot the electron trajectories. Quite random!

In [73]:
plot_2D(r,n_interpolated=1,plot_size=7)

We plot our average field.

In [74]:
E_val=np.zeros((N_E,N_E))
for i in range (N_E):
    for j in range (N_E):
        E_val[i,j]=np.linalg.norm(E_ave[:,i,j])

fig, ax = plt.subplots()

im = plt.imshow(np.flip(np.flip(E_val,0),1), cmap=plt.cm.jet, extent=box_E)
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()

According to the grid size we picked we should have eight Debye lengths across. We see that the geometrical size of field fluctuations is on the order of the Debye length.

In [ ]: