#!/usr/bin/env python # coding: utf-8 # # 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$ # 1. Solar corona plasma: $n_e=10^{15}m^{-3}$ and $T_e=100eV$ # 2. Solar wind plasma: $n_e=10^{7}m^{-3}$ and $T_e=10eV$ # 3. A fusion reactor: $n_e=10^{23}m^{-3}$ and $T_e=10keV$ # 4. A Z-pinch: $n_e=10^{25}m^{-3}$ and $T_e=500eV$ # 5. 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: # 2. 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 get_ipython().run_line_magic('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') # ## 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) # 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[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[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) # ... 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[ ]: