Ideal Gas Simulation

In this worksheet a simple 2-dimensional ideal gas simulation is performed. Particles are moving in a 2D container and collide with the walls. Since we are simulating a dilute, ideal gas, we can ignore the collisions between the gas molecules. This is also called a collisionless plasma. The effect of collisions with the walls is to couple the x and y motions, and to define the pressure of the gas in the container. The temperature of the gas is related to the total kinetic energy. Our gas is simplified in one important aspect: the molecules are treated as point particles, and thus are able to store kinetic energy only through their translational motion. In reality, diatomic (or polytaomic) molecules can convert heat into electronic excitations (eV energy scale), vibrations (0.1 eV scale), and rotations (0.01 eV scale). There are quantum manifestations of these phenomena, i.e., abrupt changes in the specific heat as a function of temperature when a storage mechanism becomes unavailable as the temperature T is lowered. It is useful to know the energy associated with temperature, i.e., kT , which at room temperature corresponds to 1/40 eV (0.025 eV).

The Maxwell-Boltzmann distribution

We begin with the Boltzmann distribution and learn how to generate a random sample with a non-uniform distribution. Our simulation cannot be completely free of assumptions, since it does not describe the relaxation process that takes place over a long time scale: much longer than the time it takes for a particle to traverse the box, i.e., the time scale of our simulation. The relaxation to a steady state from any intial state takes place as a result of collisions between the molecules. The result is the Maxwell-Boltzmann distribution which describes how the particles of the ensemble are distributed over kinetic energy (magnitude of velocity) for given temperature. First-year physics texts provide the main principle behind a derivation. Collisions between the gas molecules imply that the individual molecules can have any amount of kinetic energy, since they exchange it all the time amongst each other. The likelihood to observe a given kinetic energy depends on the number of possibilities to reach that state. To observe a very high, or a very low kinetic energy at given temperature is less likely than some value in between. The average kinetic energy is given by the following relation

> avg:=m/2*u^2=d/2*k*T;

Here d is the number of dimensions (d=3 in nature, d=2 for our simulation), k is Boltzmann's constant, T is the temperature of the gas, m the mass of the particles (molecules), and u is the average(!) velocity. The distribution has to have u for its average speed, but note that the average velocity itself (as a vector) vanishes. For d=3 the distribution reads as:

> fB:=(vx,vy,vz,T)->(m/(2*Pi*k*T))^(3/2)*exp(-m*(vx^2+vy^2+vz^2)/(2*k*T));

Note how the symmetry of the distribution will render a zero result for an expectation value of either of the components ; it is more convenient to define a function that depends on the magnitude of .

> fB:=(v,T)->(m/(2*Pi*k*T))^(3/2)*exp(-m*v^2/(2*k*T));

> int(fB(v,T)*v^2,v=0..infinity);

Definite integration: Can't determine if the integral is convergent.

Need to know the sign of --> 1/2*m/k/T

Will now try indefinite integration and then take limits.

> assume(m>0): assume(k>0): assume(T>0): interface(showassumed=0):

We perform now the calculation of the average in 1D, 2D, 3D: the volume integration in velocity space contributes a factor 1, , and respectively.

> int(fB(v,T)*v^2,v=0..infinity)/int(fB(v,T),v=0..infinity);

> int(fB(v,T)*v^3,v=0..infinity)/int(fB(v,T)*v,v=0..infinity);

> int(fB(v,T)*v^4,v=0..infinity)/int(fB(v,T)*v^2,v=0..infinity);

Now we need to make a choice of units. It is possible (albeit cumbersome) to always carry the units in Maple (multiplied onto the numerical values for the constants). One can define a natural system of units that renders the equations dimensionless. Alternatively one can fix a popular system of units (cgs or MKSA=SI), ignore the units and supply them at the end of the calculation. This is less safe than either of the previous methods, and requires the diligence to supply all variables and constants in the same unit system. We choose the MKSA system of units, for k we omit the J/K (J = Nm = Ws). The Boltzmann constant can be determined from the observation of Brownian motion (Millikan experiment to determine the unit of charge), and is given approximately as:

> k:=1.38*10^(-23);

The mass of the gas particles in kg is assumed as a multiple of the proton mass:

> m:=2*1.67*10^(-27);

We are ready to graph M-B distributions for a few temperatures. We include the volume element for a 2D gas (factor of v ) to observe the probability to find particles in the range [ v , v +d v ]. Before we graph we need to figure out ranges for x and y :

> evalf(fB(v,100));

> plot([10^8*v*fB(v,100),10^8*v*fB(v,300),10^8*v*fB(v,1000)],v=0..0.5*10^4,color=[red,blue,green]);

Note the increase in the most likely speed (in m/sec) as the temperature changes from 100 to 1000 degrees Kelvin. One can see that the distributions are not properly normalized yet, the area changes as the temperature is increased. For this reason the expectation value of was normalized by dividing through the area under the distribution curve. Note that there is a difference between the position of the maximum, and the average value (We leave it as a small exercise to determine both values for a few given temperatures).

Simulations in Maple using random numbers

We proceed with random number generation. Maple has a built-in generator of random integers which can be easily modified to produce equidistributed random numbers on the unit interval [0,1].

> # info can be obtained from ?rand

> randomize(Seed); # this can be changed to generate another sequence

> randx:=()->evalf(rand()/999999999999);

> randx();

A way to check the quality of uniform random numbers is to generate a sequence of points, to display them and observe whether patterns appear.

> for i from 1 to 500 do:

> xv[i]:=randx(): yv[i]:=randx(): od:

> i:='i': plot([[xv[i],yv[i]]$i=1..500],0..1,0..1,style=point, axes=boxed);

Ideally we should find a uniformly covered square. To assess the quality of the generator one has to observe how the square is covered as the number of randoms is increased. The development and testing of pseudorandom number generators is a science (or an art) in its own right. Note that we are using seeded reproducible (deterministic) sequences (thus, pseudorandom numbers) with certain statistical properties. More quantitative ways to test these generators involve the observation of the mean (does it approach 0.5 with a 1/sqrt(N) deviation behaviour), and more generally the quality of the histogram over the unity interval.

> histo:=proc(xv,N,xi,xf,dx) local bin,n,i; global nm;

> nm:=trunc((xf-xi)/dx); for n from 1 to nm do: bin[n]:=0; od:

> for i from 1 to N do:

> n:=trunc((xv[i]-xi)/dx)+1;

> bin[n]:=bin[n]+1; od:

> bin;

> end:

> res:=histo(xv,500,0,1,0.1);

> seq(res[i],i=1..nm);

It is not appropriate to graph these numbers directly. A bar graph is carried below. The op([,]) construct is used to place two neighboring points beside each other in the sequence construct, which is then turned into a list for graphing.

> i:='i': plot([seq(op([[i-1,res[i]],[i,res[i]]]),i=1..nm)]);

Exercise 1:

Observe the deviations from the expected average value (50 for a sample of 500) for each bin as a function of N . How does it decrease with N ?

We can turn the histogram generation into a single procedure (to which we supply the maximum count value as well), and try it out on the other chosen sample of random numbers:

> histo:=proc(xv,N,xi,xf,dx,ymax) local bin,n,i,nm;

> nm:=trunc((xf-xi)/dx); for n from 1 to nm do: bin[n]:=0; od:

> for i from 1 to N do:

> n:=trunc((xv[i]-xi)/dx)+1;

> bin[n]:=bin[n]+1; od:

> i:='i': plot([seq(op([[i-1,bin[i]],[i,bin[i]]]),i=1..nm)],0..nm,0..ymax);

> end:

> histo(yv,500,0.,1.,0.1,80);

To dial up an ensemble that is not equidistributed but follows a given distribution we use a rejection method. First we truncate the infinite range of velocities to a finite one with a cut-off appropriate for the given temperature (see the graph). For finite ensembles this constitutes a small error. Alternatively we could remap the independent variable v , using e.g. the tangent function, such that the infinite range would be mapped to a finite range (0..pi/2).

For T<1000 K we should be able to truncate the velocities at 10,000 m/s. We dial up a velocity in this range and call it v_i. Next we dial a random number in the range 0..1, and calculate the distribution function normalized such that at the maximum of fB a value of 1 is obtained. Now we compare and : we accept v_i as a member of the ensemble only if .

> T0:=1000;

> f0:=evalf(v*fB(v,T0));

> v0:=solve(diff(f0,v)=0,v);

> f0max:=evalf(subs(v=abs(v0[1]),f0));

> f0n:=f0/f0max;

> vmax:=10000;

> Np:=100;

The rejection method requires an infinite loop. For practical reasons we truncate that loop at 1000 attempts.

> #begin loop:

> for i from 1 to Np do:

> vi:=vmax*randx();

> f0i:=subs(v=vi,f0n);

> xi:=randx();

> for j from 1 to 1000 while xi>f0i do: #1000 attempts

> vi:=vmax*randx();

> f0i:=subs(v=vi,f0n);

> xi:=randx(); od:

> vB[i]:=vi;

> od:

> histo(vB,Np,0,vmax,vmax/20,Np/4);

Exercise 2:

Modify the procedure histo to display the velocity instead of bin number on the x axis, and compare the graph with the distribution. The maximum in the fourth bin appears to be in the correct place, i.e., near 2000 m/s for T=1000 K.

In our next step we should dial up velocity vectors. A randomly dialed polar angle in the x - y plane should allow us to define x and y components of the initial velocities. Theta varies between -pi and pi, and the uniform distribution is obtained by stretching the random numbers obtained on a unit interval.

> for i from 1 to Np do:

> theta:=evalf((randx()-0.5)*2*Pi);

> vx[i]:=vB[i]*cos(theta); vy[i]:=vB[i]*sin(theta); od:

We determine some averages for our discrete ensemble.

> vxavg:=0: vyavg:=0: vavg:=0:

> for i from 1 to Np do:

> vxavg:=vxavg+vx[i]: vyavg:=vyavg+vy[i]: vavg:=vavg+sqrt(vx[i]^2+vy[i]^2): od:

> vxavg:=vxavg/Np: vyavg:=vyavg/Np: vavg:=vavg/Np:

> print(`vx,vy,v: `,vxavg,vyavg,vavg);

Note that the components of the average velocity vector should go to zero as the number of test particles is increased, but that the magnitude of the average velocity vector should be given by the ideal gas estimate:

> evalf(sqrt(vxavg^2+vyavg^2),3);

> evalf(sqrt(2*k*T0/m),5); evalf(vavg,5);

A proper histogram of the velocity components can be obtained by using Maple's statistics package. On the x -axis we will find velocities (in m/sec), while the y -axis does not show counts for the bins, but the ratio of counts to total particle number (count frequency).

> with(stats):
with(stats[statplots]):

> histogram(convert(vx,list),color=cyan,area=1);

Simulation of the particle motion

The container is defined through its dimensions: - L < x < L , and - L < y < L . To propagate the ensemble of molecules we can use the solution to Newton's equation for rectilinear motion with possible reflections from one of the four walls.

W choose a size for the container in meters and pick random initial positions.

> L:=0.1;

> for i from 1 to Np do:

> x[i]:=evalf((randx()-0.5)*2*L); y[i]:=evalf((randx()-0.5)*2*L); od:

> i:='i': plot([[x[i],y[i]]$i=1..Np],-L..L,-L..L,style=point,axes=boxed);

To propagate the particles we need to determine the time at which they collide with one of the walls. A scan has to be performed and the particles are updated at each time that one of them hits a wall.

We can limit the run time of the simulation by fixing the total number of collisions with the walls (our choice for ntmax corresponds to 20 wall encounters per particle).

An average collision interval or time scale for our problem is given by:

> tscale:=2*L/vavg;

Inside the loop we determine the time tcoll at which the next particle will encounter one of the walls. This check is done independently for the x - and y -components of each particle's motion. A compact expression can be found that is valid irrespective of the particle direction: for v >0 one relates v t _c with L - x , for v <0 one relates v t _c with - L - x . This can be combined into one expression using L sign( v ).

All other particle positions are updated, and the velocity component perpendicular to the wall is reversed for the colliding particle. In addition we sum the momentum transfer to the container (mass m times twice the velocity component perpendicular to the wall - this comes from elastic collisions of objects of mass m with an infinitely heavy object). The force impulse due to the collision can be converted to an incremental pressure increase on the container (by dividing the force by the container area, and averaging over time). The variable Delv accumulates the momentum transfers (apart from the factor m ). The pressure information is accumulated in P .

> tcoll:=tscale: P[0]:=0:

> tim:=0: Delv:=0: Area:=8*L: Vol:=(2*L)^2: ntmax:=Np*20:

> for ntim from 1 to ntmax do: #step through collisions

> for i from 1 to Np do: #determine shortest collision time

> tx[i]:=(L*sign(vx[i])-x[i])/vx[i]:

> if tx[i]<tcoll then icoll:=i: tcoll:=tx[i]: fi:

> ty[i]:=(L*sign(vy[i])-y[i])/vy[i]:if ty[i]<tcoll then icoll:=i: tcoll:=ty[i]: fi: od:

> tim:=tim+tcoll: #now update positions:

> for i from 1 to Np do:

> x[i]:=evalf(x[i]+vx[i]*tcoll): y[i]:=evalf(y[i]+vy[i]*tcoll): od: tcoll:=tscale:

> #now reverse the momentum (velocity) of colliding TP:

> if tx[icoll]<ty[icoll] then vx[icoll]:=-vx[icoll]: Delv:=Delv+2.*abs(vx[icoll]): else vy[icoll]:=-vy[icoll]: Delv:=Delv+2.*abs(vy[icoll]): fi:

> P[ntim]:=m*Delv/(Area*tim):

> od:

> P[ntmax-1000],P[ntmax-500],P[ntmax];

Evidently the value of the pressure is stable in the later parts of the simulation. To display the time evolution we perform averages over a certain number of steps.

> Pa:=[seq(add(P[(i-1)*20+nu],nu=1..20)/20,i=1..ntmax/20)]:

> i:='i': plot([[i,10^(17)*Pa[i]]$i=2..nops(Pa)],style=point);

Now that we have confidence that the measurement of the gas pressure on the container walls is reasonably stable we can compare it with the expected value from the ideal gas law:

> Pr:=Np/Vol*k*T0;

We are getting a close answer from the simulation. It is possible that the truncation of the velocity interval results in a distribution that corresponds to a somewhat lower temperature. In fact, the average speed in the sample fell short of the expected value.

One should not be too surprised that the simulation comes up with a value for the pressure that agrees with the ideal gas law. The derivation of the ideal gas law is given in first-year physics texts which emphasize the connection between descriptive thermodynamics and kinetic theory (e.g., Serway, 4th ed., chapter 21). The relationship between the gas pressure and the average molecular kinetic energy (which is connected to the temperature) follows from a straightforward argument about the momentum transfer for a particle collision at the wall, and the time it takes to traverse the container between two molecule-wall collisions.

One has to sum over all molecules (thus, the dependence on the particle number, and to convert the force impulse to a pressure.

The simple programming exercise allows one to verify the ideal gas law by simulation. Note that we have not included collisions between the gas molecules (hence ideal gas), and that we have not discussed the mechanism by which the gas acquires and maintains its temperature. In a real setting the heating is provided through (some of) the walls. Maintaining the walls at some temperature implies that the walls are not stationary and impart kinetic energy on the gas molecules during the collisions. Nevertheless, we can verify a few predictions made by the ideal gas law.

Exercise 3:

Repeat the present simulation with a new ensemble of particles. How consistent is the result for the pressure? Increase/decrease the ensemble size by a factor of two, and observe the effect on the scatter in the results for the gas pressure for the chosen temperature value.

Exercise 4:

Repeat the simulation for different values of the temperature T0. Graph the result for the pressure as a function of T0, and compare with the expected result from the ideal gas law. Explain why discrepancies will appear at high T0, unless the algorithm for the sample selection is modified.

Exercise 5:

Change the volume of the container. Repeat simulations at temperature values T0 for which results for the pressure were obtained before. Check the consistency with the ideal gas law.

Exercise 6:

At constant volume and temperature observe the dependence of the pressure on the particle number.

Pressure balance

An interesting extension of the simulation is to replace the top wall of the container by a piston which experiences the force of gravity due to some mass M . The educational program package Albert (Springer-Verlag) contains a simulation of this set-up. This

What modifications do we require from our program to accomodate this change? We can imagine that the height of the piston is a direct function of the gas pressure: if the external pressure exceeds the gas pressure, the piston moves down; if the gas pressure exceeds the pressure due to the mass M pushing the piston down as a result of gravity, the piston moves up. Eventually the piston reaches a height at which the two pressures are equal (equilibrium). It is probably useful to add some friction to the equation of motion for the piston.

The main addition to the simulation is the equation of motion for the piston. This can be handled by a very naive step-by-step integration of Newton's equation for piston height (variable H0 in the program below) and piston speed ( V0 ). In addition one has to modify the collision time criterion by checking separately for collisions with the variable upper wall (piston height H instead of fixed height L ).

The Newton equation for the piston height depends on the pressure of the gas times the surface area of the piston (length L in the 2D simulation). When the external mass on the piston is increased, and the piston moves down, the gas pressure increases on account of the fact that more collisions occur in the vertical direction as a result of the reduced path length. At some point the gas pressure is sufficient to overcome the downward pressure, and some equilibrium is found.

We start the simulation at the end of a previous cycle (piston height L ). We have a well-defined value for the gas pressure.

When choosing the values of parameters we have to keep a few things in mind:

1) this is a 2D simulation, and not a real-life one

2) our sample of gas molecules is extremely small (100 instead of Avogadro's number); it might be better to think of the particles as being bunches of particles, i.e., to increase their mass appropriately. Alternatively, we can adjust the external mass to be very small [of the order of 100/N_A, where N_A=6*10^(23)]

3) the time scale over which simulations need to be run to reach equilibrium may be very long. To adjust for this we chose a strong value for the gravitational acceleration (increase by a factor of 100), and adjusted also the friction constant correspondingly.

The simulation below runs for a very long time in Maple (hours). Simulations are better performed in a compiled floating-point environment (Fortran, C, or Matlab). Nevertheless, one has to make compromises, and limit the realism of a simulation due to lack of computing power.

> M:=20*10^(-21);

> b0:=100.; # friction constant for the piston

> g:=9.81*100; #artificial increase of gravity

> Prt:=Pr;

> H0:=L; V0:=0;

> tcoll:=tscale: P[0]:=Prt:

> ntmax:=Np*50:

> for ntim from 1 to ntmax do: tim:=0: Delv:=0: Area:=4*L+4*H0: for iit from 1 to 40 do: #step through collisions

> for i from 1 to Np do: #determine shortest collision time

> tx[i]:=(L*sign(vx[i])-x[i])/vx[i]:

> if tx[i]<tcoll then icoll:=i: tcoll:=tx[i]: fi:

> if vy[i]>0 then ty[i]:=(H0-y[i])/vy[i]; else ty[i]:=(-L-y[i])/vy[i]: fi: if ty[i]<tcoll then icoll:=i: tcoll:=ty[i]: fi: od:

> tim:=tim+tcoll: #now update positions:

> for i from 1 to Np do:

> x[i]:=evalf(x[i]+vx[i]*tcoll): y[i]:=evalf(y[i]+vy[i]*tcoll): od: tcoll:=tscale:

> #now reverse the momentum (velocity) of colliding TP:

> if tx[icoll]<ty[icoll] then vx[icoll]:=-vx[icoll]: Delv:=Delv+2.*abs(vx[icoll]): else vy[icoll]:=-vy[icoll]: Delv:=Delv+2.*abs(vy[icoll]): fi:

> od: P[ntim]:=m*Delv/(Area*tim):

> H0:=evalf(H0+tim*V0); V0:=evalf(V0+tim*(-g+P[ntim]*L/M-b0*V0)); Hv[ntim]:=H0;

> od:

> P[2],P[3],P[4],P[5],P[6],P[40],P[100];

> seq(Hv[i],i=1..20);

>

> i:='i': plot([[i,10^(17)*P[i]]$i=2..ntmax],style=point,view=[0..ntmax,0..40]);

> i:='i': plot([[i,Hv[i]]$i=1..ntmax],style=point,view=[0..ntmax,-0.1..0.1]);

Please not that the horizontal axis is not labeled by real time, but by the number of collisions. To look at the evolution as a function of time one would have to store the time in a table or list (by accumulating the variable tim in the loop), and plot the pressure and displacement of the piston against this variable.

The friction constant was chosen to be sufficiently large to suppress much large-scale oscillations. It is interesting to observe that despite the erratic input from the pressure measurements the piston height turns out to be a rather smooth function. This is the result of the smoothing properties of the integration of Newton's equations.

Exercise 7:

If your computer runs fast enough (or you can leave it to run unattended for hours) look at the results for different values of the input parameters: (external mass M , friction constant b ). You can also look at runs for different initial temperatures (in this case the first part of the worksheet has to be re-run with a different value of T0 ). For the present choice of T0=1000 (Kelvin) we found the approximate equilibrium positions and pressures:

Mass Height H Pressure P

5*10^(-21) 0.06

1*10^(-20) 0.0 1*10^(-16)

2*10^(-20) -0.04 2*10^(-16)

Verify that pressure time volume (area of the rectangle in the 2D case) remains approximately constant.

>