标签：insert phi dh PIC nx ny 源码 matlab np

一、简介

基于matlab PIC模型的太阳风粒子模拟

二、源代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.854e-12;		%permittivity of free space
QE = 1.602e-19;			%elementary charge
K = 1.381e-23;			%boltzmann constant
AMU = 1.661e-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:nx-2,:) - phi(3:nx,:);  %central difference on internal nodes
    efy(:,2:ny-1) = phi(:,1:ny-2) - 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

三、运行结果

四、备注

版本：2014a
完整代码或代写加1564658423

标签：insert,phi,dh,PIC,nx,ny,源码,matlab,np
来源： https://www.cnblogs.com/homeofmatlab/p/14944005.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。