1D PIC Code for Two-Stream Plasma Instability
Explore the computational fluid dynamics method for solving the two-stream plasma instability using a 1D PIC code. Learn about the setup, methodology, and parameters involved in simulating this phenomenon to understand plasma behavior better.
Download Presentation
Please find below an Image/Link to download the presentation.
The content on the website is provided AS IS for your information and personal use only. It may not be sold, licensed, or shared on other websites without obtaining consent from the author. If you encounter any issues during the download, it is possible that the publisher has removed the file from their server.
You are allowed to download the files provided on this website for personal or commercial use, subject to the condition that they are used lawfully. All files are the property of their respective owners.
The content on the website is provided AS IS for your information and personal use only. It may not be sold, licensed, or shared on other websites without obtaining consent from the author.
E N D
Presentation Transcript
Example: 1D PIC code for Two- Stream Plasma Instability
Setup x Periodic box 0 L Electrons with position ( ) and velocity ( ) x N v Initial electron distribution function: 2 counter-propagating Maxwellian beams of mean speed bv 1 1 n + = + ( )^ 2 / 2 ( )^ 2 / 2 v v v v ( , ) 0 f x v e e b b 2 2 2 Prof. Succi "Computational Fluid Dynamics" 2
Method Solve electron EOM as coupled first order ODEs using RK4 N 2 dx= dv = ( ) v E x i i i i dt dt where d 2 ( ) ( n ) x n x ( ) d x = = 1 ( ) E x 2 dx dx 0 Prof. Succi "Computational Fluid Dynamics" 3
Method Need electric field at every time step. Solve Poisson s equation (using spectral method) (x ) E Calculate the electron number density source term using PIC approach (x ) n Prof. Succi "Computational Fluid Dynamics" 4
DriverPIC.m == 1D PIC code for the Two-Stream Plasma Instability Problem == Parameters L = 100; % domain of solution 0 <= x <= % domain of solution 0 <= x <= L L N = 20000; % number of electrons % number of electrons J = 1000; % number of grid points % number of grid points vb = 3; % beam velocity % beam velocity dt = 0.1; % time % time- -step (in inverse plasma frequencies) step (in inverse plasma frequencies) tmax = 80; % simulation run from t = 0 to t = % simulation run from t = 0 to t = tmax tmax Prof. Succi "Computational Fluid Dynamics" 5
DriverPIC.m Initialize solution t = 0; rng(42); % seed the rand # generator r = L*rand(N,1); % electron positions v = double_maxwellian(N,vb); % electron velocities % seed the rand # generator % electron positions % electron velocities Prof. Succi "Computational Fluid Dynamics" 6
DriverPIC.m Evolve solution while while (t<=tmax) % load % load r,v solution_coeffs = [r; v]; % take a 4th order % take a 4th order Runge k1 = AssembleRHS AssembleRHS(solution_coeffs,L,J); k2 = AssembleRHS AssembleRHS(solution_coeffs + 0.5*dt*k1,L,J); k3 = AssembleRHS AssembleRHS(solution_coeffs + 0.5*dt*k2,L,J); k4 = AssembleRHS AssembleRHS(solution_coeffs + dt*k3,L,J); solution_coeffs = solution_coeffs + dt/6*(k1+2*k2+2*k3+k4); % unload solution coefficients % unload solution coefficients r = solution_coeffs(1:N); v = solution_coeffs(N+1:2*N); % make sure all coordinates are in the range 0 to L % make sure all coordinates are in the range 0 to L r = r + L*(r<0) - L*(r>L); t = t + dt; end end r,v into a single vector into a single vector Runge- -Kutta Kutta timestep timestep Prof. Succi "Computational Fluid Dynamics" 7
AssembleRHS.m AssembleRHS( solution_coeffs, L, J ) r = solution_coeffs(1:N); v = solution_coeffs(N+1:2*N); r = r + L*(r<0) - L*(r>L); % Calculate electron number density % Calculate electron number density ne = GetDensity GetDensity( r, L, J ); % Solve Poisson's equati % Solve Poisson's equation n0 = N/L; rho = ne/n0 - 1; phi = Poisson1D Poisson1D( rho, L ); % Calculate electric field % Calculate electric field E = GetElectric GetElectric( phi, L ); % equations of motion % equations of motion dx = L/J; js = floor(r0/dx)+1; ys = r0/dx - (js-1); js_plus_1 = mod(js,J)+1; Efield = E(js).*(1-ys) + E(js_plus_1); rdot = v; vdot = -Efield; RHS = [rdot; vdot]; end end function function RHS = AssembleRHS Prof. Succi "Computational Fluid Dynamics" 8
GetDensity.m function function n = GetDensity % Evaluate number density n in grid of J cells, length % Evaluate number density n in grid of J cells, length L, from the electron positions r L, from the electron positions r GetDensity( r, L, J ) dx = L/J; js = floor(r/dx)+1; ys = r/dx - (js-1); js_plus_1 = mod(js,J)+1; n = accumarray(js,(1-ys)/dx,[J,1]) + ... accumarray(js_plus_1,ys/dx,[J,1]); ... end end Prof. Succi "Computational Fluid Dynamics" 9
Poisson1D.m function function u = Poisson1D % Solve 1 % Solve 1- -d Poisson equation: d Poisson equation: % % d^u d^u / dx^2 = v for 0 <= x <= L / dx^2 = v for 0 <= x <= L % using spectral method % using spectral method J = length(v); % Fourier transform source term % Fourier transform source term v_tilde = fft(v); % vector of wave numbers % vector of wave numbers k = (2*pi/L)*[0:(J/2-1) (-J/2):(-1)]'; k(1) = 1; % Calculate Fourier transform of u % Calculate Fourier transform of u u_tilde = -v_tilde./k.^2; % Inverse Fourier transform to obtain u % Inverse Fourier transform to obtain u u = real(ifft(u_tilde)); % Specify arbitrary constant by forcing corner u = 0; % Specify arbitrary constant by forcing corner u = 0; u = u - u(1); end end Poisson1D( v, L ) Prof. Succi "Computational Fluid Dynamics" 10
GetElectric.m function function E = GetElectric % Calculate % Calculate electric field from potential electric field from potential J = length(phi); dx = L/J; % E(j) = (phi(j % E(j) = (phi(j- -1) 1) - - phi(j+1)) / (2* E = (circshift(phi,1)-circshift(phi,-1))/(2*dx); end end GetElectric( phi, L ) phi(j+1)) / (2*dx dx) ) Prof. Succi "Computational Fluid Dynamics" 11
Results Prof. Succi "Computational Fluid Dynamics" 12
