charge1d.m

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% charge1d.m - 1d charging calculation of static dissipative or dielectric materials
% 
% AUTHOR: Alexander C Barrie
%
% LICENSING: This file is intended for educational use. It is licensed
% under GPL (GNU Public License)
%
% This code calculates the charging of a semiconductive strip of material
% measuring 10cm by 1cm with a 3mm thickness. There is an input flux at one
% end and a ground at the other end. This case simulates carbon infused materials
% common to space applications where a dielectric material would typically be used,
% but some level of charge dissipation is required - ESD materials, teflon tubing, etc.
% 
% USE: [dq pot]=charge1d(perm, cond, V0, flux, time, ploton)
% The code is intended to be used in two ways: Setting an initial potential
% with little or no flux will show you how long it takes for a surface
% charge to dissipate to the ground. Having a high flux with little or no
% initial voltage will show you how high that flux will charge the surface
% to. Of course, you can have a high flux and high initial charge to make
% cool movies as well. 
% 
% INPUTS:
% perm - permittivity of material
% cond - conductivity of material
% V0 - initial voltage to put on end of strip
% flux - input flux (#/sq m) on end of strip
% tol - tolerance of qin to qout, which should be equal at steady state
% ploton - 1=turn plotting on, 0=turn plotting off (plotting as in make movie)
% 
% OUTPUTS:
% dq - charge flow between elements. At steady state this will be the same at the 
%      two ends of the strip and ~ 0 in between
% pot - potential on each element of the strip
% mov - an array of frames that you can play with the movie command. This
%       will just be null if ploton is turned off. 
%
% NOTES:
% - The sim is much faster with plotting turned off
% - Modelling conductive materials is impractical due to small timestamp
%   so keep conductivity below ~ 1e-7
% - The number of cells can be changed near the top of the code if need be
%
% EXAMPLE: 
% [dq pot frames]=charge1d(1e-9, 1e-9, 0, 1e14, 1e-3, 0);
% Converged to 0.001000 tolerance with a leading edge potential of 1067.314345V after 3284.356691s and after 18814 iterations 
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [dq pot frames]=charge1d(perm, cond, V0, flux, tol, ploton)

cells=50;

pot=zeros(1,cells+1);                   %set initial potentials to 0

dx=0.1 / cells;                         %square dimensions of each cell (m)
dy=0.01;
dz=0.003;                               %cell thickness (m)

C=perm*dx*dy/dz;                         %parallel plate based capacitance of each cell

pot(1:round(cells/10))=V0;

eq=1.602e-19;                           %electron charge

%zero out some initial stuff
t=0;
i=zeros(1,length(pot));
dq=zeros(1,length(pot));
dt=zeros(1,length(pot));
dt(end)=1e5;

%set up figure
if ploton==1
    frm=0;
    hFig = figure();
    set(hFig, 'Position', [0 500 1000 100]);

    axis tight
    set(hFig,'nextplot','replacechildren');  
end

x=0:dx:cells*dx;
x=x*100;                                %cm

plotcounter=0;

it=0;

converged=0;

% main loop, finish when specified time is reached
while converged==0
    it=it+1;
    i(:)=0;
    q=C*pot;                            %get initial charges
    % loop through each element
    for cnt=1:length(pot)-1             
        E=-(pot(cnt)-pot(cnt+1))/dx;    %E field
        j=E*cond;                       %current density
        
        % now get the current - negative applied to this cell, positive
        % applied to next cell
        i(cnt)=i(cnt)+j*dy*dz;          
        i(cnt+1)=i(cnt+1)-j*dy*dz;      
        
        % calculate the timestep based on the amount of charge we want to
        % move
        if i(cnt) < 0
            % we only want to move a small amount of charge
            % so we check for 1% movement. If you make this
            % too big, you can overshoot.
            dt(cnt)=abs(0.01*q(cnt)/i(cnt));
        else
            %moving charge TO an element, we don't care so much (set it
            %high)
            dt(cnt)=1;
        end
        
    end
    
    % for the overall timestep, use the smallest value
    dtmin=min(dt);
    
    % now calculate the new charge and potential values
    dq=i*dtmin;
    dqin=flux*dx*dy*dtmin*eq;           %q in due to flux
    dq(1)=dq(1)+dqin;
    dpot=dq/C;
    pot=pot+dpot;
    pot(end)=0;
   
    t=t+dtmin;
    
    %calculate convergence based on current in vs out
    conv=abs(dqin/dq(end) -1 );
    
    % make the new movie frame every 100 timesteps
    if mod(plotcounter,25)==0
        if ploton==1
            frm=frm+1;
            contourf(x,[0 1],[pot;pot], 10);
            set(gca, 'ytick',[]);
            set(gca, 'clim', [0 1000]);
            lbl=sprintf('X (cm), t=%1.6fs',t);
            xlabel(lbl);
            h=colorbar;
            colormap(jet(11));
            set(get(h,'ylabel'),'String', 'Potential (V)');
            frames(frm)=getframe(hFig);
            
        end
        fprintf('time: %fs, dt: %f, conv: %f\n', t, dtmin, conv );
        
    end
    plotcounter=plotcounter+1;

    % check for steady state based on the input and output current
    if conv < tol
        converged=1;
    end
    
    if ploton==0
        frames=0;
    end
    
end

fprintf('Converged to %f tolerance with a leading edge potential of %fV after %fs and after %i iterations\n', conv, pot(1), t, it);


%now get analytic answer to compare
R=0.1/(cond*dz*dy);
I=flux*dx*dy*eq;
V=I*R;

fprintf('Analytical solution yields potential of %fV\n', V);