A list,
Solar wind particle Simulation based on MATLAB PIC model
Ii. Source code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Si
%
% For more,
% and
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%reset variables
clear variables
%identify globals needed by the potential solver
global EPS0 QE den A n0 phi0 phi_p Te box
%setup constants
EPS0 = 8.854 e-12; %permittivity of free space
QE = 1.602 e-19; %elementary charge
K = 1.381 e-23; %boltzmann constant
AMU = 1.661 e-27; %atomic mass unit
M = 32*AMU; %ion mass (molecular oxygen)
%input settings
n0 = 1e12; %density in #/m^3
phi0 = 0; %reference potential
Te = 1; %electron temperature in eV
Ti = 0.1; %ion velocity in eV
v_drift = 7000; %ion injection velocity, 7km/s
phi_p = - 5; %wall potential
%calculate plasma parameters
lD = sqrt(EPS0*Te/(n0*QE)); %Debye length
vth = sqrt(2*QE*Ti/M); %thermal velocity with Ti in eV
%set simulation domain
nx = 16; %number of nodes in x direction
ny = 10; %number of nodes in y direction
ts = 200; %number of time steps
dh = lD; %cell size
np_insert = (ny- 1) *15; %insert 15 particles per cell
%compute some other values
nn = nx*ny; %total number of nodes
dt = 0.1*dh/v_drift; %time step, at vdrift move 0.10dx
Lx = (nx- 1)*dh; %domain length in x direction
Ly = (ny- 1)*dh; %domain length in y direction
%specify plate dimensions
box(1, :) = [floor(nx/3) floor(nx/3) +2]; %x range
box(2, :) = [1 floor(ny/2)]; %y range
%create an object domain for visualization
object = zeros(nx,ny);
for j=box(2.1):box(2.2)
object(box(1.1):box(1.2),j)=ones(box(1.2)-box(1.1) +1.1);
end
%calculate specific weight
flux = n0*v_drift*Ly; %flux of entering particles
npt = flux*dt; %number of real particles created per timestep
spwt = npt/np_insert; %specific weight, real particles per macroparticle
mp_q = 1; %macroparticle charge
max_part=20000; %buffer size
%allocate particle array
part_x = zeros(max_part,2); %particle positions
part_v = zeros(max_part,2); %particle velocities
%set up multiplication matrix for potential solver
%here we are setting up the Finite Difference stencil
A = zeros(nn); %allocate empty nn * nn matrix
%set regular stencil on internal nodes
for j=2:ny- 1 %only internal nodes
for i=2:nx- 1
u = (j- 1)*nx+i; %unknown (row index)
A(u,u) = 4 -/(dh*dh); %phi(i,j)
A(u,u- 1) =1/(dh*dh); %phi(i- 1,j)
A(u,u+1) =1/(dh*dh); %phi(i+1,j)
A(u,u-nx)=1/(dh*dh); %phi(i,j- 1)
A(u,u+nx)=1/(dh*dh); %phi(i,j+1)
end
end
%neumann boundary on y=0
for i=1:nx
u=i;
A(u,u) = - 1/dh; %phi(i,j)
A(u,u+nx) = 1/dh; %phi(i,j+1)
end
%neumann boundary on y=Ly
for i=1:nx
u=(ny- 1)*nx+i;
A(u,u-nx) = 1/dh; %phi(i,j- 1)
A(u,u) = - 1/dh; %phi(i,j)
end
%neumann boundary on x=Lx
for j=1:ny
u=(j- 1)*nx+nx;
A(u,:)=zeros(1,nn); %clear row
A(u,u- 1) = 1/dh; %phi(i- 1,j)
A(u,u) = - 1/dh; %phi(i,j)
end
%dirichlet boundary on x=0
for j=1:ny
u=(j- 1)*nx+1;
A(u,:)=zeros(1,nn); %clear row
A(u,u) = 1; %phi(i,j)
end
%dirichlet boundary on nodes corresponding to the plate
for j=box(2.1):box(2.2)
for i=box(1.1):box(1.2)
u=(j- 1)*nx+i;
A(u,:)=zeros(1,nn); %clear row
A(u,u)=1; %phi(i,j)
end
end
%initialize
phi = ones(nx,ny)*phi0; %set initial potential to phi0
np = 0; %clear number of particles
disp(['Solving potential for the first time. Please be patient, this could take a while.']);
%%%%%%%%%%%%%%%%%%%%%%%%
% MAIN LOOP
%%%%%%%%%%%%%%%%%%%%%%%%
for it=1:ts %iterate for ts time steps
%reset field quantities
den = zeros(nx,ny); %number density
efx = zeros(nx,ny); %electric field, x-component
efy = zeros(nx,ny); %electric field, y-component
chg = zeros(nx,ny); %charge distribution
%*** 1. CALCULATE CHARGE DENSITY ***
% deposit charge to nodes
for p=1:np %loop over particles
fi = 1+part_x(p,1)/dh; %real i index of particle's cell
i = floor(fi); %integral part
hx = fi-i; %the remainder
fj = 1+part_x(p,2)/dh; %real j index of particle's cell
j = floor(fj); %integral part
hy = fj-j; %the remainder
%interpolate charge to nodes
chg(i,j) = chg(i,j) + (1-hx)*(1-hy);
chg(i+1,j) = chg(i+1,j) + hx*(1-hy);
chg(i,j+1) = chg(i,j+1) + (1-hx)*hy;
chg(i+1,j+1) = chg(i+1,j+1) + hx*hy;
end
%calculate density
den = spwt*mp_q*chg/(dh*dh);
%apply boundaries
den(1, :) = 2*den(1, :); %double density since only half volume contributing
den(nx,:) = 2*den(nx,:);
den(:,1) = 2*den(:,1);
den(:,ny) = 2*den(:,ny);
%add density floor for plotting and to help the solver
den = den + 1e4; % * * *2.CALCULATE POTENTIAL *** phi = eval_2dpot_GS(phi); % * * *3. CALCULATE ELECTRIC FIELD ***
efx(2:nx- 1, :) = phi(1:nx2 -,:) - phi(3:nx,:); %central difference on internal nodes
efy(:,2:ny- 1) = phi(:,1:ny2 -) - phi(:,3:ny); %central difference on internal nodes
efx(1, :) = 2*(phi(1,:) - phi(2:)); %forward difference on x=0
efx(nx,:) = 2*(phi(nx- 1,:) - phi(nx,:)); %backward difference on x=Lx
efy(:,1) = 2*(phi(:,1) - phi(:,2)); %forward difference on y=0
efy(:,ny) = 2*(phi(:,ny- 1) - phi(:,ny)); %forward difference on y=Ly
efx = efx / (2*dh); %divide by dominator
efy = efy / (2*dh); % * * *4. GENERATE NEW PARTICLE ***
if (np+np_insert>=max_part) %make sure we don't exceed array limits
% np_insert=max_part-np;
end
%insert particles randomly distributed in y and in the first cell
part_x(np+1:np+np_insert,1)=rand(np_insert,1)*dh; %x position
part_x(np+1:np+np_insert,2)=rand(np_insert,1)*Ly; %y position
%sample Maxwellian in x and y, add drift velocity in x
part_v(np+1:np+np_insert,1)=v_drift+(1.5+rand(np_insert,1)+rand(np_insert,1)+rand(np_insert,1))*vth;
part_v(np+1:np+np_insert,2) =0.5* (1.5+rand(np_insert,1)+rand(np_insert,1)+rand(np_insert,1))*vth;
np=np+np_insert; %increment particle counter
%*** 5. MOVE PARTICLES ***
p=1;
while(p<=np) %loop over particles
fi = 1+part_x(p)/dh; %i index of particle's cell
i = floor(fi);
hx = fi-i; %fractional x position in cell
fj = 1+part_x(p,2)/dh; %j index of particle' cell
j = floor(fj);
hy = fj-j; %fractional y position in cell
%gather electric field
E=[0 0];
E = [efx(i,j) efy(i,j)]*(1-hx)*(1-hy); %contribution from (i,j)
E = E+ [efx(i+1,j) efy(i+1,j)]*hx*(1-hy); %(i+1,j)
E = E + [efx(i,j+1) efy(i+1,j)]*(1-hx)*hy; %(i,j+1)
E = E + [efx(i+1,j+1) efy(i+1,j+1)]*hx*hy; %(i+1,j+1)
%update velocity and position
F = QE*E; %Lorentz force, F=qE
a = F/M; %acceleration
part_v(p,:) = part_v(p,:)+a*dt; %update velocity
part_x(p,:) = part_x(p,:)+part_v(p,:)*dt; %update position
Copy the code3. Operation results
Fourth, note
Version: 2014a complete code or write plus 1564658423