dsmcne.m

% dsmcne - Program to simulate a dilute gas using DSMC algorithm
% This version simulates planar Couette flow
clear;  help dsmcne;   % Clear memory and print header
% Initialize variables  (MKS units)
rand('seed',1);        % Initialize rand
npart = input('Enter number of simulation particles - ');
boltz = 1.3806e-23;    % Boltzmann's constant
mass = 6.63e-26;       % Mass of argon atom
diam = 3.66e-10;       % Effective diameter of argon atom
T = 273;               % Initial temperature
density = 1.78;        % Mass density of argon at STP
L = 1e-6;              % System size is one micron
xwall = [0 L];         % Positions of walls
Volume = L^3;          % Volume of the system
eff_num = (density/mass)*Volume/npart;
fprintf('Each simulation particle represents %g atoms\n',eff_num);
mfp = 1/(sqrt(2)*pi*diam^2*density/mass);
fprintf('System width is %g mean free paths \n',L/mfp);
mpv = sqrt(2*boltz*T/mass);  % Most probable initial velocity 
vwall_m = input('Enter wall velocity as Mach number - ');
vwall = vwall_m * sqrt(5/3 * boltz*T/mass) * [-1 1];
fprintf('Wall velocities are %g and %g m/s \n',vwall(1),vwall(2));
x = L*rand(npart,1);   % Assign random positions
% Assign thermal velocities using Gaussian random numbers
v = mpv * sqrt(-log(ones(npart,3)-rand(npart,3))) ...
                          .* sin(2*pi*rand(npart,3));
% Add velocity gradient
v(:,2) = v(:,2) + vwall(1)+(vwall(2)-vwall(1))*(x(:)/L);
ncell = 20;            % Number of cells
tau = 0.2*(L/ncell)/mpv;   % Set timestep tau
vrmax = 3*mpv*ones(ncell,1);  % Estimated max rel. speed in a cell
selxtra = zeros(ncell,1);     % Used by collision routine "hitter"
coeff = 0.5*eff_num*pi*diam^2*tau/(Volume/ncell);
ave_n = zeros(ncell,1);    % Average cell number density
ave_u = zeros(ncell,3);    % Average cell velocity
ave_T = zeros(ncell,1);    % Average cell temperature
dvtot = zeros(1,2);        % Total momentum change at a wall
dverr = zeros(1,2);        % Used to find error in dvtot
nsamp = 0;                 % Total number of samples
tsamp = 0;                 % Sample time
nstep = input('Enter total number of timesteps - ');
for istep = 1:nstep
  [x, v, strikes, delv] = mover(x,v,npart,xwall,mpv,vwall,tau);
  [cell_n, index, Xref] = sorter(x,npart,ncell,L);
  [v, vrmax, selxtra, col] = colide(v,cell_n,index,Xref,ncell, ...
                                          vrmax,tau,selxtra,coeff);
  fprintf('Finished %g of %g steps, Collisions = %g\n', ...
                                            istep,nstep,col);
  fprintf('Strikes on left wall = %g, right wall = %g \n', ...
                                            strikes(1),strikes(2));
  if(istep > nstep/10) % Start taking samples
    [ave_n,ave_u,ave_T,nsamp]=sampler(x,v,npart,ncell,L, ...
                                      ave_n,ave_u,ave_T,nsamp);
    dvtot = dvtot + delv;
    dverr = dverr + delv.^2;
    tsamp = tsamp + tau;
  end
end
% Normalize statistics
ave_n = ave_n/nsamp;    
for i=1:3
  ave_u(:,i) = ave_u(:,i)/nsamp;
end
ave_T = mass/(3*boltz) * (ave_T/nsamp);
dverr = dverr/(nsamp-1) - (dvtot/nsamp).^2; % Normalize
dverr = sqrt(dverr*nsamp);
xcell(:) = ((1:ncell)-0.5)/ncell * L;
plot(xcell,ave_n); xlabel('position');  ylabel('Number density');
pause;
plot(xcell,ave_u); xlabel('position');  ylabel('Velocities');
pause;
plot(xcell,ave_T); xlabel('position');  ylabel('Temperature');
force = (eff_num*dvtot*mass)/tsamp /L^2;
ferr = (eff_num*dverr*mass)/tsamp /L^2;
fprintf('Force per unit area is \n');
fprintf('Left wall:   %g +/- %g \n',force(1),ferr(1));  
fprintf('Right wall:  %g +/- %g \n',force(2),ferr(2));
vgrad = (vwall(2)-vwall(1))/L;  % Velocity gradient
visc = 1/2*(-force(1)+force(2))/vgrad;  % Average viscosity
viscerr = 1/2*(ferr(1)+ferr(2))/vgrad;  % Error
fprintf('Viscosity = %g +/- %g\n',visc,viscerr);
eta = 5*pi/32*density*(2/sqrt(pi)*mpv)*mfp;
fprintf('Theoretical value of viscoisty is %g \n',eta);