Thermo/Stat Mech (comments to: moloney@nextwork.rose-hulman.edu)  [*.ms = maple, *.xls = excel]

(combin1.ms)

combin1.ms 6/15/96 mjm (in packets\thermo) rev 3

This calculation goes toward showing that the number of ways of sharing n indistinguishable energy units among N distinguishable oscillators is

(n+N-1)!/[ n!(N-1)! ]

Case of N oscillators sharing 1 energy unit:

Looking toward the general method, we could separate off one oscillator and divide the problem into subproblems where the 'separate' oscillator has 0 or 1 energy units. When it has 1 unit of energy, there is only one way the rest can share 0 units. When it has 0 units there are N-1 ways the remaining oscillators can share the 1 unit of energy. Adding these sub-problems together we get 1 + N-1 = N = ways of sharing 1 energy unit among N oscillators.

Case of N oscillators sharing 2 energy units:

For 2 energy units, we have 3 sub-problems; the 'separate' oscillator has 2, 1, or 0 units of energy. When the separate oscillator has 2 units of energy, there is only 1 way for the rest to share 0 units. When the separate oscillator has 1 unit of energy, the other N-1 have to share 1 unit. When the separate oscillator has 0 units, the other N-1 have to share 2 units.

We can write an equation for N oscillators sharing 2 energy units. It will contain 3 terms, one for the separate oscillator having 0, 1, or 2 energy units. We can also write an equation also for N-1 oscillators sharing 2 units of energy, also containing 3 terms. We can keep on writing 3-term equations until we reach the case of 1 oscillator 'sharing' 2 units of energy. This 'last' equation will say there is only one way for 1 oscillator to share 2 units (or any other number of units) of energy.

Adding these equations together gives the number of ways of N oscillators to share 2 energy units. This number of ways equals the sum from 1 to N-1 of sharing 0 energy units, plus the sum from 1 to N-1 of sharing 1 energy unit, plus the number of ways 1 oscillator shares 2 energy units (only 1 way).

Case of N oscillators sharing 3, 4, etc. energy units.

The general method builds up from sharing 0 to sharing 1 to sharing 2, etc, via sums which run over 1 to N-1 of the ways of sharing smaller numbers of energy units, adding in 1 for the way 1 oscillator 'shares' n energy units.

restart; xs:=1: # ways that N share 0 units .
x1:=1+sum(xs,N=1..N-1); # ways that N-1 share 1 unit, +( 1 sharing 1 unit)
xs:=xs+x1; # number of ways of N sharing (0 units & 1 unit)
x1:=factor(1+sum(xs,N=1..N-1)); # sum-> nr of ways of N sharing 2 units

Apparently the sum goes thru ok because N-1 is not put in until the end, and is still symbolic.

xs:=factor(xs+x1); # nr of ways of N sharing (0 & 1 & 2 units)
x1:=factor(1+sum(xs, N=1..N-1)); # sum ->nr of ways of N sharing 3 units
xs:=factor(xs+x1); # nr of ways of sharing (0 & 1 & 2 & 3 units)

The loop below corresponds to a general, iterative method of finding the ways N oscillators can share n energy units. It is not fully general since n is specified at the beginning, but any small n value chosen will be seen to agree with the general formula.


xs:=1: # number of microstates with N oscillators sharing 0 unit of energy
n:=4: # number of energy units to be shared between N oscillators
for k from 1 to n do > x1:=1+sum(xs, N=1..N-1): # add up states for N-1 oscillators
xs:=xs+x1: # running sum of ways, needed for next loop iteration
od:
ans:=factor(x1); # sharing n units of energy between N oscillators

The BE formula tells the # of ways of putting n indistinguishable 'things' into N distinguishable 'boxes'. If the 'things' are oscillator quanta, and the 'boxes' are oscillators, then this corresponds to sharing n units of energy between N distinguishable oscillators.

f:=(n+N-1)!/((n)!*(N-1)!); # B-E formula for n indisting. units in N boxes
ratio:=ans/f;
simplify(ratio); # the ratio should be 1.0 for any selected value of n.

Download combin1.ms [Use'Save File']

Two different kinds of oscillators are used (one has larger energy steps than the other) in equal numbers (N of each) to model degrees of freedom 'freezing out'. The number of energy units is fixed, and the maximum of entropy is determined. At 'low' temperatures the oscillators with large steps have only a small fraction of the total energy. At 'high' temperatures, the large-step oscillators have close to half the energy (approaching equipartition of energy).

comb2osc.ms 6/16/96 mjm sharing energy among 2 oscillator types rev b

We want to model different degrees of freedom, and will do this by having two types of oscillators, N of each, but with different step sizes. Each type of oscillator is supposed to represent a different 'degree of freedom' in a set of atoms. The intent is to see where the entropy is a maximum when the two groups of oscillators have to share a fixed amount of energy. This is supposed to get at the idea of 'freezing out' degrees of freedom at low temperatures. The temperature obtained is 'reduced' in that it represents kT/step , where step represents the size of a step for a given oscillator.

The 'large' step oscillators will have an energy step L times that of the smaller- step oscillators. Nbig will be the number of larger energy units to be shared between the two groups of oscillators.

Given Nbig and N, we will find a maximum in the total system entropy. We will see what fraction of the Nbig energy units are in each group of oscillators.

restart; with(plots):

w:=(N+n-1)!/(n!*(N-1)!); ways to share n units among N oscillators
First, stirling's approx for ln(g!) [ter haar, p.444], then S = ln(w), (k=1).


ws:=ln(sqrt(2*Pi*g)*(g/E)^g);
w1:=expand(ws);
S:=subs(g=N+n-1,w1)-subs(g=n,w1)-subs(g=N-1,w1);

First check heat capacity; get T = 1/dS/dn; Then dT = T(n+1,N)-T(n,N) C = 1/dT (from dQ = TdS = C dT )

T:=simplify(1/diff(S,n));
dT:=subs(n=n+1,T)-T;
c:=1/dT;
Nosc:=100;
plot(subs(N=Nosc,c),n=2..100,title=`C vs energy (1st group)`);
plot(subs(N=Nosc,T),n=2..100,title=`t vs energy (1st group)`);

When Nosc=200, these graphs agree well with the spreadsheet version done by PRM. > Nbig:=50: # number of 'big' energy units shared among 2N oscillators > L:=4: # 2nd group steps are L times larger than 1st group steps

L is the multiplier of large energy steps over small energy steps. Do the two groups of N oscillators. 2nd group energy step = L* 1st group energy step . m is the number of big energy units residing in the 2nd group (large-step oscillators). 2nd group contains mL units of energy, 1st group contains ntotal-mL units.

2nd group has m large units, of L each. Get entropy of each group. Then plot.

s1:=subs(n=L*(Nbig-m),S); s2:=subs(n=m,S);
plot(subs(N=Nosc,s1)+subs(N=Nosc,s2),m=0.05*Nbig..0.85*Nbig,title=`S(both groups) vs big units in 2nd group`);

h:=diff(subs(N=Nosc,s1)+subs(N=Nosc,s2),m); # need derivative to solve for max in S
m1:=fsolve(h=0,m,0.05*Nbig..0.95*Nbig);
t1:=evalf(subs(N=Nosc,n=4*(Nbig-m1),T)); # temp of group 1
t2:=4*evalf(subs(N=Nosc,n=m1,T)); # temp of group 2; step2=4*step1!

Sharing 200 energy units ( 50 big units at 4 apiece) among 100 of 2 types of oscillator gives m=13.8 or so for max entropy. This is not close to equipartition, since equipartition would mean around 25 of the 50 big units.

Increase Nbig, raising the average energy per oscillator (temperature) and then find a larger fraction of energy units among the 'big' oscillators, closer to half.

Or increase N leaving Nbig the same and predict what will happen to m.

Download comb2osc.ms [Use'Save File']

(degfree4.ms)

The heat capacity of metals in the early 1800's was known to lie near 3R per mole for a lot of different metals. This the rule of Dulong and Petit. The gas constant R is boltzmann's constant k multiplied by avogadro's number. Another way of saying it is that R is the gas constant per mole, and k is the gas constant per atom. The 'Einstein solid' is a theoretical model from around 1905 which treated a solid as a collection of atoms each of which was a harmonic oscillator which could vibrate in 3 different directions. The model for N atoms had 3N harmonic oscillators (one for each of the directions). Each oscillator had equally-spaced energy levels of size e. This model gave the heat capacity per atom as

C = 3 k x^2 exp(x)/ [exp(x) -1]^2, where x = e/(kT).

The calculations below 'normalize' the heat capacity so that its maximum value is one (by omitting the factor of 3 k on the front). This makes is the heat capacity per k for a solid with only 1-D vibrations instead of 3-D vibrations. At high temperatures, x approaches zero and C approaches 1.

According to the theorem of equipartition of energy, there should be kT/2 of energy for each degree of freedom available. (A degree of freedom is represented by a squared term in the energy.) A 1-D oscillator has two degrees of freedom, one due to potential energy and the other due to kinetic energy, so each oscillator at high temperatures should have kT of energy on the average. The heat capacity is the rate of change of heat capacity with temperature, so at high temperature, the heat capacity per atom in 1-D vibration should equal boltzmann's constant k. We have divided out the k, so our heat capacity will be 1 at high temperature for each 1-D oscillator. What you will see below is that a temperature can be 'high' for one type of oscillator, and yet be 'low' for another type of oscillator.

e1 is the energy in eV of one set of vibrational levels (one type of oscillator) e2 is the energy in eV of the other set of vibrational levels (second type of oscillator)

For reference, kT = 0.0259 eV at 300 K. Thus, a vibrational degree of freedom which has an energy step of 3 eV should definitely be 'frozen out' at room temperature. This means the temperature is 'low' for that oscillator (or that kT << e for that oscillator.)

The user should try to select a value of e2 so that this degree of freedom is freezing out at around 150 K. (What will be the value of heat capacity in this version of the model if this degree of freedom is 'halfway' frozen out at 150K?)

restart; with(plots):
k:=1.3805e-23/1.602e-19; # k/q in eV/K
e1:=1.0e-4; e2:=3; # energy levels in eV
c:=x^2*exp(x)*(exp(x)-1)^(-2); # Einstein solid heat cap as a fxn of T
plot(subs(x=e1/(k*T),c)+subs(x=e2/(k*T),c),T=10..300,numpoints=300);
plot(subs(x=e1/(k*T),c)+subs(x=0.03/(k*T),c),T=10..300,numpoints=300);

Download degfree4.ms [Use'Save File']

(chempot.ms)

The chemical potential mu has to do with changing numbers of particles. At constant entropy, dU = mu dN, where N is the number of particles in the system. This problem explores the chemical potential by allowing the number of particles to vary, and the system energy to vary such that the entropy stays constant.

chmptlb.ms 6/19/96 mjm rev 4

BACKGROUND.

Consider a group of harmonic oscillators, equally spaced levels, ground state energy of zero. The number of ways of sharing n energy units among N such oscillators is

W = (n+N-1)!/[ n! (N-1)! ]

The entropy is S = k ln W . We will be considering two kinds of oscillators, 'large' and 'small' . The unit of energy for 'small' oscillators will be 'e', and for large oscillators, it will be L times larger, or L e. A small oscillator having n units of energy will have U = internal energy = n e, while a large oscillator having n 'large' energy units will have internal energy U = n L e .

If particle number doesn't change and energy units n can change, then at constant volume and particle number we have TdS = dU. For small oscillators dU = e dn, and for large oscillators dU = L e dn. The temperature is then T = 1/[ dS/dU ] . For small oscillators, T = 1/ [ (k/e) d(ln W )/dn ]. [N is held constant.].

When particle number changes, we have TdS = dU - mu dN. The idea is to set dS=0 and find an expression for dN. We write dS in terms of partials (shown here as ordinary derivatives)

dS = (dS/dn) dn + (dS/dN) dN = 0. (1)

From TdS = dU - mu dN, dS = 1/T [ e dn - mu dN ] = 0 (2)

Solving (2) for mu gives e dn/dN, and from (1) this becomes

mu = -e dS/dN / dS/dn. (3)

We deal with 'reduced' parameters in this program. S is written without the factor of k, the temperatures are written without the factor of e/k, and mu is written without the factor of e. In effect, the reduced temperature is kT/e, and the reduced mu is mu/e.

We will take N1=100 small oscillators and N2 = 100 large oscillators. The total energy available to be shared between these is 30 large units or 30L small energy units.

There are 4 groups of particles, in effect. Initially, we imagine N1 'small' oscillators which have unit energy spacing, and N2 'large' oscillators whose energy spacing is L times larger than the small oscillators. Since we are investigating the chemical potential mu, we let some of the N1 particles wander over into the group of N2 large oscillators, and some of the N2 large oscillators wander over among the N1 small oscillators and let the whole thing come to equilibrium, with each group of oscillators sharing some amount of the total available energy.

The number of 'small' oscillators which have wandered is N1w. The number of (small) energy units in this group is n1w. The large oscillators share a total number of large energy units equal to n2, and the group of N2w large oscillators which wandered share n2w large energy units.

The system entropy is a function of five parameters - N1w, N2w, n2, n1w, and n2w. We obtain the conditions for an entropy maximum with respect to these parameters, and find that all four subsystems have the same temperature when the entropy is a maximum. The groups of small oscillators have the same chemical potential, and the groups of large oscillators have the same chemical potential (different from that of the small oscillators). The total amount of energy shared among the large oscillators is much smaller than that of the small oscillators at equilibrium.

Note on Maple : be sure to use 'assign' statement so solutions are properly assigned ! ! !

restart;with(plots):
L:=5: # factor by which 2nd group energy steps larger than 1st
N1:=100: N2:= 100: Nbig:=30:
ws:=ln(sqrt(2*Pi*g)*(g/E)^g); # stirling's approx for ln(g!) [ter haar, p.444]
w1:=expand(ws);
S:=subs(g=N+n-1,w1)-subs(g=n,w1)-subs(g=N-1,w1);

Ssmlf:=subs(N=N1-N1w,n=L*(Nbig-n2)-n1w,S): # entropy of smalls on the left
Ssmrt:=subs(N=N1w,n=n1w,S): # entropy of smalls on the right (wandered)
Slgrt:=subs(N=N2-N2w,n=n2-n2w,S): # entropy of larges on the right
Slglf:=subs(N=N2w,n=n2w,S): # entropy of larges on the left (wandered)
Stotal:=Ssmrt+Ssmlf+Slgrt+Slglf: # total system entropy
eq1:= diff(Stotal,N1w) = 0, diff(Stotal,N2w) = 0, diff(Stotal,n2)=0 ;
eq2:=diff(Stotal,n1w) = 0, diff(Stotal,n2w)=0;
eqs:={eq1,eq2}: vars:={N1w,n1w,n2,N2w,n2w};
conds:={N1w=1..N1-1,n1w=1..L*Nbig,n2=1..Nbig,N2w=1..N2,n2w=1..Nbig};
sol:=fsolve(eqs,vars,conds); # solve the equations


assign(sol); # Important ! Turns the 'equalities' above into assignments for N1w, etc.


mu:=-diff(S,N)/diff(S,n); # chemical potential
musmall0:=evalf(subs(N=N1,n=L*(Nbig-n2),mu)); # mu(small), no wandering
mularge0:=evalf(subs(N=N2,n=n2,mu)); # mu(small), no wandering
musmrt:=evalf(subs(N=50,n=n1w,mu)); # chemical potential on left, small
musmlf:=evalf(subs(N=50,n=L*(Nbig-n2)-n1w,mu)); # mu rt, small
mulglf:=evalf(subs(N=50,n=n2w,mu)); # mu, on left, large
mulgrt:=evalf(subs(N=50,n=n2-n2w,mu)); # mu on right, large
T:=1/diff(S,n);
Tsmrt:=evalf(subs(N=50,n=n1w,T));
Tsmlf:=evalf(subs(N=50,n=L*(Nbig-n2)-n1w,T)); # T on rt, small
Tlglf:=evalf(L*subs(N=50,n=n2w,T)); # T on left, large
Tlgrt:=evalf(L*subs(N=50,n=n2-n2w,T)); # T on right, large

Temperatures are the same. Mu and T values are close to those with no wandering.

Download chempot.ms [Use'Save File']

(match.ms) (maple version of statme1.xls)

This program matches the statistical model of N oscillators to the einstein model of a solid, and shows the two models agree very well for the specific heat

BACKGROUND. Consider a group of harmonic oscillators, equally spaced levels, ground state energy of zero. The number of ways of sharing n energy units among N such oscillators is

W = (n+N-1)!/[ n! (N-1)! ]

The entropy is S = k ln W . We will be considering two kinds of oscillators, 'large' and 'small' . The unit of energy for 'small' oscillators will be 'e', and for large oscillators, it will be L times larger, or L e. A small oscillator having n units of energy will have U = internal energy = n e, while a large oscillator having n 'large' energy units will have internal energy U = n L e .

If particle number doesn't change and energy units n can change, then at constant volume and particle number we have TdS = dU. For small oscillators dU = e dn, and for large oscillators dU = L e dn. The temperature is then T = 1/[ dS/dU ] . For small oscillators, T = 1/ [ (k/e) d(ln W )/dn ]. [N is held constant.].

When a unit of heat is added, the internal energy changes at constant volume: dQ = TdS = dU =Cv dT. The heat capacity Cv is then dU/dT. dU will be one unit of energy, and dT will be the difference in temperatures.

When particle number changes, we have TdS = dU - mu dN. The idea is to set dS=0 and find an expression for dN. We write dS in terms of partials (shown here as ordinary derivatives)

(1)    dS = (dS/dn) dn + (dS/dN) dN = 0.

From TdS = dU - mu dN,

(2)    dS = 1/T [ e dn - mu dN ] = 0 . 

Solving (2) for mu gives e dn/dN, and from (1) this becomes

(3)    mu = -e dS/dN / dS/dn.

We deal with 'reduced' parameters in this program. S is written without the factor of k, the temperatures are written without the factor of e/k, and mu is written without the factor of e. In effect, the reduced temperature is kT/e, and the reduced mu is mu/e.

Statistical model idea taken from Thomas Moore via a spreadsheet of Dan Schroeder.

restart;with(plots):
w:=(N+n-1)!/(n!*(N-1)!); # nr of ways for n energy steps among N oscs
ws:=ln(sqrt(2*Pi*g)*(g/E)^g); # stirling's approx for ln(g!) [ter haar, p.444]
w1:=expand(ws);
evalf(ln(5!)/subs(g=5,w1)); # ( check the formula at small g; outstanding ! )
S:=subs(g=n+N-1,w1)-subs(g=n,w1)-subs(g=N-1,w1); # log of nr of ways
T:=simplify(1/diff(S,n)); # deriv of reduced entropy wrt energy
dT:=subs(n=n+1,T)-T;
c:=1/dT;
Nosc:=100;
plot(subs(N=Nosc,c),n=10..100,title=`heat capacity vs energy`);
plot(subs(N=Nosc,T),n=10..100,title=`temperature vs energy`);
Ce:=x^2*exp(x)/(exp(x)-1)^2; # Einstein heat capacity
Ce:=subs(x=1/t,Ce); # put it in terms of reduced temperature.
cvs:=plot( [subs(N=Nosc,T),subs(N=Nosc,c),n=10..100],color=red,title=`Cv (stat. model) vs t`):
display(cvs);

This plot tells the temperature and heat capacity based on n the number of energy units. Use it to guide the temperature range of the next plot coming up.

With 100 oscillators, both models are approaching a heat capacity of 1 per oscillator at n=100 energy units. Would this also be true if we had 500 oscillators and 100 energy units? Make a prediction before you try it!

Set t range to that of last plot.

cve:=plot(Nosc*Ce,t=0.426..1.46,color=blue,title=`Cv (Einstein model) vs t`):
display(cve);
display({cvs,cve}, title=`Cv stat(red) and Einstein(blue)`); # compare the two models

Download match.ms [Use'Save File']

statme1.ms 1/1/97 mjm sharing energy among groups of oscillators

Consider a group of harmonic oscillators, equally spaced levels, ground state energy of zero. The number of ways of sharing n energy units (1 energy unit is the spacing between oscillator energy levels) among N such oscillators is

W = (n+N-1)!/[ n! (N-1)! ]

The entropy is S = k ln W .

We will be sharing a total of n energy units between a group of N1 oscillators, and another group of N2 oscillators. In the implementation below, n=100, N1=300, and N2 = 200. n1 is the number of energy units shared among the group of N1=300 oscillators. We will let n1 vary from 0 to n (n=100), and plot individual entropies s1 and s2 of the subgroups, and the total system entropy s1+s2. We will also plot the total number of ways (the total microstates) by plotting the exponential of S=s1+s2.

This is a reduced system of units, further discussed in xxxx.ms .

From this notion, we can go on to temperature as the derivative of entropy with respect to energy (ds1/dn1, ds2/dn2, etc.).

From there we can examine heat capacity, and chemical potential.

When the groups of oscillators are different, (one group has one spacing, the other group has a different spacing), we can discuss 'freezing out' of degrees of freedom.

restart;with(plots):

# w:=(N+n-1)!/(n!*(N-1)!); ways to share n units among N oscillators

ws:=ln(sqrt(2*Pi*g)*(g/E)^g); # stirling's approx for ln(g!) [ter haar, p.444]

w1:=expand(ws);

S:=subs(g=N+n-1,w1)-subs(g=n,w1)-subs(g=N-1,w1);

s1:=subs(N=300,n=n1,S); # entropy w/ n1 energy units among 300 oscillators

s2:=subs(N=200,n=100-n1,S); # entropy w/ 100-n1 units among 200 oscillators

plot({s1,s2,s1+s2},n1=0..100,title=`Entropy vs. n1`);

plot(exp(s1+s2),n1=0..100,title=`Total microstates vs. n1`);

Download statme1.ms [Use'Save File']

(vander.ms)

VAN DER WAALS EQUATION OF STATE

The equation of state for an ideal gas is Pv = RT, where v is the specific volume (volume per mole). The model on which this equation is based assumes negligible forces between the molecules and a negligible amount of space taken up by the molecules themselves. Both of assumptions are reasonable except when the volume and temperature are low enough that the molecules spend a significant amount of time in the neighborhood of other molecules, where they might experience forces exerted by them. If these forces are attractive it is reasonable to assume that the pressure on the walls of a container would be reduced, since a molecule that is about to collide with a wall will be attracted by its nearest neighbors in the interior and momentarily slowed down. In the van der Waals model the pressure is reduced by an amount a/(v)^2 . Also the van der Waals model assumes that the volume in which molecules can move around is less than the volume of the container because of the space taken up by the molecules themselves, so that the v in the equation of state should be replaced by a smaller volume v-b. The van der Waals equation of state for a (non-ideal) gas then becomes:

(P+a/(v)^2)*(v-b) = R*T

(a) Write down the van der Waals equation and solve for P. Substitute the value of the ideal gas constant R = 8314 Joules per kilomole, and the appropriate constants a and b for the gas.

(b) Choose three temperatures (about 40 degree intervals) for which you think that graphs of P vs v might look non-ideal in appropriate ranges of specific volume, and substitute them into the expression for P.

(c) Plot on the same graph three P vs. v curves for the three temperatures and specific volume ranges from roughly 0.06 to about 0.4 cubic meters per kilomole . Your lowest temperature line should have a temperature such that there is a region of slope zero, but not so low that there is a region of positive slope (positive slope regions occur in the liquid region where the van der Waals equation is not valid).

(d) Repeat for an ideal gas (P=RT/v) and compare the graph with that obtained using van der Waals' equation of state.

Solution:

restart; with(plots):

Part A Write down the van der Waals equation and solve for P.

eq1:=(p+a/v^2)*(v-b)=R*T;
Pressure:=solve(eq1,p);

Substitute the ideal gas constant R = 8314 (J/kilomole K)

R:=8314:

Substitute van der Waals constants for CO2 a = 366000 (J m^3/kilomole) and b = .0429 (m^3/kilomole)

a:=366000: b:=.0429:

---

Substitute van der Waals constants for oxygen (O2) (these will supersede those of CO2) a = 13800 (J m^3/kilomole) and b = .0318 (m^3/kilomole)

a:=13800: b:=.0318: R:=8314:

---

Part B:

Calculate three expressions relating pressure to specific volume for three particular temperatures in the range where deviations from ideal gas behavior might be expected.

Pressure1:=subs(T=300,Pressure);
Pressure2:=subs(T=340,Pressure);
Pressure3:=subs(T=380,Pressure);

Part C:

Plot P vs v for the above three temperatures:

plot({Pressure1,Pressure2,Pressure3},v=0.068..0.4,p=0..2*10^7,title=`Pressure vs Volume`);
plot3d(Pressure,v=0.068..0.4,T=250..380,view=0..2e+007);

Part D: Compare the above graph with that obtained from the ideal gas equation of state for the same temperatures and the same range of specific volumes.

Pr:=R*T/v;
Pr1:=subs(T=300,Pr);
Pr2:=subs(T=340,Pr);
Pr3:=subs(T=380,Pr);
plot({Pr1,Pr2,Pr3},v=.068..0.4);

Download vander.ms [Use'Save File']

(vdwcrit.ms) (Could be assigned as several problems)

VAN DER WAALS EQUATION OF STATE - reduced equation values at the critical point

The van der Waals equation of state for a (non-ideal) gas is:

(P+a/(v)^2)*(v-b) = R*T

a) Find the pressure, temperature, and volume at the critical point. [ dP/dv)T = 0, etc.

b) Obtain the 'universal' van der Waals equation by using 'reduced' variables p/p(critical), v/v(critical) and T/T(critical).

c) Show that the pressure in a van der Waals gas falls to zero when T = 27/32 Tc.

d) Use the given van der Waals parameters a and b for CO2 and O2 to obtain values of the critical temperature, volume and pressure of O2 and CO2. Compare these values to the experimental values and explain how the given van der Waals constants were probably selected.

Solution:

restart; with(plots):
p1:=-a/v^2+ R*T/(v-b);
pdv1:=diff(p1,v)=0; pdv2:=diff(p1,v,v)=0;
sol:=solve({subs(v=vc,T=tc,pdv1),subs(v=vc,T=tc,pdv2)},{vc,tc});
assign(sol);                    # vc, tc = critical volume, temperature
pc:=subs(v=vc,T=tc,p1); # pressure at the critical point
eqvdw:= p = -a/v^2+ R*T/(v-b);
eq3:=(p+3/v^2)*(3*v-1) = 8*t; # reduced equation of state v/vcr, p/pcr, t/tcr
R:=8314:
eq4:=subs(p=0,eq3); # eqn for volume when p=0; 2 real solutions => p went <0
v1:=solve(eq4,v); # this says pressure goes negative if 81>96t. t<81/96=27/32 => p<0.

Substitute oxygen van der Waals constants a=138000, b= 0.0318 O2 (Sears-Salinger) critical constants are tc=154.8 K, vc = 0.0725 (m^3/kmole) pc= 5.02x10^6 N/m^2

a1:=138000: b1:=0.0318:

Substitute carbon dioxide van der Waals constants a=13800, b= 0.0318 CO2 (Sears-Salinger) critical constants are tc=304.2 K, vc = 0.094 (m^3/kmole) pc=7.3x10^6 N/m^2

a2:=366000: b2:=0.0429:
vcrox:=evalf(subs(a=a1,b=b1,vc));
tcrox:=evalf(subs(a=a1,b=b1,tc));     pcrox:=evalf(subs(a=a1,b=b1,pc));
vcrcd:=evalf(subs(a=a2,b=b2,vc));
tcrcd:=evalf(subs(a=a2,b=b2,tc));
pcrcd:=evalf(subs(a=a2,b=b2,pc));

pox:=subs(a=a1,b=b1,p1);
plot(subs(T=155,pox),v=1.5*b1..5*b1);
plot(subs(T=26*155/32,pox)/1000000,v=1.5*b1..5*b1,title=`Yoicks! The pressure went below zero.`);

Download vdwcrit.ms [Use'Save File']

(joulthom.ms)

JOULE-THOMSON COOLING

This problem presents some background on the Joule-Thomson effect, and proceeds to find the Joule-Thomson coefficient for a van der Waals gas. This coefficient enables one to say if a gas will heat or cool upon expansion through a porous plug (or nozzle). The inversion curve for helium is obtained, showing in what range of parameters helium can be cooled upon expansion. (Follows Sommerfeld's text.)

When a stream of gas at pressure p1 and temperature T1 is forced through a throttling valve, it emerges at a lower pressure and a temperature which may be lower or higher, depending on the gas and the initial conditions. For this effect to be observed the gas must be nonideal; there must be forces between the molecules. The process increases the time average distance between molecules. If conditions are such that the average intermolecular force are attractive (rather than repulsive) as the average intermolecular distances increase, then the process increases the potential energy of the system; we might expect at least some of this increase to come at the expense of kinetic energy. A drop in kinetic energy means a lower temperature. However, if conditions are such that average intermolecular forces are repulsive as intermolecular distances increase, potential energy will decrease and we might expect an increase in kinetic energy and temperature. Whether the temperature increases or decreases depends on the intial state of the gas (T1, p1), although in all cases the enthalpy h=u+pv is constant. Generally, if the initial pressure is high (pushing the molecules close together) repulsive forces will dominate and the process will produce a heating effect, but if the initial pressure is lower and the temperature is in the right range there can be a cooling effect. This cooling effect has been used in the process of liquefying gases.

The Joule-Thomson coefficient is the indicator of whether the throttling process produces cooling or heating and how effectively it does so. It is defined as the rate of change of temperature with respect to pressure as the enthalpy is held constant (partial derivative of T w.r.t. p, constant h). From the first and second laws of thermodynamics, the Joule-Thomson coefficient (dT/dp) = - (1/cp)*(dh/dp) = (1/cp)[T(dv/dT)-v]. If this coefficient is positive, reductions in pressure will cause reductions in temperature; this is what we would expect to happen at lower initial pressures (in the proper temperature range). If the coefficient is negative, then reductions in pressure cause increases in temperature; we would expect this to happen at higher initial pressures. Therefore, for a gas which has some particular enthalpy h, decreasing the initial pressure (in the proper temperature range) causes the effect to change from one that causes heating to one that causes cooling at some "inversion" initial pressure and temperature. The "inversion curve" is the locus of these inversion points on a p-T graph..

The following exercise uses the van der Waals model of a gas, which takes account of intermolecular forces by adding a/(v)^2 to the pressure term in the equation of state (see the van der Waals exercise):

(p+a/(v)^2)*(v-b)=R*T

The calculations use several relationships that come from the first and second laws of thermodynamics.

(a) Solve van der Waals' equation for pressure and calculate partial derivatives dp/dT, dp/dv, and dv/dT=-(dp/dT)/(dp/dv). Also find the volume expansion coefficient beta=(1/v)*(dv/dT) and the isothermal compressibility kappa=-(1/v)*(dv/dp). Simplify these expressions as much as possible.

(b) Find an expression for the rate of change of internal energy with respect to changes in volume when the temperature is held constant:

du/dv=T*(dp/dT)-p Simplify this expression.

(c) Find an expression for the Joule-Thomson coefficient (see the discussion above) and simplify it. Set the coefficient equal to zero and solve for T to find an expression for the inversion temperature in terms of volume v. Solve this equation for v and substitute into the van der Waals equation to get an expression relating the inversion temperature to pressure. Solve for the inversion pressure and simplify.

(d) Substitute the value for R (8314 J/kilomole K) and the appropriate constants a (3440 J m^3/kilomole^2) and b (0.0234 m^3/kilomole) for Helium. Plot the inversion curve on a p-T graph for the temperature range from 4 to 35K. For what part of the graph would the throttling process produce a cooling effect. In order to use the throttling process to cool Helium for the purpose of liquefying it, to what temperature must you get it first by other methods?

(e) Repeat part (d) for Oxygen (a = 138000 Jm^3/kilomole^2, b = 0.0318 m^3/kilomole) and find the temperature-pressure range for which the Joule-Thomson effect results in cooling. In the graph, eliminate negative pressures. [There are temperatues and pressures that fit van der Waals equation that are completely "unphysical"; the van der Waals model doesn't fit reality everywhere.]

__________________________________________________________________________________

Part A:

restart;
eq1:=(p+a/(v)^2)*(v-b)=R*T
Press:=solve(eq1,p);
dpdT:=diff(Press,T);
dpdv:=diff(Press,v);
dpdvs:=simplify(dpdv);
dvdT:=-dpdT/dpdv;
dvdTs:=simplify(dvdT);
beta:=(1/v)*(dvdTs);
kappa:=-(1/v)/dpdvs;

Part B:

dudv:=T*dpdT-Press;
dudvs:=simplify(dudv);

Part C:

dpdTh:=(1/cp)*(T*(dvdTs)-v);
pThs:=simplify(dpdTh);
Tinv:=solve(pThs=0,T);
Tinv:=factor(Tinv);
v2:=solve(Tinv=T,v)[1];
sol:=subs(v=v2,eq1);
pinv:=solve(sol,p);
pinv2:=simplify(pinv);

Part D:

pinv3:=subs(R=8314,a=3440,b=.0234,pinv2);
plot(pinv3,T=4..35,title=`Inversion Pressure vs Temperature (He)`);

Part E:

pinv4:=subs(R=8314,a=138000,b=.0318,pinv2);
plot(pinv4,T=115..1045,title=`Inversion Pressure vs. Temperature (O2)`);

Download joulthom.ms [Use'Save File']

(maxwell.ms)

In the Maxwell theory of an ideal gas the fraction of molecules with speeds between v and v+dv is given by f(v)dv = 4*Pi*(m/(2*Pi* k*T))^(3/2)*v^2*exp(-m*v^2/(2*k*T)), where m is the mass of the gas molecule, k is the Boltzmann constant, and T is the absolute (Kelvin) temperature. The most probable speed is the one that that gives a maximum value for f(v). The average speed is obtained by integrating the product vf(v) over the entire range of speeds. Similarly the mean square speed is found by integrating the v^2*f(v) over the entire range; the root mean square is the square root of the mean square speed. The average molecular kinetic energy is (1/2) m (vrms)^2 . (a) For oxygen molecules at 300K, calculate kT/m. Then make a plot of the speed distribution function f(v) versus the speed v for a speed range from 0 to 1200 m/s.

(b) In order to find the most probable speed v, differentiate f(v) with respect to v and set it equal to zero (look for a solution in the range v=1..1200). Compare your result with your graph of part (a) and with (2*k*T/m)^(1/2).

(c) Find the average molecular speed and the rms speed and locate them on your graph. Compare your average v with (8*k*T/(Pi*m))^(1/2) the rms speed with (3*k*T/m)^(1/2) .

(d) Calculate the fraction of the molecules having speeds which are less than the average speed .

(e) Calculate the average kinetic energy and compare with 3kT/2 . ________________________________________________________________________

Solution

restart;

Part A Calculate kT/m

con:=(k*T)/m;
k:=138*10^(-25);
T:=300;
m:=5316*10^(-29);
c:=subs(m=1329/25000000000000000000000000000,k=69/5000000000000000000000000,T=300,con);

Speed Distribution Function

fv:=4*Pi*(1/(2*Pi*c))^(3/2)*v^2*exp(-v^2/(2*c));

Plot f(v) vs v

plot(fv,v=0..1200,title=`Maxwell Speed Distribution Function`);

Part B:

Find the derivative of fv with respect to v.

der:=diff(fv,v);

Set the derivative equal to zero and solve for v.

sol:=fsolve(der=0,v=1..1200);
vm:=(2*c)^(1/2);
evalf(vm);

Part C:

Average speed

vav:= int(v*fv,v=0..infinity);
evalf(vav);
va:=(8*c/Pi)^(1/2);
evalf(va);

Root mean square speed

rms:=(int((v^2)*fv,v=0..infinity))^(1/2);
evalf(rms);
vrm:=(3*c)^(1/2);
evalf(vrm);

Part D:

Find the fraction of the molecules with speeds less than vav:

fra:=int(fv,v=0..vav);
evalf(fra);

Part E:

Find the average kinetic energy

ke:=(1/2)*m*(rms)^2;
evalf(ke);
kin:=3*k*T/2;
evalf(kin);

Download maxwell.ms [Use'Save File']

boltzmann distribution (rotsteps.ms)

Given the distance between H and Br in HBr is 1.4 A = 0.14 nm, (m(H) = 1 amu, and m(Br) = 80 amu, (1 amu = 1.66 x 10^-27 kg)

a) calculate the moment of inertia of HBr with respect to its center of mass.

b) Since for a rotating body the kinetic energy is 1/2 I w^2 = 1/2 L^2/I , and angular momentum L is quantized in steps of h/(2 pi), calculate the size of the smallest rotational 'energy step'.

c) Find the ratio of HBr molecules in the first rotational excited state (one unit of angular momentum) to the number in the ground state (zero units of angular momentum) at 500K.

d) Plot the ratio in part c) from 100K to 600K e) Repeat part d) for the 10th excited rotational state, plotting from 100 K to 500K

restart; with(plots):
l1:=ltot*m2/(m1+m2); # distance from m1 to CM
l2:=ltot*m1/(m1+m2); # ditto for m2
RotInertia:=simplify(m1*l1^2 + m2*l2^2); # moment of inertia about CM
Eunit:=(hbar^2)/(2*RotInertia);
A:=1*10^(-10): # angstrom
M:=1.66*10^(-27): # amu
Boltzfactor:=exp(-n*(n+1)*Eunit/(k*T)); # n is the rotational state quantum number
boltz:=subs(hbar=6.626*10^(-34)/(2*Pi),k=1.3805*10^(-23),m1=80*M,m2=M,ltot=14*A/10,Boltzfactor);
b500:=evalf(subs(T=500,n=1,boltz)); # ratio at 500 K
plot(subs(n=1,boltz),T=100..600,title=`HBr populations: n_1/n_0`); # plot of n_1 / n_0
plot(subs(n=10,boltz),T=100..600,title=`HBr populations: n_10 / n_0`); # plot of n_10 / n_0

Download rotsteps.ms [Use'Save File']

Classic cycle of two isotherms and two adiabats using ideal gas. Easily modified for other kinds of cycles.

cycle1.mws Plot and calculate a thermodynamic cycle. Determine efficiency.

isotherm is now 1->2 and not adiabat! This is correct order of cycles.

restart;with(plots):
adiabat:= pa*Va^(5/3)= pb*Vb^(5/3);
isotherm:= pa*Va = pb*Vb;
idealgas:= pa*Va = 831/100*T;
p1:=2.5*10^5; V1:=0.0226; V2:=0.05; V3:=0.07;
startpV:={pa=p1,Va=V1};
T1:=solve(subs(startpV,idealgas),T);
p2:= solve(subs(pa=p1,Va=V1,Vb=V2,isotherm),pb);
p3:=solve(subs(pa=p2,Va=V2,Vb=V3,adiabat),pb); sol:=solve({subs(pa=p3,Va=V3,isotherm),subs(pa=p1,Va=V1,adiabat)},{pb,Vb});
assign(sol);
p4:=pb; V4:=Vb;
p_adiabat:=p0*(V0/V)^(5/3);
p_isothermal:=p0*V0/V;

plot12:=plot(subs(p0=p1,V0=V1,p_isothermal),V=V1..V2): work12:=int(subs(p0=p1,V0=V1,p_isothermal),V=V1..V2); plot23:=plot(subs(p0=p2,V0=V2,p_adiabat),V=V2..V3):
work23:=int(subs(p0=p2,V0=V2,p_adiabat),V=V2..V3);
display({plot12,plot23});

plot34:=plot(subs(p0=p3,V0=V3,p_isothermal),V=V3..V4): work34:=int(subs(p0=p3,V0=V3,p_isothermal),V=V3..V4);
display({plot12,plot23,plot34});

plot41:=plot(subs(p0=p4,V0=V4,p_adiabat),V=V4..V1):
work41:=int(subs(p0=p4,V0=V4,p_adiabat),V=V4..V1);
display({plot12,plot23,plot34,plot41});
work:=work12+work23+work34+work41;
T2:=solve(subs(pa=p3,Va=V3,idealgas),T);

Maximum efficiency in a Carnot cycle
EffMax :=1-T2/T1;

Efficiency is net work done divided by heat input (from 1->2 at higher temperature)
EffEngine := work/work12;

Download cycle1.mws [Use'Save File']